# An Improved Equilibrium Optimizer Algorithm for Solving Optimal Power Flow Problem with Penetration of Wind and Solar Energy

Authors: Ngoc Anh Nguyen; Dieu Ngoc Vo; Thuan Thanh Nguyen; Thanh Long Duong
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

The paper is proposed an improved equilibrium optimizer (IEO) algorithm to solve the optimal power flow (OPF) problem with the participation of a renewable energy source (RES). In the proposed IEO method, the “exponential term” is replaced by a function that does not dependent on the number of iterations. This modification of the IEO algorithm increases exploration ability compared to EO algorithm. In addition, the exploration of the proposed IEO algorithm will not decrease according to the number of iterations which avoids to get stuck at the local optimal solution. The IEO algorithm is tested on two IEEE 30-bus and IEEE 118-bus systems with three different objective functions. The performance of the proposed IEO method is compared with equilibrium optimizer (EO), artificial ecosystem optimization (AEO), cuckoo search algorithm (CSA), teaching-learning-based optimization (TLBO), artificial bee colony (ABC), and many other existing methods. Besides, a simple probabilistic formula for calculating RES output power based on the Monte-Carlo simulation model is proposed in this paper to reduce the computation time for the OPF problem with RES. The simulation results obtained show that the proposed IEO algorithm has better quality of the solution as well as stability level compared to the original EO algorithm and other algorithm. Thus, the proposed IEO algorithm is also one of effective and reliable algorithms for handling OPF problem with RES.


---

## Body

## 1. Introduction

Optimal power flow (OPF) is an optimizing tool for power system planning and operation. The OPF has been received significant interest of power system optimization research group since first introduced in 1962 [1]. The main aim of the OPF problem is to determine the optimal setting control variables which optimize chosen objective functions such as fuel cost, emission cost, power loss, and voltage profile while satisfying a set of operational and physical constraints. The control variables of the OPF problem are that the actual power at the generator buses, excluding the slack bus, the voltage magnitude at all the generator buses, tap changer of transformer, and shunt compensators. Many mathematical programming methods have been deployed to deal with OPF problems, such as linear programming (LP) [2], nonlinear programming (NLP) [3], Newton-based techniques [4], quadratic programming (QP) [5], and interior-point (IP) methods [6]. However, the objective functions of the OPF problem, which was solved by these conventional methods, are simple and differentiable. In addition, the energy sources of the conventional OPF problem only involved thermal power sources. In fact, the OPF problem in modern power systems is always a nonlinear optimization problem and may be a nondifferentiable one; thus, it is an actual challenge for optimization methods for dealing with, especially, the conventional methods. In order to overcome the limitations of classical methods, heuristic methods have been considered as alternative approaches to solve the OPF problem with the advantages of obtaining nearly optimum solution whether the problem is differentiable or not.In recent decades, heuristic optimization algorithms have been developed based on inspired from physical phenomena, animal behavior, or evolutionary concepts. They provide straightforward and effective solutions for the OPF problem. The typical algorithms such as particle swarm optimization (PSO) [7], gravitational search algorithm (GSA) [8, 9], differential evolutionary (DE) [10], krill herd algorithm (KHA) [11], artificial bee colony (ABC) [12], biogeography-based optimization (BBO) [13], Jaya algorithm [14], harmony search (HS) [15], teaching-learning-based-optimization technique (TLBO) [16], Sine-Cosine (SC) [17], grey wolf optimizer (GWO) [18], and moth swarm approach (MSA) [19]. The objective functions of the OPF problem usually are considered consisting of three different types of fuel functions, namely, quadratic cost, piecewise quadratic cost, quadratic cost curve, and power loss, emission cost and voltage deviation. In addition, the objective of maximum social welfare [20] and available transfer capability [21] are also used when solving the OPF problem by heuristic optimization algorithms. In general, these methods are successfully applied in the OPF problem. However, the solution of these algorithms often falls into the local optimum for the large-scale power system.So, in order to improve the solution quality and efficiency of the metaheuristic methods, the hybrid and improved algorithms have been proposed to solve the OPF problem [22–30]. The goal of the hybrid algorithms is to achieve optimal balance between exploration and exploitation strategy. In [22], a hybrid (SFLA-SA) algorithm of simulated annealing algorithms (SFLAs) with the probabilistic shuffle jump characteristic of the shuffled frog algorithm (SA) is proposed to avoid stagnation at the local solution of SFLA. In [23], the hybrid (PSO-GSA) algorithm is proposed for the OPF problem by combining the social thinking ability in PSO and the local search ability of GSA to reach the better solution quality compared to PSO. In [24], the hybrid algorithm between cuckoo search algorithm (CSA) and sunflower algorithm (SFO) (HCSA-SFO) is proposed, wherein the Lévy flight function of CSA has been replaced by the mutation and selection mechanism of SFO to reach the higher solution quality and shorter executed time over CSA. Many studies have been proposed in recent years with aims of improving exploration performance or the solution quality by using random jumps of Lévy flight or mutation vector to increase diversity of population such as the improved grey wolf optimizer (DGWO) [25], modified honey bee mating optimization (MHBMO) [26], and skip shuffle modification (MSFL) [27]. An improved moth-flame optimization algorithm (IMFA) for solving optimal power flow problem is proposed in [28]. In [29], the improved collision object optimization (ICBO) algorithm is proposed by using three colliding objects instead of two objects collision as the original algorithm. A modified imperialist competitive algorithm (MICA) with a review of the imperial power mechanism to promote global optimization has been proposed in [30] to improve the simulation time as well as the optimal solution. These hybrid or improved algorithms in general have yielded significant performance compared to the original version in terms of obtained solution quality as solving the OPF problems which have many local optimizations.Nowadays, renewable energy sources (RESs) including wind power and solar energy are gradually replacing conventional thermal power plants due to advantages of these energy sources. Integrating wind and solar power on the traditional grid can have significant impact on the efficiency of power system operation. However, they have posed many difficulties in planning and operating the power system due to their uncertain and discontinuous nature. Therefore, solving the optimal power system (OPF) problem to minimize the objective functions in the power system containing traditional power plants and renewable energy power plants is one of the important tasks for the independent operator system.Several published papers related to the OPF problem in the power system including conventional power plants and renewable energy power plants have been presented in recent years. Optimization methods are widely used to solve problems related to penetration of renewable energy. A hybrid Harris hawks optimizer for integration of renewable energy sources considering stochastic behavior of energy sources is proposed in [31]. In [32], authors have been proposed differential evolutionary particle swarm optimization to deal with optimal power flow of power systems with controllable wind-photovoltaic energy systems. A barnacle mating optimizer (BMO) has been proposed to solve the optimal power flow (OPF) problem with and without renewable energy sources [33]. The effectiveness of the proposed BMO in solving the OPF is tested on a modified IEEE-30 bus system that is integrated with solar PV farms. In [34], authors have been proposed an algorithm hybrid based on combination of phasor particle swarm optimization and a gravitational search algorithm (PPSOGSA) for the OPF in power systems with an integrated wind turbine (WT) and solar photovoltaic (PV) generators. The forecasted active power of WT and PV are considered as dependent variables in the OPF formulation, while the voltage magnitude at WT and PV buses is considered as control variables. Forecasting the output power of WT and PV generators is based on the real-time measurements and the probabilistic models of wind speed and solar irradiance. In [35], the authors have been used the genetic algorithm (GA) and the two-point estimation method to solve the OPF problem integrating wind and solar energy. A hybrid method of moth swarm algorithm and gravitational search algorithm (HMSAGSA) is used to solve OPF problems including wind energy [36]. The improved two-point estimation method is proposed to solve the OPF problem combining wind and photovoltaic cells [37]. A self-adaptive evolutionary programming (SAEP) [38] is applied to solve OPF with combined wind energy. A number of other algorithms are proposed to solve OPF problems with integrated renewable energy sources as introduced in [39–46]. In general, these methods are applied in the OPF problem with RES. However, the simulation results only are tested on IEEE 30 bus system with RES while the IEEE 118 bus system is tested without RES.Recently, an algorithm equilibrium optimizer (EO) is developed by Faramarzi et al. [47]. The EO has fast convergence ability, but it is also easy to fall in local minima. In order to overcome the limit, this paper is proposed an improved equilibrium optimizer algorithm (IEO) for solving the OPF problem with integrated renewable energy sources. In the IEO, the exponential term is replaced by a function that does not dependent on the number of iterations. This modification of the IEO algorithm increases exploration ability compared to the EO algorithm. In addition, the exploration of the IEO algorithm will not decrease according to the number of iterations that avoid to get stuck at the local optimal solution. The proposed IEO method is tested on the modified IEEE 30-bus and 118-bus systems with RESs, and simulation results are compared with many other methods. The simulation results obtained show that the proposed IEO algorithm provides quality solution and stability ability with other objective functions. Thus, the proposed method is an alternative approach quite promising for solving optimal power flow problems. The main contributions in this paper are as follows:A new improved IEO algorithm for solving OPF problem with RES, which improves quality solution and stability level, is proposed in this paper.A simple probabilistic formula for calculating RES output power based on the Monte-Carlo simulation model is proposed to reduce the computation time for the OPF with RES.A modified IEEE 118-bus system with wind and solar power plants is introduced in this paper.The remaining organization of the paper is done as follows. Section2 includes the mathematical model related to the OPF problem. In Section 3, uncertainty modelling of wind and solar power outputs is presented. Section 4 describes EO and IEO algorithms and the application of the IEO algorithm for the OPF problem. Section 5 gives simulation results. Finally, the conclusion is given in Section 6.


## 2. Problem Formulation

The two IEEE 30-bus and IEEE 118-bus systems integrating of wind and solar energy are considered in this study. The wind and solar power output are scheduled as like other fossil energy sources. However, the uncertainty of the RES output should require a combination of all generator outputs and reserve outputs for the system to operate in a balanced manner under all circumstances. Therefore, the total cost of generation includes fuel costs for fossil fuel generators, direct costs that ISO buys electricity from RES suppliers, penalty costs, and reserve costs. The generation cost models are described in the following sections.


### 2.1. Cost of Wind and Solar Photovoltaic Power

Wind or solar energy is a random quantity. Therefore, it is difficult to establish an accurate cost model for these energy sources. The authors in [41] have been introduced a model of costs for RES energy sources which is owned by ISO, no fuel costs. However, the excess or lack of this energy at a time is related to the overall operation cost of the power system. Thus, the correctly assessment of the scheduled capacity compared to the actual capacity of the RES sources will reduce the operation costs of the power system. There are two cases when constructing of the cost function for the RES sources. The actual capacity of the RES is less than the scheduled capacity, and the independent system operator (ISO) must dispatch power from other plants to make up for the shortfall of these sources. The cost of committing the reserve generating plants to meet over estimated amount is termed as reserve cost. Opposite, the scheduled capacity is less than the actual power of the RES, and the surplus power of RES will be wasted. This case needs to reduce power output from conventional generators if not possible to utilize RES. The ISO is required to pay a penalty cost corresponding to the surplus amount of RES. The steps to calculate the production cost of RES are described as shown in Figure 1.Figure 1
Cost calculation model for RES source.


#### 2.1.1. Direct Cost of Wind and Solar Photovoltaic Power

Normally, wind or solar farms are owned by private operators. Therefore, the grid operator ISO incurs the cost of scheduled purchases of electricity from these private operators, and this cost is treated as a direct cost.The direct cost of thejth wind power plant in terms of scheduled power is modelled as follows:(1)Cw,jPs.wind,j=gjPs.wind,j,wheregj is coefficient of direct cost and Ps.wind,j is scheduled power of the jth wind power plant.The direct cost of thekth solar power plant in terms of scheduled power is modelled as follows:(2)Cs,kPs.solar,k=hkPs.solar,k,wherehk is the direct cost coefficient and Ps.solar,k is the scheduled power of the kth solar power plant.


#### 2.1.2. Penalty Costs of Wind and Solar Plants

The overestimation of RES occurs when the actual power output is lower than the scheduled power output. Therefore, ISO needs to mobilize from other energy sources and additional costs are incurred. The incurred costs for this mobilization are shown as follows:Reserve cost for thejth wind plant is(3)CwL,j=KwL,j⋅fw,jPwind,j<Ps.wind,j⋅Ps.wind,j−EPwind,j<Ps.wind,jPwind,j,whereCwL,j is the incurred cost for the shortfall of wind power output, Pwind,j is the actual capacity of the wind farm that is less than the scheduled capacityPs.wind,j, fw,jPwind,j<Ps.wind,j is probability of wind power shortage, EPwind,j<Ps.wind,jPwind,j is the scheduled wind power output according to Pwind,j<Ps.wind,j, the scheduled value of the left half-plane in Figure 2, and KwL,j is the coefficient of the incurred cost.Reserve cost for thekth solar plant is(4)CsL,k=KsL,k⋅fs,kPsolar,k<Ps.solar,k⋅Ps.solar,k−EPsolar,k<Ps.solar,kPsolar,k,whereCsL,k is the incurred cost for the shortfall of solar power Ps.solar,k,Psolar,k is the scheduled and actual solar power output, fs,kPsolar,k<Ps.solar,k is the probability of a solar power shortageEPsolar,k<Ps.solar,kPsolar,k is the scheduled power output of solar below Psolar,k<Ps.solar,k, the scheduled value of the left half-plane in Figure 3, and KsL,k is the coefficient of the incurred costs.Figure 2
Wind power distribution for wind farm#1 at bus 5 (C = 9, k = 2).Figure 3
Power output distribution from solar power plant at bus 13 (µ = 6, σ = 0.6).


#### 2.1.3. Opportunity Cost of Wind and Solar Power Surplus

Similarly, the underestimation of RES occurs when the actual energy is higher than the estimated value. As a result, ISO incurs a penalty fee. This cost can be expressed as follows:Penalty cost for the underestimation ofjth wind plant is(5)CwH,j=KwH,j⋅fw,jPwind,j>Ps.wind,j⋅EPwind,j>Ps.wind,jPwind,j−Ps.wind,j,whereCwH,j is the penalty cost of wind power surplus, fw,jPwind,j>Ps.wind,j is the probability of a wind power surplus, EPwind,j>Ps.wind,jPwind,j is the scheduled of power wind output under Pwind,j>Ps.wind,j, the scheduled value of right half-plane in Figure 2 (in kW), and KwH,j is penalty cost of wind power surplus.Penalty cost for the underestimation ofkth solar plant is(6)CsH,k=KsH,k⋅fs,kPsolar,k>Ps.solar,k⋅EPsolar,k>Ps.solar,kPsolar,k−Ps.solar,k,whereCsH,k is the penalty cost of the solar surplus, fs,kPsolar,k>Ps.solar,k is the probability of a solar surplus; EPsolar,k>Ps.solar,kPsolar,k is the scheduled of solar power output under Psolar,k>Ps.solar,k, the scheduled value of right half-plane in Figure 3 (in kW), and KsH,k is penalty cost of solar surplus.The total cost of wind power generated from a wind farm and solar power is described as follows:(7)Cwind=∑j=1NwindCw,jPws,j+CwL,j+CwH,j,Csolar=∑k=1NsolarCs,kPss,k+CsL,k+CsH,k.


### 2.2. Objective Function

In the study, three different objective functions including generation cost, total emission, and power losses are considered in the OPF problem integrating wind and solar power on the traditional grid.


#### 2.2.1. Generation Cost

(8)FC=FC+Cwind+Csolar⋅$h,where FC is cost function of thermal power generators. Cwind,Csolar are the costs of wind and solar energy, respectively.Fuel costs for thermal power plants using fossil fuels can be approximated according to the quadratic function as follows:(9)FC=∑i=1NTGai+biPTGi+ciPTGi2⋅$h.In fact, thermal power generators with multivalve steam turbines exhibit a greater variation in the fuel-cost functions. Therefore, the valve-point effect in terms of absolute value of the sinusoidal function needs to be considered in formula (9):(10)FC=∑i=1NTGai+biPTGi+ciPTGi2+di+sinei⋅PTGimin−PTGi⋅$h,where ai,bi,ci,di, and ei are fuel cost coefficients of the ith generator, while NTG represents the total number of thermal power plants.


#### 2.2.2. Total Emission

The total emissionsFE(ton/h) of thermal power plants are given by(11)FE=∑i=1NTGαi+βiPTGi+γiPTGi2+ξiexpλiPTGi⋅tonh,where αi,βi,γi,ξi, and λi are emission factors of ith thermal power plant.


#### 2.2.3. Total Transmission Loss

The total power lossPloss in the transmission network can be calculated according to the following formula:(12)Ploss=∑i=1NGPGi−∑j=1NLPLjMW.


### 2.3. Constraints

#### 2.3.1. Equality Constraints

Constraints on real and reactive power balance:(13)PGi−PLi=Vi∑j=1NBYijVjcosδi−δj−θij,QGi−QLi=Vi∑j=1NBYijVjsinδi−δj−θij.


#### 2.3.2. Inequality Constraints

The limits of power generation:(14)PGi,min≤PGi≤PGi,max,i=1,2,…,NG,(15)QGi,min≤QGi≤QGi,max,i=1,2,…,NG.The limits of generator voltage bus and load voltage bus:(16)VGi,min≤VGi≤VGi,max,i=1,2,…,NG,(17)VLj,min≤VLj≤VLj,max,j=1,2,…,NL.The limits of transmission line:(18)Sl≤Sl,max,l=1,2,…,NLine.The limits of switchable capacitor capacity:(19)Qci,min≤Qci≤Qci,max,i=1,2,…,NC.The limits of transformer tap:(20)Tk,min≤Tk≤Tk,max,k=1,2,…,NT.


## 2.1. Cost of Wind and Solar Photovoltaic Power

Wind or solar energy is a random quantity. Therefore, it is difficult to establish an accurate cost model for these energy sources. The authors in [41] have been introduced a model of costs for RES energy sources which is owned by ISO, no fuel costs. However, the excess or lack of this energy at a time is related to the overall operation cost of the power system. Thus, the correctly assessment of the scheduled capacity compared to the actual capacity of the RES sources will reduce the operation costs of the power system. There are two cases when constructing of the cost function for the RES sources. The actual capacity of the RES is less than the scheduled capacity, and the independent system operator (ISO) must dispatch power from other plants to make up for the shortfall of these sources. The cost of committing the reserve generating plants to meet over estimated amount is termed as reserve cost. Opposite, the scheduled capacity is less than the actual power of the RES, and the surplus power of RES will be wasted. This case needs to reduce power output from conventional generators if not possible to utilize RES. The ISO is required to pay a penalty cost corresponding to the surplus amount of RES. The steps to calculate the production cost of RES are described as shown in Figure 1.Figure 1
Cost calculation model for RES source.


### 2.1.1. Direct Cost of Wind and Solar Photovoltaic Power

Normally, wind or solar farms are owned by private operators. Therefore, the grid operator ISO incurs the cost of scheduled purchases of electricity from these private operators, and this cost is treated as a direct cost.The direct cost of thejth wind power plant in terms of scheduled power is modelled as follows:(1)Cw,jPs.wind,j=gjPs.wind,j,wheregj is coefficient of direct cost and Ps.wind,j is scheduled power of the jth wind power plant.The direct cost of thekth solar power plant in terms of scheduled power is modelled as follows:(2)Cs,kPs.solar,k=hkPs.solar,k,wherehk is the direct cost coefficient and Ps.solar,k is the scheduled power of the kth solar power plant.


### 2.1.2. Penalty Costs of Wind and Solar Plants

The overestimation of RES occurs when the actual power output is lower than the scheduled power output. Therefore, ISO needs to mobilize from other energy sources and additional costs are incurred. The incurred costs for this mobilization are shown as follows:Reserve cost for thejth wind plant is(3)CwL,j=KwL,j⋅fw,jPwind,j<Ps.wind,j⋅Ps.wind,j−EPwind,j<Ps.wind,jPwind,j,whereCwL,j is the incurred cost for the shortfall of wind power output, Pwind,j is the actual capacity of the wind farm that is less than the scheduled capacityPs.wind,j, fw,jPwind,j<Ps.wind,j is probability of wind power shortage, EPwind,j<Ps.wind,jPwind,j is the scheduled wind power output according to Pwind,j<Ps.wind,j, the scheduled value of the left half-plane in Figure 2, and KwL,j is the coefficient of the incurred cost.Reserve cost for thekth solar plant is(4)CsL,k=KsL,k⋅fs,kPsolar,k<Ps.solar,k⋅Ps.solar,k−EPsolar,k<Ps.solar,kPsolar,k,whereCsL,k is the incurred cost for the shortfall of solar power Ps.solar,k,Psolar,k is the scheduled and actual solar power output, fs,kPsolar,k<Ps.solar,k is the probability of a solar power shortageEPsolar,k<Ps.solar,kPsolar,k is the scheduled power output of solar below Psolar,k<Ps.solar,k, the scheduled value of the left half-plane in Figure 3, and KsL,k is the coefficient of the incurred costs.Figure 2
Wind power distribution for wind farm#1 at bus 5 (C = 9, k = 2).Figure 3
Power output distribution from solar power plant at bus 13 (µ = 6, σ = 0.6).


### 2.1.3. Opportunity Cost of Wind and Solar Power Surplus

Similarly, the underestimation of RES occurs when the actual energy is higher than the estimated value. As a result, ISO incurs a penalty fee. This cost can be expressed as follows:Penalty cost for the underestimation ofjth wind plant is(5)CwH,j=KwH,j⋅fw,jPwind,j>Ps.wind,j⋅EPwind,j>Ps.wind,jPwind,j−Ps.wind,j,whereCwH,j is the penalty cost of wind power surplus, fw,jPwind,j>Ps.wind,j is the probability of a wind power surplus, EPwind,j>Ps.wind,jPwind,j is the scheduled of power wind output under Pwind,j>Ps.wind,j, the scheduled value of right half-plane in Figure 2 (in kW), and KwH,j is penalty cost of wind power surplus.Penalty cost for the underestimation ofkth solar plant is(6)CsH,k=KsH,k⋅fs,kPsolar,k>Ps.solar,k⋅EPsolar,k>Ps.solar,kPsolar,k−Ps.solar,k,whereCsH,k is the penalty cost of the solar surplus, fs,kPsolar,k>Ps.solar,k is the probability of a solar surplus; EPsolar,k>Ps.solar,kPsolar,k is the scheduled of solar power output under Psolar,k>Ps.solar,k, the scheduled value of right half-plane in Figure 3 (in kW), and KsH,k is penalty cost of solar surplus.The total cost of wind power generated from a wind farm and solar power is described as follows:(7)Cwind=∑j=1NwindCw,jPws,j+CwL,j+CwH,j,Csolar=∑k=1NsolarCs,kPss,k+CsL,k+CsH,k.


## 2.1.1. Direct Cost of Wind and Solar Photovoltaic Power

Normally, wind or solar farms are owned by private operators. Therefore, the grid operator ISO incurs the cost of scheduled purchases of electricity from these private operators, and this cost is treated as a direct cost.The direct cost of thejth wind power plant in terms of scheduled power is modelled as follows:(1)Cw,jPs.wind,j=gjPs.wind,j,wheregj is coefficient of direct cost and Ps.wind,j is scheduled power of the jth wind power plant.The direct cost of thekth solar power plant in terms of scheduled power is modelled as follows:(2)Cs,kPs.solar,k=hkPs.solar,k,wherehk is the direct cost coefficient and Ps.solar,k is the scheduled power of the kth solar power plant.


## 2.1.2. Penalty Costs of Wind and Solar Plants

The overestimation of RES occurs when the actual power output is lower than the scheduled power output. Therefore, ISO needs to mobilize from other energy sources and additional costs are incurred. The incurred costs for this mobilization are shown as follows:Reserve cost for thejth wind plant is(3)CwL,j=KwL,j⋅fw,jPwind,j<Ps.wind,j⋅Ps.wind,j−EPwind,j<Ps.wind,jPwind,j,whereCwL,j is the incurred cost for the shortfall of wind power output, Pwind,j is the actual capacity of the wind farm that is less than the scheduled capacityPs.wind,j, fw,jPwind,j<Ps.wind,j is probability of wind power shortage, EPwind,j<Ps.wind,jPwind,j is the scheduled wind power output according to Pwind,j<Ps.wind,j, the scheduled value of the left half-plane in Figure 2, and KwL,j is the coefficient of the incurred cost.Reserve cost for thekth solar plant is(4)CsL,k=KsL,k⋅fs,kPsolar,k<Ps.solar,k⋅Ps.solar,k−EPsolar,k<Ps.solar,kPsolar,k,whereCsL,k is the incurred cost for the shortfall of solar power Ps.solar,k,Psolar,k is the scheduled and actual solar power output, fs,kPsolar,k<Ps.solar,k is the probability of a solar power shortageEPsolar,k<Ps.solar,kPsolar,k is the scheduled power output of solar below Psolar,k<Ps.solar,k, the scheduled value of the left half-plane in Figure 3, and KsL,k is the coefficient of the incurred costs.Figure 2
Wind power distribution for wind farm#1 at bus 5 (C = 9, k = 2).Figure 3
Power output distribution from solar power plant at bus 13 (µ = 6, σ = 0.6).


## 2.1.3. Opportunity Cost of Wind and Solar Power Surplus

Similarly, the underestimation of RES occurs when the actual energy is higher than the estimated value. As a result, ISO incurs a penalty fee. This cost can be expressed as follows:Penalty cost for the underestimation ofjth wind plant is(5)CwH,j=KwH,j⋅fw,jPwind,j>Ps.wind,j⋅EPwind,j>Ps.wind,jPwind,j−Ps.wind,j,whereCwH,j is the penalty cost of wind power surplus, fw,jPwind,j>Ps.wind,j is the probability of a wind power surplus, EPwind,j>Ps.wind,jPwind,j is the scheduled of power wind output under Pwind,j>Ps.wind,j, the scheduled value of right half-plane in Figure 2 (in kW), and KwH,j is penalty cost of wind power surplus.Penalty cost for the underestimation ofkth solar plant is(6)CsH,k=KsH,k⋅fs,kPsolar,k>Ps.solar,k⋅EPsolar,k>Ps.solar,kPsolar,k−Ps.solar,k,whereCsH,k is the penalty cost of the solar surplus, fs,kPsolar,k>Ps.solar,k is the probability of a solar surplus; EPsolar,k>Ps.solar,kPsolar,k is the scheduled of solar power output under Psolar,k>Ps.solar,k, the scheduled value of right half-plane in Figure 3 (in kW), and KsH,k is penalty cost of solar surplus.The total cost of wind power generated from a wind farm and solar power is described as follows:(7)Cwind=∑j=1NwindCw,jPws,j+CwL,j+CwH,j,Csolar=∑k=1NsolarCs,kPss,k+CsL,k+CsH,k.


## 2.2. Objective Function

In the study, three different objective functions including generation cost, total emission, and power losses are considered in the OPF problem integrating wind and solar power on the traditional grid.


### 2.2.1. Generation Cost

(8)FC=FC+Cwind+Csolar⋅$h,where FC is cost function of thermal power generators. Cwind,Csolar are the costs of wind and solar energy, respectively.Fuel costs for thermal power plants using fossil fuels can be approximated according to the quadratic function as follows:(9)FC=∑i=1NTGai+biPTGi+ciPTGi2⋅$h.In fact, thermal power generators with multivalve steam turbines exhibit a greater variation in the fuel-cost functions. Therefore, the valve-point effect in terms of absolute value of the sinusoidal function needs to be considered in formula (9):(10)FC=∑i=1NTGai+biPTGi+ciPTGi2+di+sinei⋅PTGimin−PTGi⋅$h,where ai,bi,ci,di, and ei are fuel cost coefficients of the ith generator, while NTG represents the total number of thermal power plants.


### 2.2.2. Total Emission

The total emissionsFE(ton/h) of thermal power plants are given by(11)FE=∑i=1NTGαi+βiPTGi+γiPTGi2+ξiexpλiPTGi⋅tonh,where αi,βi,γi,ξi, and λi are emission factors of ith thermal power plant.


### 2.2.3. Total Transmission Loss

The total power lossPloss in the transmission network can be calculated according to the following formula:(12)Ploss=∑i=1NGPGi−∑j=1NLPLjMW.


## 2.2.1. Generation Cost

(8)FC=FC+Cwind+Csolar⋅$h,where FC is cost function of thermal power generators. Cwind,Csolar are the costs of wind and solar energy, respectively.Fuel costs for thermal power plants using fossil fuels can be approximated according to the quadratic function as follows:(9)FC=∑i=1NTGai+biPTGi+ciPTGi2⋅$h.In fact, thermal power generators with multivalve steam turbines exhibit a greater variation in the fuel-cost functions. Therefore, the valve-point effect in terms of absolute value of the sinusoidal function needs to be considered in formula (9):(10)FC=∑i=1NTGai+biPTGi+ciPTGi2+di+sinei⋅PTGimin−PTGi⋅$h,where ai,bi,ci,di, and ei are fuel cost coefficients of the ith generator, while NTG represents the total number of thermal power plants.


## 2.2.2. Total Emission

The total emissionsFE(ton/h) of thermal power plants are given by(11)FE=∑i=1NTGαi+βiPTGi+γiPTGi2+ξiexpλiPTGi⋅tonh,where αi,βi,γi,ξi, and λi are emission factors of ith thermal power plant.


## 2.2.3. Total Transmission Loss

The total power lossPloss in the transmission network can be calculated according to the following formula:(12)Ploss=∑i=1NGPGi−∑j=1NLPLjMW.


## 2.3. Constraints

### 2.3.1. Equality Constraints

Constraints on real and reactive power balance:(13)PGi−PLi=Vi∑j=1NBYijVjcosδi−δj−θij,QGi−QLi=Vi∑j=1NBYijVjsinδi−δj−θij.


### 2.3.2. Inequality Constraints

The limits of power generation:(14)PGi,min≤PGi≤PGi,max,i=1,2,…,NG,(15)QGi,min≤QGi≤QGi,max,i=1,2,…,NG.The limits of generator voltage bus and load voltage bus:(16)VGi,min≤VGi≤VGi,max,i=1,2,…,NG,(17)VLj,min≤VLj≤VLj,max,j=1,2,…,NL.The limits of transmission line:(18)Sl≤Sl,max,l=1,2,…,NLine.The limits of switchable capacitor capacity:(19)Qci,min≤Qci≤Qci,max,i=1,2,…,NC.The limits of transformer tap:(20)Tk,min≤Tk≤Tk,max,k=1,2,…,NT.


## 2.3.1. Equality Constraints

Constraints on real and reactive power balance:(13)PGi−PLi=Vi∑j=1NBYijVjcosδi−δj−θij,QGi−QLi=Vi∑j=1NBYijVjsinδi−δj−θij.


## 2.3.2. Inequality Constraints

The limits of power generation:(14)PGi,min≤PGi≤PGi,max,i=1,2,…,NG,(15)QGi,min≤QGi≤QGi,max,i=1,2,…,NG.The limits of generator voltage bus and load voltage bus:(16)VGi,min≤VGi≤VGi,max,i=1,2,…,NG,(17)VLj,min≤VLj≤VLj,max,j=1,2,…,NL.The limits of transmission line:(18)Sl≤Sl,max,l=1,2,…,NLine.The limits of switchable capacitor capacity:(19)Qci,min≤Qci≤Qci,max,i=1,2,…,NC.The limits of transformer tap:(20)Tk,min≤Tk≤Tk,max,k=1,2,…,NT.


## 3. Uncertainty and Power Model of Wind and Solar PV Power Plants

### 3.1. Uncertainty and Power Model of Wind Power Plants

Unlike conventional thermal power generators, the characteristic of the output power generated by the wind turbine is a random variable depending on wind speed. The wind speed distribution follows the Weibull probability density function with a scale factor (C) and a shape factor (k). The probability of wind speed v (m/s) is given by(21)fvv=kCvCk−1e−v/Ckfor0≤v<+∞.The average wind speed is determined by(22)Meanv=C⋅Γ1+k−1,where Γx is gamma function and is described as follows:(23)Γx=∫0∞e−ttx−1dt.Wind speed following the Weibull distribution can be calculated as follows:(24)v=C⋅−lnr1/k,where r is a random number uniformly distributed in the range [0,1].Figures4 and 5 show the best fit of the Weibull distribution and the frequency of wind output at bus 5 and bus 11 using 10000 Monte-Carlo scenarios, respectively. In the study, the value of the shape parameter (k) and scale parameter (C) for the simulation of wind farms is chosen as reference [41]. For the modified IEEE-30 bus system, C = 9 and C = 10 are assumed for two wind farms located at bus 5 and bus 11. The average wind speed of the two farms is 7,976 (m/s) and 8,862 (m/s), respectively, which is completely consistent with the design standards for wind turbine installations specified in [48], when the maximum allowable average wind speed is 10 m/s.Figure 4
Wind speed distribution for wind farm#1 at bus 5 (C = 9, k = 2).Figure 5
Wind speed distribution for wind farm#2 at bus 11 (C = 10, k = 2).Relationship between power output of wind turbine and wind speed (v) is described as(25)Pwindv=0,forv<vinandv>vout,Pwrv−vinvr−vin,forvin<v<vr,Pwr,forvr<v<vout,where vin, vout, and vr are the cut-in wind speed, cut-out wind speed, and rated wind speed of the turbine, respectively. The rated output power of the wind turbine is Pwr.According to formula (25) and the frequency distribution of wind speed, the distribution of power output of the wind farm at bus 5 can be described by the frequency distribution graph as shown in Figure 2. The dotted line shows the scheduled output power needed to supply to the grid.


### 3.2. Uncertainty and Power Model of Solar PV Power Plants

Similarly, the power model of wind plants, the random nature of solar radiation (H) is modelled according to a log-normal distribution function. The solar irradiance probability depends on the standard deviation (µ) and the mean (σ) as follows:(26)fHH=1Hσ2πe−logH−μ2/2σ2.The average solar irradiance is determined by(27)MeanH=eμ+σ2/2.Figure6 shows the fit between the log-normal distribution curve and the solar radiation frequency using 10000 Monte-Carlo scenarios. The solar irradiance (H) to energy conversion for solar PV is given by(28)PsolarH=PsrH2HstdRc,for0<H<Rc,PsrHHstd,forH≥Rc,where Hstd is the solar irradiance in standard environment, Rc is a certain irradiance point, and Psr is the rated output power of the solar PV unit. In this study, Hstd andRc are set to 800 (W/m2) and 120 (W/m2), respectively [41].Figure 6
Solar irradiance distribution from solar power plant at bus 13 (µ = 6, σ = 0.6).Similar to the wind farm, the distribution of power output from the solar power plant at bus 13 can be depicted by a histogram as shown in Figure3.


### 3.3. Conditional Probabilities of Wind and Solar Energy Sources

In this paper, the calculation probability of wind power is proposed as equation (29). Based on the number of RES output power samples obtained from Monte-Carlo simulation, we can calculate the probability of the RES output power when the RES output power is smaller or larger than the scheduled capacity. The probability of RES output power is described as follows:(29)fxA<KorA>K=∑PA<KofA>KNMontecarlo,where X set is the output power of the RES determined through Monte-Carlo simulation, K is the scheduled power obtained from the algorithm’s solutions, fxA<KofA>K is the probability that the RES power is less than or greater than the scheduled power value K, and NMontecarlo is the number of Monte-Carlo samples.


## 3.1. Uncertainty and Power Model of Wind Power Plants

Unlike conventional thermal power generators, the characteristic of the output power generated by the wind turbine is a random variable depending on wind speed. The wind speed distribution follows the Weibull probability density function with a scale factor (C) and a shape factor (k). The probability of wind speed v (m/s) is given by(21)fvv=kCvCk−1e−v/Ckfor0≤v<+∞.The average wind speed is determined by(22)Meanv=C⋅Γ1+k−1,where Γx is gamma function and is described as follows:(23)Γx=∫0∞e−ttx−1dt.Wind speed following the Weibull distribution can be calculated as follows:(24)v=C⋅−lnr1/k,where r is a random number uniformly distributed in the range [0,1].Figures4 and 5 show the best fit of the Weibull distribution and the frequency of wind output at bus 5 and bus 11 using 10000 Monte-Carlo scenarios, respectively. In the study, the value of the shape parameter (k) and scale parameter (C) for the simulation of wind farms is chosen as reference [41]. For the modified IEEE-30 bus system, C = 9 and C = 10 are assumed for two wind farms located at bus 5 and bus 11. The average wind speed of the two farms is 7,976 (m/s) and 8,862 (m/s), respectively, which is completely consistent with the design standards for wind turbine installations specified in [48], when the maximum allowable average wind speed is 10 m/s.Figure 4
Wind speed distribution for wind farm#1 at bus 5 (C = 9, k = 2).Figure 5
Wind speed distribution for wind farm#2 at bus 11 (C = 10, k = 2).Relationship between power output of wind turbine and wind speed (v) is described as(25)Pwindv=0,forv<vinandv>vout,Pwrv−vinvr−vin,forvin<v<vr,Pwr,forvr<v<vout,where vin, vout, and vr are the cut-in wind speed, cut-out wind speed, and rated wind speed of the turbine, respectively. The rated output power of the wind turbine is Pwr.According to formula (25) and the frequency distribution of wind speed, the distribution of power output of the wind farm at bus 5 can be described by the frequency distribution graph as shown in Figure 2. The dotted line shows the scheduled output power needed to supply to the grid.


## 3.2. Uncertainty and Power Model of Solar PV Power Plants

Similarly, the power model of wind plants, the random nature of solar radiation (H) is modelled according to a log-normal distribution function. The solar irradiance probability depends on the standard deviation (µ) and the mean (σ) as follows:(26)fHH=1Hσ2πe−logH−μ2/2σ2.The average solar irradiance is determined by(27)MeanH=eμ+σ2/2.Figure6 shows the fit between the log-normal distribution curve and the solar radiation frequency using 10000 Monte-Carlo scenarios. The solar irradiance (H) to energy conversion for solar PV is given by(28)PsolarH=PsrH2HstdRc,for0<H<Rc,PsrHHstd,forH≥Rc,where Hstd is the solar irradiance in standard environment, Rc is a certain irradiance point, and Psr is the rated output power of the solar PV unit. In this study, Hstd andRc are set to 800 (W/m2) and 120 (W/m2), respectively [41].Figure 6
Solar irradiance distribution from solar power plant at bus 13 (µ = 6, σ = 0.6).Similar to the wind farm, the distribution of power output from the solar power plant at bus 13 can be depicted by a histogram as shown in Figure3.


## 3.3. Conditional Probabilities of Wind and Solar Energy Sources

In this paper, the calculation probability of wind power is proposed as equation (29). Based on the number of RES output power samples obtained from Monte-Carlo simulation, we can calculate the probability of the RES output power when the RES output power is smaller or larger than the scheduled capacity. The probability of RES output power is described as follows:(29)fxA<KorA>K=∑PA<KofA>KNMontecarlo,where X set is the output power of the RES determined through Monte-Carlo simulation, K is the scheduled power obtained from the algorithm’s solutions, fxA<KofA>K is the probability that the RES power is less than or greater than the scheduled power value K, and NMontecarlo is the number of Monte-Carlo samples.


## 4. Proposed Optimization Algorithm

### 4.1. Equilibrium Optimizer Algorithm

The equilibrium optimizer (EO) algorithm [47] is based on determining the dynamic equilibrium for the mass equilibrium model. The mathematical formula is described as follows:(30)GdXdt=QXeq−QX+P,where matter concentration (X) represents the optimized search space. In the EO algorithm, the new particle concentration is updated based on the equilibrium concentration group Xeq, the exponential E, and the rate of generation P. The four main stages of EO can be summarized as follows.


#### 4.1.1. Initialization

In this stage, the concentration of a particle is randomly initialized for the D-dimensional problem and the number ofNpop samples as follows:(31)Xi,j=rand0,1⋅Xmax,j−Xmin,j+Xmin,j,where iε1,2,…,Npop,jε1,2,…,Dim and rand [0, 1] is a random number between 0 and 1.


#### 4.1.2. Determine the Equilibrium Group

The four candidates in theXeq,i equilibrium are considered to contribute to a better EO for exploration, while the average is for mining. The candidates at equilibrium including the four best solutions in the current population are represented in the following equation:(32)Xpool=Xeq1,Xeq2,Xeq3,Xeq4,Xavg,where Xeq1,Xeq2,Xeq3,Xeq4 is the best, second best, third best, and fourth best solution in the population; Xavg is the average candidate (average solution) of the four best individuals and obtained by the following equation:(33)Xavg=Xeq1+Xeq2+Xeq3+Xeq44.


#### 4.1.3. Exponential Term (E)

The exponential term (E) is one of several that changes over time that will help EO strike the right balance between exploration and exploitation:(34)E=α1signr−0.5⋅e−ψt−1.Here, time is defined as a change function of repetition (It) and maximum number of iterations (Max_It), andψ has a value in the range [0, 1]:(35)t=1−ItMax_Itα2It/Max_It.α1,α2 are the two parameters that manage the exploration and exploitation capabilities of the algorithm.


#### 4.1.4. Generation Rate P

As one of the most important terms in the proposed algorithm to provide the correct solution by improving the exploitation phase,(36)P=P0E,where P0 is the initial value defined as follows:(37)P0=PMXeq−ψX,(38)PM=0.5r1,r2≥EP,0,r2≤EP,in which, r1,r2 are random values in the range [0, 1] and EP is a parameter that adjusts the exploration or exploitation ability of the algorithm EP = [0,1]. EP = 1, the algorithm only implements the exploration phase, leading to the difficulty of finding a globally optimal solution. When EP = 0, the algorithm implements the exploitation phase, which easily leads to a local optimal search. Therefore, to balance exploration and exploitation, the EP value is chosen to be 0.5.Finally, updating the concentration in the algorithm (EO) will be as follows:(39)Xk+1=Xeq+Xk−XeqE+PψG1−E.Xeq is the random choice in equilibrium pool and G is considered as unit index, which is normally set to 1.


### 4.2. Improved Equilibrium Optimizer Algorithm (IEO)

The exploration and exploitation in the optimal algorithm greatly impact the solution result. The unbalance between the exploration and exploitation strategy of the algorithm easily leads to fast convergence (nonoptimal solution quality), or slow convergence, and it takes a long time to find a good solution. Therefore, the balance between exploration and exploitation during the optimal search of each algorithm is extremely important. In the EO algorithm, the exponential termE is defined as balance between exploration and exploitation.However, this balance is not uniform throughout the optimization process. The way to generate the exponential term (E) from formulas (34) and (35) shows that the term (E) is dominated by the number of iterations. The relationship between the value of (E) and the number of iterations is shown in Figure 7. The size of the received limit region of E gradually decreases to zero as it approaches Max_It. This means that the algorithm’s exploration will decrease according to the number of iterations. So, the result gets stuck at the local optimal solution. Although the coefficient (ψ) in equation (39) plays a role in helping (EO) to induce the mutation, the probability of mutation is very low since the coefficient ψ must be relatively small. Therefore, to strike a balance between the exploration and exploitation capabilities of the EO, this paper is proposed two adjustments for creating the new solution. First, replace the exponential term (E) in EO with the random value given by equation (40). The periodicity of the sine function allows the creation of a solution around the previous solution. This helps exploit the space between the two solutions. In addition, when the component of 3/2.rand.sign (rand-0.5) has a value >1, it helps the algorithm to explore the space outside two solutions:(40)E=3/2⋅rand⋅signrand−0.5sinrand.Figure 7
Exponential range (E).Furthermore, to increase the exploitation ability around the best solution of IEO, equation (41) is proposed to generate solutions that locate near the best equilibrium solution. The probability of a solution that is generated by this technique or the technique described from equations (34) to (38) is the same:(41)Xk+1=Xeq1k+rand⋅randp1Xpoolk−randp2Xpoolk,where rand is a vector of size [1 dim] containing random values distributed over the interval [0, 1], and randp1Xpoolk and randp2Xpoolk are random individuals in the equilibrium pool population of kth loop.The implementation of the proposed IEO algorithm to solve the OPF problem can be summarized as follows.To apply the algorithm (IEO) to the OPF problem, each candidate is a vector consisting of control variables such as the scheduled power of the generators (excluding the slack generator) and the voltage at generator buses:(42)Xi,j=PG2,…,PGNG,VG1,…,VGNG,i=1,…,Npopandj=1,…,D.Initially, each candidate is randomly generated in the upper and lower values as expression (31). The upper and lower values in expression (31) are taken from expressions (14) and (16).Step 1.
Determine the system data including branch and bus data as well as the limits of buses voltage, branches power, and generators power.Step 2.
Set the parameters for the algorithm IEO including population sizeNpop, maximum number of iterations Max-It, and number of simulation samples Monte-Carlo.Step 3.
CalculateNMontecarlo sample Monte-Carlo, for wind speed and solar radiation according to Weibull and Lognormal distributions respectively. Then, using equations (25) and (28), calculate the power sample set PWind and PSolar, respectively.Step 4.
Randomly initialize the Npop initial candidates using equation (29). With the control variable found in equation (42).Step 5.
Set up four equilibrium candidates and assume each candidate’s target is extremely large.Step 6.
Use the Matpower toolbox to run the power flow for each candidate, and evaluate the corresponding objective function for each candidate as in equation (8), (11), or (12). In this study, the penalty coefficient method is used for inequality constraints to reject unfeasible solutions. The new objective function will be redefined through the objective functions combined with the penalty function given by (43). The penalty function for all state variables in the OPF problem is as follows:(43)Fobj=FFC,FE,Ploss+Penalty,(44)Penalty=λPslPsl−Psllim2+λQ∑i=1NGQGi−QGilim2+λV∑j=1NLVLj−VLjlim2+λL∑l=1NLineSl−Sllim2,where λPsl,λQ,λV,λL are the penalty coefficients for the real power at the slack busPsl, reactive power at all generator buses QGi, voltage magnitude at all load bus VLj, and apparent power transmitted on branchesSl, respectively.Psllim,QGilim,VLjlim,Sllim are the limits of the above variables, respectively.Step 7.
Update candidates in the balance group using equations (32) and (33)Step 8.
Update the respective candidate for the best solutionXBest = Xeq1.Step 9.
Update the new candidates as follows:(45)Xk+1=using to Eq.36from Eq.40,if rand<0.5,using Eq.41,otherwise.Step 10.
Update the candidate limit according to(46)Xik+1Xmin,i,ifXik+1<Xmin,i,Xmax,i,ifXik+1>Xmax,i,Xik+1,otherwise.i=1,…,D.Step 11.
The new candidate will be kept to replace the old candidate, if the objective function value is better, otherwise, it will be discarded.Step 12.
Repeat Steps from 6 to 11 until the Stopping Condition Is SatisfiedStep 13.
The candidate with the smallest objective function value is considered the optimal solution.


## 4.1. Equilibrium Optimizer Algorithm

The equilibrium optimizer (EO) algorithm [47] is based on determining the dynamic equilibrium for the mass equilibrium model. The mathematical formula is described as follows:(30)GdXdt=QXeq−QX+P,where matter concentration (X) represents the optimized search space. In the EO algorithm, the new particle concentration is updated based on the equilibrium concentration group Xeq, the exponential E, and the rate of generation P. The four main stages of EO can be summarized as follows.


### 4.1.1. Initialization

In this stage, the concentration of a particle is randomly initialized for the D-dimensional problem and the number ofNpop samples as follows:(31)Xi,j=rand0,1⋅Xmax,j−Xmin,j+Xmin,j,where iε1,2,…,Npop,jε1,2,…,Dim and rand [0, 1] is a random number between 0 and 1.


### 4.1.2. Determine the Equilibrium Group

The four candidates in theXeq,i equilibrium are considered to contribute to a better EO for exploration, while the average is for mining. The candidates at equilibrium including the four best solutions in the current population are represented in the following equation:(32)Xpool=Xeq1,Xeq2,Xeq3,Xeq4,Xavg,where Xeq1,Xeq2,Xeq3,Xeq4 is the best, second best, third best, and fourth best solution in the population; Xavg is the average candidate (average solution) of the four best individuals and obtained by the following equation:(33)Xavg=Xeq1+Xeq2+Xeq3+Xeq44.


### 4.1.3. Exponential Term (E)

The exponential term (E) is one of several that changes over time that will help EO strike the right balance between exploration and exploitation:(34)E=α1signr−0.5⋅e−ψt−1.Here, time is defined as a change function of repetition (It) and maximum number of iterations (Max_It), andψ has a value in the range [0, 1]:(35)t=1−ItMax_Itα2It/Max_It.α1,α2 are the two parameters that manage the exploration and exploitation capabilities of the algorithm.


### 4.1.4. Generation Rate P

As one of the most important terms in the proposed algorithm to provide the correct solution by improving the exploitation phase,(36)P=P0E,where P0 is the initial value defined as follows:(37)P0=PMXeq−ψX,(38)PM=0.5r1,r2≥EP,0,r2≤EP,in which, r1,r2 are random values in the range [0, 1] and EP is a parameter that adjusts the exploration or exploitation ability of the algorithm EP = [0,1]. EP = 1, the algorithm only implements the exploration phase, leading to the difficulty of finding a globally optimal solution. When EP = 0, the algorithm implements the exploitation phase, which easily leads to a local optimal search. Therefore, to balance exploration and exploitation, the EP value is chosen to be 0.5.Finally, updating the concentration in the algorithm (EO) will be as follows:(39)Xk+1=Xeq+Xk−XeqE+PψG1−E.Xeq is the random choice in equilibrium pool and G is considered as unit index, which is normally set to 1.


## 4.1.1. Initialization

In this stage, the concentration of a particle is randomly initialized for the D-dimensional problem and the number ofNpop samples as follows:(31)Xi,j=rand0,1⋅Xmax,j−Xmin,j+Xmin,j,where iε1,2,…,Npop,jε1,2,…,Dim and rand [0, 1] is a random number between 0 and 1.


## 4.1.2. Determine the Equilibrium Group

The four candidates in theXeq,i equilibrium are considered to contribute to a better EO for exploration, while the average is for mining. The candidates at equilibrium including the four best solutions in the current population are represented in the following equation:(32)Xpool=Xeq1,Xeq2,Xeq3,Xeq4,Xavg,where Xeq1,Xeq2,Xeq3,Xeq4 is the best, second best, third best, and fourth best solution in the population; Xavg is the average candidate (average solution) of the four best individuals and obtained by the following equation:(33)Xavg=Xeq1+Xeq2+Xeq3+Xeq44.


## 4.1.3. Exponential Term (E)

The exponential term (E) is one of several that changes over time that will help EO strike the right balance between exploration and exploitation:(34)E=α1signr−0.5⋅e−ψt−1.Here, time is defined as a change function of repetition (It) and maximum number of iterations (Max_It), andψ has a value in the range [0, 1]:(35)t=1−ItMax_Itα2It/Max_It.α1,α2 are the two parameters that manage the exploration and exploitation capabilities of the algorithm.


## 4.1.4. Generation Rate P

As one of the most important terms in the proposed algorithm to provide the correct solution by improving the exploitation phase,(36)P=P0E,where P0 is the initial value defined as follows:(37)P0=PMXeq−ψX,(38)PM=0.5r1,r2≥EP,0,r2≤EP,in which, r1,r2 are random values in the range [0, 1] and EP is a parameter that adjusts the exploration or exploitation ability of the algorithm EP = [0,1]. EP = 1, the algorithm only implements the exploration phase, leading to the difficulty of finding a globally optimal solution. When EP = 0, the algorithm implements the exploitation phase, which easily leads to a local optimal search. Therefore, to balance exploration and exploitation, the EP value is chosen to be 0.5.Finally, updating the concentration in the algorithm (EO) will be as follows:(39)Xk+1=Xeq+Xk−XeqE+PψG1−E.Xeq is the random choice in equilibrium pool and G is considered as unit index, which is normally set to 1.


## 4.2. Improved Equilibrium Optimizer Algorithm (IEO)

The exploration and exploitation in the optimal algorithm greatly impact the solution result. The unbalance between the exploration and exploitation strategy of the algorithm easily leads to fast convergence (nonoptimal solution quality), or slow convergence, and it takes a long time to find a good solution. Therefore, the balance between exploration and exploitation during the optimal search of each algorithm is extremely important. In the EO algorithm, the exponential termE is defined as balance between exploration and exploitation.However, this balance is not uniform throughout the optimization process. The way to generate the exponential term (E) from formulas (34) and (35) shows that the term (E) is dominated by the number of iterations. The relationship between the value of (E) and the number of iterations is shown in Figure 7. The size of the received limit region of E gradually decreases to zero as it approaches Max_It. This means that the algorithm’s exploration will decrease according to the number of iterations. So, the result gets stuck at the local optimal solution. Although the coefficient (ψ) in equation (39) plays a role in helping (EO) to induce the mutation, the probability of mutation is very low since the coefficient ψ must be relatively small. Therefore, to strike a balance between the exploration and exploitation capabilities of the EO, this paper is proposed two adjustments for creating the new solution. First, replace the exponential term (E) in EO with the random value given by equation (40). The periodicity of the sine function allows the creation of a solution around the previous solution. This helps exploit the space between the two solutions. In addition, when the component of 3/2.rand.sign (rand-0.5) has a value >1, it helps the algorithm to explore the space outside two solutions:(40)E=3/2⋅rand⋅signrand−0.5sinrand.Figure 7
Exponential range (E).Furthermore, to increase the exploitation ability around the best solution of IEO, equation (41) is proposed to generate solutions that locate near the best equilibrium solution. The probability of a solution that is generated by this technique or the technique described from equations (34) to (38) is the same:(41)Xk+1=Xeq1k+rand⋅randp1Xpoolk−randp2Xpoolk,where rand is a vector of size [1 dim] containing random values distributed over the interval [0, 1], and randp1Xpoolk and randp2Xpoolk are random individuals in the equilibrium pool population of kth loop.The implementation of the proposed IEO algorithm to solve the OPF problem can be summarized as follows.To apply the algorithm (IEO) to the OPF problem, each candidate is a vector consisting of control variables such as the scheduled power of the generators (excluding the slack generator) and the voltage at generator buses:(42)Xi,j=PG2,…,PGNG,VG1,…,VGNG,i=1,…,Npopandj=1,…,D.Initially, each candidate is randomly generated in the upper and lower values as expression (31). The upper and lower values in expression (31) are taken from expressions (14) and (16).Step 1.
Determine the system data including branch and bus data as well as the limits of buses voltage, branches power, and generators power.Step 2.
Set the parameters for the algorithm IEO including population sizeNpop, maximum number of iterations Max-It, and number of simulation samples Monte-Carlo.Step 3.
CalculateNMontecarlo sample Monte-Carlo, for wind speed and solar radiation according to Weibull and Lognormal distributions respectively. Then, using equations (25) and (28), calculate the power sample set PWind and PSolar, respectively.Step 4.
Randomly initialize the Npop initial candidates using equation (29). With the control variable found in equation (42).Step 5.
Set up four equilibrium candidates and assume each candidate’s target is extremely large.Step 6.
Use the Matpower toolbox to run the power flow for each candidate, and evaluate the corresponding objective function for each candidate as in equation (8), (11), or (12). In this study, the penalty coefficient method is used for inequality constraints to reject unfeasible solutions. The new objective function will be redefined through the objective functions combined with the penalty function given by (43). The penalty function for all state variables in the OPF problem is as follows:(43)Fobj=FFC,FE,Ploss+Penalty,(44)Penalty=λPslPsl−Psllim2+λQ∑i=1NGQGi−QGilim2+λV∑j=1NLVLj−VLjlim2+λL∑l=1NLineSl−Sllim2,where λPsl,λQ,λV,λL are the penalty coefficients for the real power at the slack busPsl, reactive power at all generator buses QGi, voltage magnitude at all load bus VLj, and apparent power transmitted on branchesSl, respectively.Psllim,QGilim,VLjlim,Sllim are the limits of the above variables, respectively.Step 7.
Update candidates in the balance group using equations (32) and (33)Step 8.
Update the respective candidate for the best solutionXBest = Xeq1.Step 9.
Update the new candidates as follows:(45)Xk+1=using to Eq.36from Eq.40,if rand<0.5,using Eq.41,otherwise.Step 10.
Update the candidate limit according to(46)Xik+1Xmin,i,ifXik+1<Xmin,i,Xmax,i,ifXik+1>Xmax,i,Xik+1,otherwise.i=1,…,D.Step 11.
The new candidate will be kept to replace the old candidate, if the objective function value is better, otherwise, it will be discarded.Step 12.
Repeat Steps from 6 to 11 until the Stopping Condition Is SatisfiedStep 13.
The candidate with the smallest objective function value is considered the optimal solution.


## 5. Simulation Results

The proposed IEO algorithm is tested on the IEEE 30-bus and IEEE 118-bus systems [49, 50]. The simulation results are done using MATLAB software and performed on a computer with Intel Core i5 CPU @ 2.7 GHz and 4 GB RAM. The study cases and control parameters of algorithms are described in Table 1, 2, respectively.Table 1
The description of study cases.
SectionImplemented methodTest systemObjective function5.1EO, AEO, ABC, CSA, TLBO,proposed IEO methodIEEE 30-bus systemMinimize with RES(1) generation cost(2) total emission(3) power loss5.2EO, AEO, ABC, CSA, TLBO,proposed IEO methodIEEE 118-bus systemMinimize generation cost(1) without RES(2) with RESTable 2
The control parameters of the algorithms for the study cases.
MethodsControl parameters valueValueDescriptionAEO——ABCFoodNumber = Np/2Number of food sourcesTLBOTF selected randomly between 1 and 2—CSAPa = 0.25Mutation probabilityEOa1 = 2, a2 = 1 GP = 0.5—IEOGP = 0.5Mutation probabilityCommon parametersPopulation size = 30 for all cases of IEEE 30-bus system.Population size = 50 for all cases of IEEE 118-bus system.Maximum number of iterations = 300 for all cases of IEEE 30-bus system.Maximum number of iterations = 1000 for all cases IEEE 118-bus system.


### 5.1. The Modified IEEE-30 Bus Test System

The data of the modified IEEE 30-bus system are summarized in Table3. The modified IEEE 30-bus system with RES is shown as Figure 8. The thermal power plants at bus 5, bus 11, and bus 13 are replaced by two wind farms and a solar power plant, respectively. The data of wind farms and solar power plant are given in Table 4 [41]. In addition, the cost factors of RESs for direct costs, penalty costs, and reserve costs are shown as Table 5. The cost and emission coefficients of thermal generators for the IEEE 30-bus system are reported in Table 6. The simulation results of the proposed IEO method are compared to AEO, ABC, CSA, TLBO, and EO methods and the exiting other methods with three objective functions including generator cost, total emission, and power loss. To compare and evaluate methods, the study is based on criteria as follows:(1)
Compared the minimum value obtained from the methods with aims of evaluating the quality of the best solutions. To clearly compare, an improvement index (IF) of percentage is used:(47)IF%method=min.of another method−minof IEOmin.of another method.100%,where theIF index has a (+) sign which means that the proposed method is better than other methods. Opposite, the IF index has a (−) sign which means that the proposed method is not better than other methods.(2)
Compared the standard deviation value with aims of evaluating the stability of the methods(3)
Compared the simulation time. The simulation of methods depends on the two factors: the first, the level of complexity in the way of creating new solutions of methods. Each method has its own way to generate a new solution, so the calculation time is also different. With this factor, there has not been a specific study that offers a suitable comparative solution. The second, comparing the total number of solutions generated after the method meets the stopping criterion. The total number of solutions generated depends on the number of populations and the maximum number of iterations. A study in the literature [51] has shown that number of higher population and iterations can significantly improve quality solution, but simulation time will be longer. Especially, for the OPF problem, each new solution is performed by solving the power flow problem through iterative methods. This increases the simulation time when the number of solutions is more. In general, the total number of solutions of the methods is calculated according to one of the following three formulations:(48)TNFs=Npop∗Max−It,(49)TNFs=2∗Npop∗Max−It,(50)TNFs=Npop+Nb∗Max−It,where formulation (48) is applied to calculate TNFs for methods that generate a new solution in each iteration such as EO or IEO. Equation (49) is applied to methods that generate 2 new solutions in each iteration such as CSA, AEO, TLBO, or ABC. Formulation (50) applies to methods with probability of creating 2nd solution in each iteration such as ISSO. Nb is the number of second solutions generated in each iteration, and it is not a fixed value in each iteration. The methods, which calculate TNFs according to formulation (50), will normally calculate the average value of Nb in the total number of runs, from which the TNFs value will be calculated.Table 3
Summary of the modified IEEE 30-bus system.
TermsQuantityDetailsBuses30[3, 49, 50]Branches41[3, 49, 50]Shunts2At bus 5 and bus 24Transformers4At branch (6–9, 6–10, 4–12, 27–28)Thermal generators3Buses: 1 (slack), 2, and 8Wind generators2Buses: 5 and 11Solar PV unit1Bus: 13Control variables11Scheduled real power for 5 generators excluding theone at the slack bus. Bus voltage of all generator buses (6 nos.)Connected load283.4 +j126.2 (MVA)Generator bus voltage limits6[0.95–1.1] p.u.Load bus voltage limits24[0.95–1.05] p.u.Figure 8
Modified IEEE 30-bus system with RES.Table 4
PDF parameters of wind and solar power plants for the IEEE 30-bus system.
Wind power generating plantsSolar PV plantWind farmNo. ofturbinesRated power,Pwr (MW)Weibull PDFparametersWeibull meanRated power,Psr (MW)Lognormal PDFparametersLognormal mean1-(bus 5)2575C = 9, k = 2v¯  = 7.976 m/s50-(bus 13)µ = 6, σ = 0.6H¯=483W/m22-(bus11)2060C = 10, k = 2v¯  = 8.862 m/sTable 5
Energy source cost factor (RES) for the IEEE 30-bus system.
Cost factor for RESWind power generating plantsSolar PV plant(bus 5)(bus 11)(bus 13)Direct cost factor1.61.751.6Penalty cost factor1.51.51.5Reserve cost factor333Table 6
Cost and emission coefficients of thermal generators for the IEEE 30-bus system.
GeneratorbusabcdeαβγξλTG11020.00375180.0374.091−5.5546.490.00026.667TG2201.750.0175160.0382.543−6.0475.6380.00053.333TG3803.250.00834120.0455.326−3.553.380.0022Case1.
Generation cost minimization
In this case, the fuel cost for the thermal power generators is minimized based on formula (10). Table 7 presents the control parameters and optimal value obtained using the AEO, CSA, TLBO, ABC, and EO algorithm and the proposed IEO method with objective of fuel cost. Table 8 provides the best, worst, average, and standard deviation results of the AEO, CSA, TLBO, ABC, and EO and the proposed IEO method after 50 independent runs. It can be noted that from Table 8, the total fuel cost minimization of methods SHADE-SF [41], JS [42], GWO [43], FPA [44], GWO [45], and MFO [46] give better results than the IEO method in terms of the best result. However, the voltage at the load buses is violated voltage limit. The voltage limit at the load buses of these methods is considered in [0.95 p.u.–1.1 p.u.], while the voltage limit at the load buses of the proposed IEO method is [0.95 p.u.–1.05 p.u.]. Furthermore, to obtain the best value, these methods are used the number of solutions lager than the number of solutions of the proposed IEO method from 1,667 to 16,667 times, equivalent to 15000 solutions for GWO [45] and 150000 solutions for SHADE-SF [41] and FPA [44]. The five remaining methods, which are listed in Table 8, have the same system data. It is a necessary condition for the evaluation of solutions to be fair. The result in Table 8 shows that the total fuel cost is achieved 782.0343 ($/h) using the IEO algorithm, which are better than EO (782.0367 $/h), AEO (782,082 $/h), CSA (782.2001 $/h), TLBO (782.0373 $/h), and ABC (782.0560 $/h) methods. IF (%) of all methods are (+) sign. It indicates that the proposed IEO method has better improvement in the quality of the optimal solution with shorter time from 2 to 2.67 times compared to EO, AEO, CSA, TLBO, and ABC methods. In addition, the evaluation of the stability level of the methods is also carried out. The results of the fuel cost target with the minimum value obtained in 50 runs using IEO and other methods are also given in Table 8 and Figure 9. From Figure 9 and the standard deviation values in Table 8, the IEO has higher stability than other methods. The standard deviation value of the IEO (0.0396) is lower than EO (1.4318), AEO (0.1982), CSA (0.1943), TLBO (0.5949), and ABC (0.5875). The convergence speed of the IEO and other methods is shown in Figure 10. As observed, the convergence speed curve of IEO and EO shows that the EO converges earlier than the IEO. However, the EO tends to fall into the local optimum. It can be seen in Figure 10 that from the 250th iteration to the end of the iteration, the EO could not generate a better new solution. Meanwhile, the proposed IEO method still creates many new and better solutions from the 250th iteration to the end of the iteration. In addition, it is noted that using formula (31) in this study to calculate the indirect penalty cost and the reserve cost of the RES greatly reduces the optimal calculation time. As observed from Table 8, the run time of the proposed IEO method is 14.73 s which is smaller than the AEO (53.01 s), CSA (47.18 s), TLBO (52.95 s), ABC (44.48 s), and EO (17.21 s) methods.Table 7
The control parameters and optimal value obtained of the proposed IEO algorithm and other methods for the IEEE 30-bus system with case 1.
Control parametersLimitsCase 1MinMaxAEOCSATLBOABCEOProposed IEOP1 (MW)50200134.9051134.9148134.9086134.9079134.9079134.9079P2 (MW)208028.546427.775227.634027.310827.340427.8087P5 (MW)155043.747044.203644.217644.411844.512344.0873P8 (MW)103510.003410.029410.004110.000010.000010.0000P11 (MW)103036.274136.682836.458336.760736.417836.2702P13 (MW)124035.928535.801236.145035.964236.176136.3030V1 (p.u.)0.951.11.06301.08171.06471.07041.06221.0766V2 (p.u.)0.951.11.04621.06361.04691.05331.04571.0592V5 (p.u.)0.951.11.02061.03951.02201.02871.02231.0338V8 (p.u.)0.951.11.02161.02831.02161.02621.02151.0299V11(p.u.)0.951.11.09981.08501.10001.09431.09991.0887V13(p.u.)0.951.11.05491.04211.05321.04971.05531.0494Fuel cost——440.7292438.3081437.7190436.6351436.7325438.2770Wind cost——243.3587246.3213245.6025247.3081246.484544.5129Solar cost——97.994097.570798.715898.112898.819799.2444Total cost ($/h)——782.082782.2001782.0373782.0560782.0367782.0343Total emissions0.16010.16030.16040.16040.16040.1603PLoss (MW)——6.00446.00705.96775.95545.95455.9771Table 8
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 1.
MethodMinAverageMaxSDNpopMax-itTNFsCPU times (s)IF (%)SHADE-SF [41]782.503—————24000—+0.059897JS [42]781.6387———1001000100000—−0.0506116GWO [43]781.40———30100030000429−0.0811748CSA [43]784.77———30100030000668+0.348599ABC [43]783.81———30100030000704+0.226547FPA [44]777.3298———305000150000—−0.6052129GWO [45]780.995———3050015000—−0.1330738MFO [46]780.485———4050020000—−0.1985048AEO782.082782.3551782.89420.1982303001800053.01+6.099E-03CSA782.2001782.5745783.07200.1943303001800047.18+0.021196TLBO782.0373782.3001785.12280.5949303001800052.95+3.8361E-04ABC782.0560782.3940784.80910.5875303001800044.48+2.7747E-03EO782.0367782.7271789.34571.431830300900017.21+3.0689E-04Proposed IEO782.0343782.0918782.28780.039630300900014.730Figure 9
Convergence characteristic of the proposed IEO and other exiting methods in case 1.Figure 10
The generator cost after 50 runs of IEO and other methods for IEEE 30-bus systems.Case 2.
Generation emission minimization
The best, worst, and average value, and standard deviation after 50 independent runs are shown as Table9. Table 10 presents the control parameters and optimal values of the AEO, CSA, TLBO, ABC, and EO methods and the proposed IEO method. As observed from Table 10, the AEO, TLBO, and EO algorithm, and the IEO method obtained total emission of 0.091059953 (ton/h). Meanwhile, the CSA, TLBO, and ABC algorithms have the best emissions of 0.091060033 (ton/h), 0.091059966 (ton/h), and 0.091063397 (ton/h), respectively. In general, for total emissions target, the IEO method approximates the other methods. However, the IEO algorithm with standard deviation 5.38 × 10−7 that is more stability than other methods is shown in Table 9 and Figure 11.Table 9
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 2.
MethodMinAverageMaxSDTNFsCPU times (s)IF (%)AEO0.0910599530.0910611970.0910843598.42e − 021800054.61+9.598e − 9CSA0.0910600330.0910618690.0910669501.93e − 061800045.58+8.684e − 5TLBO0.0910599660.0910601190.0910610271.79e − 071800054.65+1.396e − 5ABC0.0910633970.0912361060.0918246011.65e − 041800046.61+3.782e − 3EO0.0910599530.0913990580.0995337321.68e − 03900016.64+7.964e − 9Proposed IEO0.0910599530.0910602880.0910621675.38e − 07900013.280Table 10
The control parameters and optimal values obtained from the algorithms for the IEEE 30-bus system with case 2.
Control parametersLimitsCase 2MinMaxAEOCSATLBOABCEOIEOP1 (MW)5020050.000050.000750.000150.027050.000050.0000P2 (MW)208046.629846.621946.622746.438546.761346.6378P5 (MW)155058.699171.295974.959846.992652.105964.9722P8 (MW)103535.000035.000035.000035.000035.000035.0000P11 (MW)103058.703442.008960.000058.252953.706356.0415P13 (MW)124038.026642.452420.300450.000049.979635.5051V1 (p.u.)0.951.11.05541.07451.05481.06151.07541.0989V2 (p.u.)0.951.11.04111.05961.03041.04591.03641.0554V5 (p.u.)0.951.11.00931.03640.98341.01860.99781.0266V8 (p.u.)0.951.10.99001.00390.98251.01521.00340.9982V11 (p.u.)0.951.11.02850.95001.05461.09261.06091.0095V13 (p.u.)0.951.11.07851.07521.04461.05821.00981.0403Fuel cost——377.3874377.3601377.3616376.7612377.8749377.4170Wind cost——388.7138373.3715463.5692341.7264341.0717403.1322Solar cost——105.1171120.992657.6198150.2402150.158396.5911Total cost ($/h)——871.2182871.7242898.5506868.7278869.1050877.1403Total emissions0.091059950.091060030.091059960.091063390.091059950.09105995PLoss (MW)——3.65893.97983.48303.31104.15324.7565Figure 11
The emission after 50 runs of IEO and other methods for the IEEE 30-bus system.Case 3.
Power loss minimization
The control parameters and results of power loss obtained using the proposed IEO algorithm and other methods are shown in Table11. The results of the solutions obtained after 50 independent runs are presented in Table 12. Figure 12 presents stability level of the IEO method compared with other algorithms after 50 independent runs. From Table 11, it can be seen that the power loss of the IEO (2.2211 MW) is lower than AEO (2.2352 MW), CSA(2.2530 MW), TLBO (2.2214 MW), ABC (2.2986 MW), and EO (2.2242 MW). Moreover, the standard deviation of IEO algorithm is better than compared to other methods as shown in Figure 12. Besides, Figure 13 and 14 present the branch power flow and reactive power of generator buses for three cases. From these figures, it can be noted that the branch power flow and reactive power of generator buses are within allowable limit.Table 11
The control parameters and optimal values obtained from the IEO and other methods for the IEEE 30-bus system with case 3.
Control parametersLimitsCase 3MinMaxAEOCSATLBOABCEOIEOP1 (MW)5020050.000051.228650.000050.216250.000050.0010P2 (MW)208029.494023.763426.557426.745225.842126.3951P5 (MW)155074.943074.891074.999875.000074.996975.0000P8 (MW)103534.860935.000034.999535.000034.989435.0000P11 (MW)103057.043158.549257.940160.000058.296257.9242P13 (MW)124039.294142.220841.124638.737341.499541.3009V1 (p.u.)0.951.11.04071.04711.04351.03491.04591.0430V2 (p.u.)0.951.11.03501.03921.03701.02881.03931.0366V5 (p.u.)0.951.11.02361.02931.02571.01581.02691.0256V8 (p.u.)0.951.11.02411.02651.02591.01561.02721.0257V11 (p.u.)0.951.11.10001.09631.09991.07661.09831.1000V13 (p.u.)0.951.11.05921.05611.05741.05381.05611.0579Fuel cost——316.0890301.6589306.9304308.2045304.5591306.4052Wind cost——450.4822456.8491454.6500463.7457456.2010454.5809Solar cost——109.5512120.1371116.1236107.5912117.4901116.7652Total cost ($/h)——876.1223878.6450877.7040879.5414878.2502877.7513Total emissions0.09300.09470.09380.09370.09400.0938PLoss (MW)——2.23522.25302.22142.29862.22422.2211Table 12
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 3.
MethodMinAverageMaxSDTNFsCPU times (s)IF (%)AEO2.23522.30882.49700.06211800050.11+0.630816CSA2.25302.29262.37120.02121800042.38+1.41589TLBO2.22142.25172.30720.02611800053.89+0.013505ABC2.29862.48522.89100.11531800043.17+3.371618EO2.22422.31793.00750.1197900014.84+0.139376Proposed IEO2.22112.23002.26140.0085900012.470Figure 12
Power loss after 50 runs of IEO and other methods for the IEEE 30-bus system.Figure 13
Power flow in branch of three cases for the IEEE 30-bus system.Figure 14
Reactive power of generator buses of the three cases for the IEEE 30-bus system.


### 5.2. The IEEE 118-Bus System

The IEEE 118-bus system is used to test the performance of the IEO in terms of global optimal value and stability ability for a large system. Besides, a modified IEEE 118-bus system with the participation of RES is also considered in this case. The IEEE118-bus system includes 54 thermal units, 99 load buses, 186 branches, 9 transformers, and 12 reactive power compensators [3, 49, 50]. In the modified IEEE 118-bus system, the thermal power plants are replaced by wind farms and solar power plants. The thermal power plants at buses 1, 15, 27, 34, and 36 are replaced by wind farms, while three solar power plants are replaced for thermal power plants at buses 54, 55, and 56. The data of IEEE 118 bus system and the parameters of wind farm and solar plants are given in Table 13 and 14, respectively. The direct cost coefficients, penalty costs, and storage costs of the RES are given in Table 15. In this case, the fuel cost for the thermal power generators is minimized based on formula (9).Table 13
The data of the modified IEEE-118 bus system.
TermsQuantityDetailsBuses118[3, 49, 50]Branches186[3, 49, 50]Shunts14Buses 5, 34, 37, 44, 45, 46, 48, 74, 79, 82, 83, 105 107, and 110Transformers11Branches (8, 32, 36, 51, 93, 95, 102, 107, 127, 134, and 183)Thermal generators46Bus: 69 (slack)Wind generators5Buses: 1, 15, 27, 34, and 36Solar PV unit3Buses: 54, 55, and 56Control variables107Scheduled real power for 54 generators excludingthe one at the slack bus. Bus voltage of all generator buses (54 nos.)Connected load4242.0 +j1438.0 (MVA)Generator bus voltage limits54[0.94–1.06] p.u.Load bus voltage limits64[0.94–1.06] p.u.Table 14
PDF parameters of wind and solar power plants for the IEEE-118 bus system.
Wind power generating plantsSolar PV plantWind farmNo. ofturbinesRated power,Pwr (in MW)Weibull PDFparametersWeibullmeanRated power,Psr (in MW)Lognormal PDFparametersLognormal mean1—(bus 1)2575C = 9, k = 2v¯  = 7.976 m/s50—(bus 54)µ = 6, σ = 0.6H¯=483W/m22—(bus 15)257550—(bus 55)3—(bus 27)257550—(bus 56)4—(bus 34)2060C = 10, k = 2v¯  = 8.862 m/s5—(bus 35)2060Table 15
Energy source cost factor (RES) with the IEEE118-bus system.
Cost factor for energy source RESWind power generating plantsSolar PV plant(bus 1)(bus 15)(bus 27)(bus 34)(bus 35)(bus 54)(bus 55)(bus 56)Direct cost factor1515151515151515Penalty cost factor2020202222182020Reserve cost factor4040404040404040Case 4.
The IEEE 118-bus system without RES
The control parameter and fuel cost obtained from the proposed IEO algorithm for the IEEE 118-bus system without RES are presented in Table16. Table 17 shows the best, average, and worst results after 50 independent runs, as well as the simulation time and number of new solutions to obtain the best result of the proposed IEO algorithm and other methods for case without RES. Figure 15 presents convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system without RES. Besides, the results obtained after 50 independent runs of the proposed IEO method and other methods are shown in Figure 16. From Table 17, it can be seen that the IF index of the methods is signed (+). This shows that the proposed IEO method has the best obtained value of (129820.7252) with only 50000 new solutions which is smaller than AEO, CSA, TLBO, ABC, EO, and BBO [50], ICBO [29], and BSA [50]. The standard deviation achieved from the proposed IEO (245.13772) is smaller than the original EO algorithm (1566.9101). In addition, observed average simulation time of the methods in Table 17 also shows the robustness of the IEO method compared to other methods. The proposed IEO method takes only 132.04(s) time per run while AEO (541.21 s), CSA (512.76 s) TLBO (550.43 s), ABC (478.02 s), and EO (168.82 s).Table 16
The results obtained after 50 independent runs for the 118-bus system without RES of the proposed IEO method and other algorithms.
MethodMinAverageMaxSDNpopMax-ItTNFsCPU times (s)IF (%)AEO129821.5088130675.3582131305.6505442.83023501000100000541.21+0.000604CSA129869.6693131081.6882132295.4056690.62946501000100000512.76+0.037687TLBO129834.0084130263.1053130820.2169300.04403501000100000550.43+0.010231ABC131035.1913132761.3158135201.91931216.3549501000100000478.02+0.926824BBO [50]135263.7289135622.6172136611.2731335.0166902500225000—+4.023994ICBO [29]135121.570135175.672——902500225000—+3.923019BSA [50]135333.4743135502.6493135689.127593.1975902500225000—+4.073456EO129876.0705130728.6536136398.25411566.910150100050000168.82+0.042614Proposed IEO129820.7252130025.2172131153.2047245.1377250100050000132.040Table 17
The results obtained after 50 independent runs for the 118-bus system with RES of the proposed IEO method and other algorithms.
MethodMinAverageMaxSDNpopMax-itTNFCPU times (s)IF (%)AEO129493.8171138800.0481148865.34916091.0705501000100000628.13+0.23759CSA129757.5576140077.2324148943.71475701.9366501000100000579.64+0.44036TLBO129409.2005133942.2353139174.87862796.3609501000100000605.16+0.17235ABC132472.8857149873.3151170494.009710438.8178501000100000562.11+2.48106EO129251.2277138703.0949184497.413714185.905250100050000216.47+0.05034IEO129186.1520130831.7454139865.87951900.886550100050000178.840Figure 15
Convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system without RES.Figure 16
The generator cost after 50 runs of the proposed IEO and other methods for IEEE 118-bus systems without RES.Case 5.
Modified IEEE 118-bus system with RES
A modified IEEE 118 bus system with RES is used to test the proposed IEO algorithm. The control parameters and fuel cost obtained from the IEEE 118-bus system with RES using the IEO method are presented in Table16. Figure 17 shows convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system with RES. Besides, the results obtained after 50 independent runs from the proposed IEO method and other algorithms also are given in Table 18 and Figure 18. From Table 18, it can be seen that the total fuel cost of IEO (129186.1520) is better than EO (129251.2277), AEO (129493.8171), CSA (129757.5576), TLBO (129409.2005), and ABC (132472.8857). As observed from Table 18 and Figure 17 and 18, the IEO algorithm can obtain better values with smaller standard deviations than the EO, AEO, CSA, TLBO, and ABC methods. The simulation results demonstrate the effectiveness and robustness of the IEO method in solving the OPF problem with RES for large systems. In addition, from Table 18, it can be seen that the run time of the proposed IEO method for solving the OPF problem with RES is 178.84(s) which is smaller than EO (216.47 s), AEO (628.13 s), CSA (579.64 s), TLBO (605.16 s), and ABC (562.11 s) methods. It is confirmed that using proposed probability equation (31) for calculating the cost indirectly of RES can reduce the simulation time compared to the other approaches.Figure 17
Convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system with RES.Table 18
The control parameters and fuel cost obtained from the IEEE 118-bus system without RES and with RES using the proposed IEO method.
Bus no.Control parameterWithout RESWith RESControl parameterWithout RESWith RES1PG01 (MW)29.975434.2786V01 (p.u.)0.97111.00074PG04 (MW)0.00020.0000V04 (p.u.)0.99931.02746PG06 (MW)0.00000.4252V06 (p.u.)0.98851.01788PG08 (MW)0.00020.0004V08 (p.u.)1.02071.029510PG010 (MW)400.0368394.7797V10 (p.u.)1.06001.060012PG012 (MW)85.758784.5533V12 (p.u.)0.98481.017215PG015 (MW)20.263235.2944V15 (p.u.)0.98181.030218PG018 (MW)0.00000.0175V18 (p.u.)0.98681.034919PG019 (MW)27.430211.8650V19 (p.u.)0.97911.028924PG24 (MW)0.00000.0000V24 (p.u.)1.02351.052525PG25 (MW)193.5861190.2270V25 (p.u.)1.05691.060026PG26 (MW)280.3155275.1830V26 (p.u.)1.03691.024627PG27 (MW)11.627234.4155V27 (p.u.)1.01441.045131PG31 (MW)7.23037.1524V31 (p.u.)0.99951.032632PG32 (MW)16.41800.0133V32 (p.u.)1.00951.042134PG34 (MW)9.443729.7483V34 (p.u.)1.00551.028736PG36 (MW)10.806029.6131V36 (p.u.)1.00131.023140PG40 (MW)51.275159.9937V40 (p.u.)0.99451.016442PG42 (MW)44.072937.4122V42 (p.u.)1.00141.015246PG46 (MW)19.063318.8211V46 (p.u.)1.01871.011949PG49 (MW)192.7917191.8848V49 (p.u.)1.03631.026554PG54 (MW)49.996933.7510V54 (p.u.)1.01920.995955PG55 (MW)31.883132.3740V55 (p.u.)1.01860.992956PG56 (MW)34.256032.3598V56 (p.u.)1.01850.994259PG59 (MW)149.2530150.1249V59 (p.u.)1.03291.011461PG61 (MW)148.5864147.0525V61 (p.u.)1.03081.018662PG62 (MW)0.00200.0581V62 (p.u.)1.02921.015865PG65 (MW)349.3146348.1722V65 (p.u.)0.99770.994866PG66 (MW)348.6932339.5984V66 (p.u.)1.05061.038169PG69 (MW)453.3021450.1949V69 (p.u.)1.04471.045070PG70 (MW)0.00000.0004V70 (p.u.)1.00691.022372PG72 (MW)0.00000.2237V72 (p.u.)1.01011.060073PG73 (MW)00.0578V73 (p.u.)1.00651.028974PG74 (MW)16.607714.7015V74 (p.u.)0.99230.998176PG76 (MW)19.398717.8025V76 (p.u.)0.98230.980677PG77 (MW)0.00000.2165V77 (p.u.)1.01971.011580PG80 (MW)431.6642426.6092V80 (p.u.)1.03421.019985PG85 (MW)0.00050.0000V85 (p.u.)1.02631.030987PG87 (MW)3.65553.6864V87 (p.u.)1.02731.059889PG89 (MW)500.0391502.9658V89 (p.u.)1.04181.045890PG90 (MW)0.00460.0001V90 (p.u.)1.01831.021191PG91 (MW)0.00000.0022V91 (p.u.)1.01961.022792PG92 (MW)0.00050.0000V92 (p.u.)1.02601.024299PG99 (MW)0.00000.0181V99 (p.u.)1.02321.0044100PG100 (MW)230.7860228.6403V100 (p.u.)1.02581.0109103PG103 (MW)38.225638.4942V103 (p.u.)1.02011.0026104PG104 (MW)0.00000.0000V104 (p.u.)1.01080.9970105PG105 (MW)4.61840V105 (p.u.)1.00920.9968107PG107 (MW)30.303038.6618V107 (p.u.)1.00300.9973110PG110 (MW)8.77538.0400V110 (p.u.)1.01710.9907111PG111 (MW)35.092535.3071V111 (p.u.)1.02950.9939112PG112 (MW)38.049736.1964V112 (p.u.)1.01090.9831113PG113 (MW)0.10550.0326V113 (p.u.)1.00211.0431116PG116 (MW)0.00050.0000V116 (p.u.)0.99450.9942Fuel cost ($/h)—119790.1728Wind cost ($/h)—6163.7140Solar cost ($/h)—3232.2652Total cost ($/h)129820.7252129186.1520Figure 18
The generator cost after 50 runs of the proposed IEO method and other methods for the IEEE 118-bus system with RES.


## 5.1. The Modified IEEE-30 Bus Test System

The data of the modified IEEE 30-bus system are summarized in Table3. The modified IEEE 30-bus system with RES is shown as Figure 8. The thermal power plants at bus 5, bus 11, and bus 13 are replaced by two wind farms and a solar power plant, respectively. The data of wind farms and solar power plant are given in Table 4 [41]. In addition, the cost factors of RESs for direct costs, penalty costs, and reserve costs are shown as Table 5. The cost and emission coefficients of thermal generators for the IEEE 30-bus system are reported in Table 6. The simulation results of the proposed IEO method are compared to AEO, ABC, CSA, TLBO, and EO methods and the exiting other methods with three objective functions including generator cost, total emission, and power loss. To compare and evaluate methods, the study is based on criteria as follows:(1)
Compared the minimum value obtained from the methods with aims of evaluating the quality of the best solutions. To clearly compare, an improvement index (IF) of percentage is used:(47)IF%method=min.of another method−minof IEOmin.of another method.100%,where theIF index has a (+) sign which means that the proposed method is better than other methods. Opposite, the IF index has a (−) sign which means that the proposed method is not better than other methods.(2)
Compared the standard deviation value with aims of evaluating the stability of the methods(3)
Compared the simulation time. The simulation of methods depends on the two factors: the first, the level of complexity in the way of creating new solutions of methods. Each method has its own way to generate a new solution, so the calculation time is also different. With this factor, there has not been a specific study that offers a suitable comparative solution. The second, comparing the total number of solutions generated after the method meets the stopping criterion. The total number of solutions generated depends on the number of populations and the maximum number of iterations. A study in the literature [51] has shown that number of higher population and iterations can significantly improve quality solution, but simulation time will be longer. Especially, for the OPF problem, each new solution is performed by solving the power flow problem through iterative methods. This increases the simulation time when the number of solutions is more. In general, the total number of solutions of the methods is calculated according to one of the following three formulations:(48)TNFs=Npop∗Max−It,(49)TNFs=2∗Npop∗Max−It,(50)TNFs=Npop+Nb∗Max−It,where formulation (48) is applied to calculate TNFs for methods that generate a new solution in each iteration such as EO or IEO. Equation (49) is applied to methods that generate 2 new solutions in each iteration such as CSA, AEO, TLBO, or ABC. Formulation (50) applies to methods with probability of creating 2nd solution in each iteration such as ISSO. Nb is the number of second solutions generated in each iteration, and it is not a fixed value in each iteration. The methods, which calculate TNFs according to formulation (50), will normally calculate the average value of Nb in the total number of runs, from which the TNFs value will be calculated.Table 3
Summary of the modified IEEE 30-bus system.
TermsQuantityDetailsBuses30[3, 49, 50]Branches41[3, 49, 50]Shunts2At bus 5 and bus 24Transformers4At branch (6–9, 6–10, 4–12, 27–28)Thermal generators3Buses: 1 (slack), 2, and 8Wind generators2Buses: 5 and 11Solar PV unit1Bus: 13Control variables11Scheduled real power for 5 generators excluding theone at the slack bus. Bus voltage of all generator buses (6 nos.)Connected load283.4 +j126.2 (MVA)Generator bus voltage limits6[0.95–1.1] p.u.Load bus voltage limits24[0.95–1.05] p.u.Figure 8
Modified IEEE 30-bus system with RES.Table 4
PDF parameters of wind and solar power plants for the IEEE 30-bus system.
Wind power generating plantsSolar PV plantWind farmNo. ofturbinesRated power,Pwr (MW)Weibull PDFparametersWeibull meanRated power,Psr (MW)Lognormal PDFparametersLognormal mean1-(bus 5)2575C = 9, k = 2v¯  = 7.976 m/s50-(bus 13)µ = 6, σ = 0.6H¯=483W/m22-(bus11)2060C = 10, k = 2v¯  = 8.862 m/sTable 5
Energy source cost factor (RES) for the IEEE 30-bus system.
Cost factor for RESWind power generating plantsSolar PV plant(bus 5)(bus 11)(bus 13)Direct cost factor1.61.751.6Penalty cost factor1.51.51.5Reserve cost factor333Table 6
Cost and emission coefficients of thermal generators for the IEEE 30-bus system.
GeneratorbusabcdeαβγξλTG11020.00375180.0374.091−5.5546.490.00026.667TG2201.750.0175160.0382.543−6.0475.6380.00053.333TG3803.250.00834120.0455.326−3.553.380.0022Case1.
Generation cost minimization
In this case, the fuel cost for the thermal power generators is minimized based on formula (10). Table 7 presents the control parameters and optimal value obtained using the AEO, CSA, TLBO, ABC, and EO algorithm and the proposed IEO method with objective of fuel cost. Table 8 provides the best, worst, average, and standard deviation results of the AEO, CSA, TLBO, ABC, and EO and the proposed IEO method after 50 independent runs. It can be noted that from Table 8, the total fuel cost minimization of methods SHADE-SF [41], JS [42], GWO [43], FPA [44], GWO [45], and MFO [46] give better results than the IEO method in terms of the best result. However, the voltage at the load buses is violated voltage limit. The voltage limit at the load buses of these methods is considered in [0.95 p.u.–1.1 p.u.], while the voltage limit at the load buses of the proposed IEO method is [0.95 p.u.–1.05 p.u.]. Furthermore, to obtain the best value, these methods are used the number of solutions lager than the number of solutions of the proposed IEO method from 1,667 to 16,667 times, equivalent to 15000 solutions for GWO [45] and 150000 solutions for SHADE-SF [41] and FPA [44]. The five remaining methods, which are listed in Table 8, have the same system data. It is a necessary condition for the evaluation of solutions to be fair. The result in Table 8 shows that the total fuel cost is achieved 782.0343 ($/h) using the IEO algorithm, which are better than EO (782.0367 $/h), AEO (782,082 $/h), CSA (782.2001 $/h), TLBO (782.0373 $/h), and ABC (782.0560 $/h) methods. IF (%) of all methods are (+) sign. It indicates that the proposed IEO method has better improvement in the quality of the optimal solution with shorter time from 2 to 2.67 times compared to EO, AEO, CSA, TLBO, and ABC methods. In addition, the evaluation of the stability level of the methods is also carried out. The results of the fuel cost target with the minimum value obtained in 50 runs using IEO and other methods are also given in Table 8 and Figure 9. From Figure 9 and the standard deviation values in Table 8, the IEO has higher stability than other methods. The standard deviation value of the IEO (0.0396) is lower than EO (1.4318), AEO (0.1982), CSA (0.1943), TLBO (0.5949), and ABC (0.5875). The convergence speed of the IEO and other methods is shown in Figure 10. As observed, the convergence speed curve of IEO and EO shows that the EO converges earlier than the IEO. However, the EO tends to fall into the local optimum. It can be seen in Figure 10 that from the 250th iteration to the end of the iteration, the EO could not generate a better new solution. Meanwhile, the proposed IEO method still creates many new and better solutions from the 250th iteration to the end of the iteration. In addition, it is noted that using formula (31) in this study to calculate the indirect penalty cost and the reserve cost of the RES greatly reduces the optimal calculation time. As observed from Table 8, the run time of the proposed IEO method is 14.73 s which is smaller than the AEO (53.01 s), CSA (47.18 s), TLBO (52.95 s), ABC (44.48 s), and EO (17.21 s) methods.Table 7
The control parameters and optimal value obtained of the proposed IEO algorithm and other methods for the IEEE 30-bus system with case 1.
Control parametersLimitsCase 1MinMaxAEOCSATLBOABCEOProposed IEOP1 (MW)50200134.9051134.9148134.9086134.9079134.9079134.9079P2 (MW)208028.546427.775227.634027.310827.340427.8087P5 (MW)155043.747044.203644.217644.411844.512344.0873P8 (MW)103510.003410.029410.004110.000010.000010.0000P11 (MW)103036.274136.682836.458336.760736.417836.2702P13 (MW)124035.928535.801236.145035.964236.176136.3030V1 (p.u.)0.951.11.06301.08171.06471.07041.06221.0766V2 (p.u.)0.951.11.04621.06361.04691.05331.04571.0592V5 (p.u.)0.951.11.02061.03951.02201.02871.02231.0338V8 (p.u.)0.951.11.02161.02831.02161.02621.02151.0299V11(p.u.)0.951.11.09981.08501.10001.09431.09991.0887V13(p.u.)0.951.11.05491.04211.05321.04971.05531.0494Fuel cost——440.7292438.3081437.7190436.6351436.7325438.2770Wind cost——243.3587246.3213245.6025247.3081246.484544.5129Solar cost——97.994097.570798.715898.112898.819799.2444Total cost ($/h)——782.082782.2001782.0373782.0560782.0367782.0343Total emissions0.16010.16030.16040.16040.16040.1603PLoss (MW)——6.00446.00705.96775.95545.95455.9771Table 8
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 1.
MethodMinAverageMaxSDNpopMax-itTNFsCPU times (s)IF (%)SHADE-SF [41]782.503—————24000—+0.059897JS [42]781.6387———1001000100000—−0.0506116GWO [43]781.40———30100030000429−0.0811748CSA [43]784.77———30100030000668+0.348599ABC [43]783.81———30100030000704+0.226547FPA [44]777.3298———305000150000—−0.6052129GWO [45]780.995———3050015000—−0.1330738MFO [46]780.485———4050020000—−0.1985048AEO782.082782.3551782.89420.1982303001800053.01+6.099E-03CSA782.2001782.5745783.07200.1943303001800047.18+0.021196TLBO782.0373782.3001785.12280.5949303001800052.95+3.8361E-04ABC782.0560782.3940784.80910.5875303001800044.48+2.7747E-03EO782.0367782.7271789.34571.431830300900017.21+3.0689E-04Proposed IEO782.0343782.0918782.28780.039630300900014.730Figure 9
Convergence characteristic of the proposed IEO and other exiting methods in case 1.Figure 10
The generator cost after 50 runs of IEO and other methods for IEEE 30-bus systems.Case 2.
Generation emission minimization
The best, worst, and average value, and standard deviation after 50 independent runs are shown as Table9. Table 10 presents the control parameters and optimal values of the AEO, CSA, TLBO, ABC, and EO methods and the proposed IEO method. As observed from Table 10, the AEO, TLBO, and EO algorithm, and the IEO method obtained total emission of 0.091059953 (ton/h). Meanwhile, the CSA, TLBO, and ABC algorithms have the best emissions of 0.091060033 (ton/h), 0.091059966 (ton/h), and 0.091063397 (ton/h), respectively. In general, for total emissions target, the IEO method approximates the other methods. However, the IEO algorithm with standard deviation 5.38 × 10−7 that is more stability than other methods is shown in Table 9 and Figure 11.Table 9
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 2.
MethodMinAverageMaxSDTNFsCPU times (s)IF (%)AEO0.0910599530.0910611970.0910843598.42e − 021800054.61+9.598e − 9CSA0.0910600330.0910618690.0910669501.93e − 061800045.58+8.684e − 5TLBO0.0910599660.0910601190.0910610271.79e − 071800054.65+1.396e − 5ABC0.0910633970.0912361060.0918246011.65e − 041800046.61+3.782e − 3EO0.0910599530.0913990580.0995337321.68e − 03900016.64+7.964e − 9Proposed IEO0.0910599530.0910602880.0910621675.38e − 07900013.280Table 10
The control parameters and optimal values obtained from the algorithms for the IEEE 30-bus system with case 2.
Control parametersLimitsCase 2MinMaxAEOCSATLBOABCEOIEOP1 (MW)5020050.000050.000750.000150.027050.000050.0000P2 (MW)208046.629846.621946.622746.438546.761346.6378P5 (MW)155058.699171.295974.959846.992652.105964.9722P8 (MW)103535.000035.000035.000035.000035.000035.0000P11 (MW)103058.703442.008960.000058.252953.706356.0415P13 (MW)124038.026642.452420.300450.000049.979635.5051V1 (p.u.)0.951.11.05541.07451.05481.06151.07541.0989V2 (p.u.)0.951.11.04111.05961.03041.04591.03641.0554V5 (p.u.)0.951.11.00931.03640.98341.01860.99781.0266V8 (p.u.)0.951.10.99001.00390.98251.01521.00340.9982V11 (p.u.)0.951.11.02850.95001.05461.09261.06091.0095V13 (p.u.)0.951.11.07851.07521.04461.05821.00981.0403Fuel cost——377.3874377.3601377.3616376.7612377.8749377.4170Wind cost——388.7138373.3715463.5692341.7264341.0717403.1322Solar cost——105.1171120.992657.6198150.2402150.158396.5911Total cost ($/h)——871.2182871.7242898.5506868.7278869.1050877.1403Total emissions0.091059950.091060030.091059960.091063390.091059950.09105995PLoss (MW)——3.65893.97983.48303.31104.15324.7565Figure 11
The emission after 50 runs of IEO and other methods for the IEEE 30-bus system.Case 3.
Power loss minimization
The control parameters and results of power loss obtained using the proposed IEO algorithm and other methods are shown in Table11. The results of the solutions obtained after 50 independent runs are presented in Table 12. Figure 12 presents stability level of the IEO method compared with other algorithms after 50 independent runs. From Table 11, it can be seen that the power loss of the IEO (2.2211 MW) is lower than AEO (2.2352 MW), CSA(2.2530 MW), TLBO (2.2214 MW), ABC (2.2986 MW), and EO (2.2242 MW). Moreover, the standard deviation of IEO algorithm is better than compared to other methods as shown in Figure 12. Besides, Figure 13 and 14 present the branch power flow and reactive power of generator buses for three cases. From these figures, it can be noted that the branch power flow and reactive power of generator buses are within allowable limit.Table 11
The control parameters and optimal values obtained from the IEO and other methods for the IEEE 30-bus system with case 3.
Control parametersLimitsCase 3MinMaxAEOCSATLBOABCEOIEOP1 (MW)5020050.000051.228650.000050.216250.000050.0010P2 (MW)208029.494023.763426.557426.745225.842126.3951P5 (MW)155074.943074.891074.999875.000074.996975.0000P8 (MW)103534.860935.000034.999535.000034.989435.0000P11 (MW)103057.043158.549257.940160.000058.296257.9242P13 (MW)124039.294142.220841.124638.737341.499541.3009V1 (p.u.)0.951.11.04071.04711.04351.03491.04591.0430V2 (p.u.)0.951.11.03501.03921.03701.02881.03931.0366V5 (p.u.)0.951.11.02361.02931.02571.01581.02691.0256V8 (p.u.)0.951.11.02411.02651.02591.01561.02721.0257V11 (p.u.)0.951.11.10001.09631.09991.07661.09831.1000V13 (p.u.)0.951.11.05921.05611.05741.05381.05611.0579Fuel cost——316.0890301.6589306.9304308.2045304.5591306.4052Wind cost——450.4822456.8491454.6500463.7457456.2010454.5809Solar cost——109.5512120.1371116.1236107.5912117.4901116.7652Total cost ($/h)——876.1223878.6450877.7040879.5414878.2502877.7513Total emissions0.09300.09470.09380.09370.09400.0938PLoss (MW)——2.23522.25302.22142.29862.22422.2211Table 12
The obtained results of the proposed IEO method and other exiting methods for the IEEE 30-bus system after 50 runs with case 3.
MethodMinAverageMaxSDTNFsCPU times (s)IF (%)AEO2.23522.30882.49700.06211800050.11+0.630816CSA2.25302.29262.37120.02121800042.38+1.41589TLBO2.22142.25172.30720.02611800053.89+0.013505ABC2.29862.48522.89100.11531800043.17+3.371618EO2.22422.31793.00750.1197900014.84+0.139376Proposed IEO2.22112.23002.26140.0085900012.470Figure 12
Power loss after 50 runs of IEO and other methods for the IEEE 30-bus system.Figure 13
Power flow in branch of three cases for the IEEE 30-bus system.Figure 14
Reactive power of generator buses of the three cases for the IEEE 30-bus system.


## 5.2. The IEEE 118-Bus System

The IEEE 118-bus system is used to test the performance of the IEO in terms of global optimal value and stability ability for a large system. Besides, a modified IEEE 118-bus system with the participation of RES is also considered in this case. The IEEE118-bus system includes 54 thermal units, 99 load buses, 186 branches, 9 transformers, and 12 reactive power compensators [3, 49, 50]. In the modified IEEE 118-bus system, the thermal power plants are replaced by wind farms and solar power plants. The thermal power plants at buses 1, 15, 27, 34, and 36 are replaced by wind farms, while three solar power plants are replaced for thermal power plants at buses 54, 55, and 56. The data of IEEE 118 bus system and the parameters of wind farm and solar plants are given in Table 13 and 14, respectively. The direct cost coefficients, penalty costs, and storage costs of the RES are given in Table 15. In this case, the fuel cost for the thermal power generators is minimized based on formula (9).Table 13
The data of the modified IEEE-118 bus system.
TermsQuantityDetailsBuses118[3, 49, 50]Branches186[3, 49, 50]Shunts14Buses 5, 34, 37, 44, 45, 46, 48, 74, 79, 82, 83, 105 107, and 110Transformers11Branches (8, 32, 36, 51, 93, 95, 102, 107, 127, 134, and 183)Thermal generators46Bus: 69 (slack)Wind generators5Buses: 1, 15, 27, 34, and 36Solar PV unit3Buses: 54, 55, and 56Control variables107Scheduled real power for 54 generators excludingthe one at the slack bus. Bus voltage of all generator buses (54 nos.)Connected load4242.0 +j1438.0 (MVA)Generator bus voltage limits54[0.94–1.06] p.u.Load bus voltage limits64[0.94–1.06] p.u.Table 14
PDF parameters of wind and solar power plants for the IEEE-118 bus system.
Wind power generating plantsSolar PV plantWind farmNo. ofturbinesRated power,Pwr (in MW)Weibull PDFparametersWeibullmeanRated power,Psr (in MW)Lognormal PDFparametersLognormal mean1—(bus 1)2575C = 9, k = 2v¯  = 7.976 m/s50—(bus 54)µ = 6, σ = 0.6H¯=483W/m22—(bus 15)257550—(bus 55)3—(bus 27)257550—(bus 56)4—(bus 34)2060C = 10, k = 2v¯  = 8.862 m/s5—(bus 35)2060Table 15
Energy source cost factor (RES) with the IEEE118-bus system.
Cost factor for energy source RESWind power generating plantsSolar PV plant(bus 1)(bus 15)(bus 27)(bus 34)(bus 35)(bus 54)(bus 55)(bus 56)Direct cost factor1515151515151515Penalty cost factor2020202222182020Reserve cost factor4040404040404040Case 4.
The IEEE 118-bus system without RES
The control parameter and fuel cost obtained from the proposed IEO algorithm for the IEEE 118-bus system without RES are presented in Table16. Table 17 shows the best, average, and worst results after 50 independent runs, as well as the simulation time and number of new solutions to obtain the best result of the proposed IEO algorithm and other methods for case without RES. Figure 15 presents convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system without RES. Besides, the results obtained after 50 independent runs of the proposed IEO method and other methods are shown in Figure 16. From Table 17, it can be seen that the IF index of the methods is signed (+). This shows that the proposed IEO method has the best obtained value of (129820.7252) with only 50000 new solutions which is smaller than AEO, CSA, TLBO, ABC, EO, and BBO [50], ICBO [29], and BSA [50]. The standard deviation achieved from the proposed IEO (245.13772) is smaller than the original EO algorithm (1566.9101). In addition, observed average simulation time of the methods in Table 17 also shows the robustness of the IEO method compared to other methods. The proposed IEO method takes only 132.04(s) time per run while AEO (541.21 s), CSA (512.76 s) TLBO (550.43 s), ABC (478.02 s), and EO (168.82 s).Table 16
The results obtained after 50 independent runs for the 118-bus system without RES of the proposed IEO method and other algorithms.
MethodMinAverageMaxSDNpopMax-ItTNFsCPU times (s)IF (%)AEO129821.5088130675.3582131305.6505442.83023501000100000541.21+0.000604CSA129869.6693131081.6882132295.4056690.62946501000100000512.76+0.037687TLBO129834.0084130263.1053130820.2169300.04403501000100000550.43+0.010231ABC131035.1913132761.3158135201.91931216.3549501000100000478.02+0.926824BBO [50]135263.7289135622.6172136611.2731335.0166902500225000—+4.023994ICBO [29]135121.570135175.672——902500225000—+3.923019BSA [50]135333.4743135502.6493135689.127593.1975902500225000—+4.073456EO129876.0705130728.6536136398.25411566.910150100050000168.82+0.042614Proposed IEO129820.7252130025.2172131153.2047245.1377250100050000132.040Table 17
The results obtained after 50 independent runs for the 118-bus system with RES of the proposed IEO method and other algorithms.
MethodMinAverageMaxSDNpopMax-itTNFCPU times (s)IF (%)AEO129493.8171138800.0481148865.34916091.0705501000100000628.13+0.23759CSA129757.5576140077.2324148943.71475701.9366501000100000579.64+0.44036TLBO129409.2005133942.2353139174.87862796.3609501000100000605.16+0.17235ABC132472.8857149873.3151170494.009710438.8178501000100000562.11+2.48106EO129251.2277138703.0949184497.413714185.905250100050000216.47+0.05034IEO129186.1520130831.7454139865.87951900.886550100050000178.840Figure 15
Convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system without RES.Figure 16
The generator cost after 50 runs of the proposed IEO and other methods for IEEE 118-bus systems without RES.Case 5.
Modified IEEE 118-bus system with RES
A modified IEEE 118 bus system with RES is used to test the proposed IEO algorithm. The control parameters and fuel cost obtained from the IEEE 118-bus system with RES using the IEO method are presented in Table16. Figure 17 shows convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system with RES. Besides, the results obtained after 50 independent runs from the proposed IEO method and other algorithms also are given in Table 18 and Figure 18. From Table 18, it can be seen that the total fuel cost of IEO (129186.1520) is better than EO (129251.2277), AEO (129493.8171), CSA (129757.5576), TLBO (129409.2005), and ABC (132472.8857). As observed from Table 18 and Figure 17 and 18, the IEO algorithm can obtain better values with smaller standard deviations than the EO, AEO, CSA, TLBO, and ABC methods. The simulation results demonstrate the effectiveness and robustness of the IEO method in solving the OPF problem with RES for large systems. In addition, from Table 18, it can be seen that the run time of the proposed IEO method for solving the OPF problem with RES is 178.84(s) which is smaller than EO (216.47 s), AEO (628.13 s), CSA (579.64 s), TLBO (605.16 s), and ABC (562.11 s) methods. It is confirmed that using proposed probability equation (31) for calculating the cost indirectly of RES can reduce the simulation time compared to the other approaches.Figure 17
Convergence characteristic of the proposed IEO method and other methods for the IEEE118-bus system with RES.Table 18
The control parameters and fuel cost obtained from the IEEE 118-bus system without RES and with RES using the proposed IEO method.
Bus no.Control parameterWithout RESWith RESControl parameterWithout RESWith RES1PG01 (MW)29.975434.2786V01 (p.u.)0.97111.00074PG04 (MW)0.00020.0000V04 (p.u.)0.99931.02746PG06 (MW)0.00000.4252V06 (p.u.)0.98851.01788PG08 (MW)0.00020.0004V08 (p.u.)1.02071.029510PG010 (MW)400.0368394.7797V10 (p.u.)1.06001.060012PG012 (MW)85.758784.5533V12 (p.u.)0.98481.017215PG015 (MW)20.263235.2944V15 (p.u.)0.98181.030218PG018 (MW)0.00000.0175V18 (p.u.)0.98681.034919PG019 (MW)27.430211.8650V19 (p.u.)0.97911.028924PG24 (MW)0.00000.0000V24 (p.u.)1.02351.052525PG25 (MW)193.5861190.2270V25 (p.u.)1.05691.060026PG26 (MW)280.3155275.1830V26 (p.u.)1.03691.024627PG27 (MW)11.627234.4155V27 (p.u.)1.01441.045131PG31 (MW)7.23037.1524V31 (p.u.)0.99951.032632PG32 (MW)16.41800.0133V32 (p.u.)1.00951.042134PG34 (MW)9.443729.7483V34 (p.u.)1.00551.028736PG36 (MW)10.806029.6131V36 (p.u.)1.00131.023140PG40 (MW)51.275159.9937V40 (p.u.)0.99451.016442PG42 (MW)44.072937.4122V42 (p.u.)1.00141.015246PG46 (MW)19.063318.8211V46 (p.u.)1.01871.011949PG49 (MW)192.7917191.8848V49 (p.u.)1.03631.026554PG54 (MW)49.996933.7510V54 (p.u.)1.01920.995955PG55 (MW)31.883132.3740V55 (p.u.)1.01860.992956PG56 (MW)34.256032.3598V56 (p.u.)1.01850.994259PG59 (MW)149.2530150.1249V59 (p.u.)1.03291.011461PG61 (MW)148.5864147.0525V61 (p.u.)1.03081.018662PG62 (MW)0.00200.0581V62 (p.u.)1.02921.015865PG65 (MW)349.3146348.1722V65 (p.u.)0.99770.994866PG66 (MW)348.6932339.5984V66 (p.u.)1.05061.038169PG69 (MW)453.3021450.1949V69 (p.u.)1.04471.045070PG70 (MW)0.00000.0004V70 (p.u.)1.00691.022372PG72 (MW)0.00000.2237V72 (p.u.)1.01011.060073PG73 (MW)00.0578V73 (p.u.)1.00651.028974PG74 (MW)16.607714.7015V74 (p.u.)0.99230.998176PG76 (MW)19.398717.8025V76 (p.u.)0.98230.980677PG77 (MW)0.00000.2165V77 (p.u.)1.01971.011580PG80 (MW)431.6642426.6092V80 (p.u.)1.03421.019985PG85 (MW)0.00050.0000V85 (p.u.)1.02631.030987PG87 (MW)3.65553.6864V87 (p.u.)1.02731.059889PG89 (MW)500.0391502.9658V89 (p.u.)1.04181.045890PG90 (MW)0.00460.0001V90 (p.u.)1.01831.021191PG91 (MW)0.00000.0022V91 (p.u.)1.01961.022792PG92 (MW)0.00050.0000V92 (p.u.)1.02601.024299PG99 (MW)0.00000.0181V99 (p.u.)1.02321.0044100PG100 (MW)230.7860228.6403V100 (p.u.)1.02581.0109103PG103 (MW)38.225638.4942V103 (p.u.)1.02011.0026104PG104 (MW)0.00000.0000V104 (p.u.)1.01080.9970105PG105 (MW)4.61840V105 (p.u.)1.00920.9968107PG107 (MW)30.303038.6618V107 (p.u.)1.00300.9973110PG110 (MW)8.77538.0400V110 (p.u.)1.01710.9907111PG111 (MW)35.092535.3071V111 (p.u.)1.02950.9939112PG112 (MW)38.049736.1964V112 (p.u.)1.01090.9831113PG113 (MW)0.10550.0326V113 (p.u.)1.00211.0431116PG116 (MW)0.00050.0000V116 (p.u.)0.99450.9942Fuel cost ($/h)—119790.1728Wind cost ($/h)—6163.7140Solar cost ($/h)—3232.2652Total cost ($/h)129820.7252129186.1520Figure 18
The generator cost after 50 runs of the proposed IEO method and other methods for the IEEE 118-bus system with RES.


## 6. Conclusion

Finding the best solution for OPF problem with RES is one of challenges for many algorithms, especially in complex systems. In this paper, the IEO algorithm is proposed to deal with OPF problem with RES for three other objective functions. The suggested IEO technique is tested on IEEE 30 and IEEE 118 bus system. The obtained results of suggested approach are compared with EO, AEO, CSA, TLBO, and ABC algorithms and other exiting methods. For the IEEE 30-bus system with RES, the total fuel cost is achieved 782.0343 ($/h) using the IEO algorithm, which are better than EO (782.0367 $/h), AEO (782,082 $/h), CSA (782.2001 $/h), TLBO (782.0373 $/h), and ABC (782.0560). Besides, simulation time of the IEO is smaller from 2 to 2.67 compared to EO, AEO, CSA, TLBO, and ABC methods. In addition, the standard deviation value of the IEO (0.0396) is lower than EO (1.4318), AEO (0.1982), CSA (0.1943), TLBO (0.5949), and ABC (0.5875). For large-scale IEEE 118-bus with RES systems, using proposed probability equation (29) in the proposed IEO method decreases the run time higher compared to other methods. The simulation results show that, the proposed IEO algorithm also is one of effective and reliable methods for dealing problem of OPF with RES in large-scale and complex systems.



---

*Source: 7827164-2022-06-30.xml*