#### 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.

#### 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. Section 2 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.

###### 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 the *j*^{th} wind power plant in terms of scheduled power is modelled as follows: where is coefficient of direct cost and is scheduled power of the *j*^{th} wind power plant. The direct cost of the *k*^{th} solar power plant in terms of scheduled power is modelled as follows: where is the direct cost coefficient and is the scheduled power of the *k*^{th} 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 the *j*^{th} wind plant is where is the incurred cost for the shortfall of wind power output, is the actual capacity of the wind farm that is less than the scheduled capacity, is probability of wind power shortage, is the scheduled wind power output according to , the scheduled value of the left half-plane in Figure 2, and is the coefficient of the incurred cost. Reserve cost for the *k*^{th} solar plant is where is the incurred cost for the shortfall of solar power is the scheduled and actual solar power output, is the probability of a solar power shortage is the scheduled power output of solar below , the scheduled value of the left half-plane in Figure 3, and is the coefficient of the incurred costs.

###### 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 of *j*^{th} wind plant is where is the penalty cost of wind power surplus, is the probability of a wind power surplus, is the scheduled of power wind output under , the scheduled value of right half-plane in Figure 2 (in kW), and is penalty cost of wind power surplus. Penalty cost for the underestimation of *k*^{th} solar plant is where is the penalty cost of the solar surplus, is the probability of a solar surplus; is the scheduled of solar power output under , the scheduled value of right half-plane in Figure 3 (in kW), and is penalty cost of solar surplus. The total cost of wind power generated from a wind farm and solar power is described as follows:

##### 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

where is cost function of thermal power generators. 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:

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):where , and are fuel cost coefficients of the *i*^{th} generator, while represents the total number of thermal power plants.

###### 2.2.2. Total Emission

The total emissions *(ton/h)* of thermal power plants are given bywhere , and are emission factors of *i*^{th} thermal power plant.

###### 2.2.3. Total Transmission Loss

The total power loss in the transmission network can be calculated according to the following formula:

##### 2.3. Constraints

###### 2.3.1. Equality Constraints

Constraints on real and reactive power balance:###### 2.3.2. Inequality Constraints

The limits of power generation: The limits of generator voltage bus and load voltage bus: The limits of transmission line: The limits of switchable capacitor capacity: The limits of transformer tap:#### 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 (m/s) is given by

The average wind speed is determined bywhere is gamma function and is described as follows:

Wind speed following the Weibull distribution can be calculated as follows:where *r* is a random number uniformly distributed in the range [0,1].

Figures 4 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.

Relationship between power output of wind turbine and wind speed () is described aswhere , , and 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 *P*_{wr}.

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:

The average solar irradiance is determined by

Figure 6 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 bywhere is the solar irradiance in standard environment, is a certain irradiance point, and is the rated output power of the solar PV unit. In this study, and are set to 800 (W/m^{2}) and 120 (W/m^{2}), respectively [41].

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 Figure 3.

##### 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: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, is the probability that the RES power is less than or greater than the scheduled power value *K*, and 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: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 *X*_{eq}, 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 of *N*_{pop} samples as follows:where and rand [0, 1] is a random number between 0 and 1.

###### 4.1.2. Determine the Equilibrium Group

The four candidates in the 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:where is the best, second best, third best, and fourth best solution in the population; is the average candidate (average solution) of the four best individuals and obtained by the following equation:

###### 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:

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]:

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,where is the initial value defined as follows:in which, 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:

*X*_{eq} 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 term *E* 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:

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:where rand is a vector of size [1 dim] containing random values distributed over the interval [0, 1], and and are random individuals in the equilibrium pool population of *k*^{th} 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:

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 size *N*_{pop}, maximum number of iterations Max-It, and number of simulation samples Monte-Carlo.

*Step 3. *Calculate *N*_{Montecarlo} 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 *P*_{Wind} and *P*_{Solar}, 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:where are the penalty coefficients for the real power at the slack bus, reactive power at all generator buses , voltage magnitude at all load bus , and apparent power transmitted on branches, respectively. 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 solution *X*_{Best} *=* *X*eq1.

*Step 9. *Update the new candidates as follows:

*Step 10. *Update the candidate limit according to

*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 Satisfied

*Step 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.

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

The data of the modified IEEE 30-bus system are summarized in Table 3. 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: where the *IF* 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: 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. *N*_{b} 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 *N*_{b} in the total number of runs, from which the *TNFs* value will be calculated.

*Case1. *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 250^{th} 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 250^{th} 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.

*Case 2. *Generation emission minimization

The best, worst, and average value, and standard deviation after 50 independent runs are shown as Table 9. 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.

*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 Table 11. 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.

##### 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).

*Case 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 Table 16. 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).

*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 Table 16. 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.

#### 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.

#### Nomenclature

: | Total generation cost ($/h) of thermal units |

: | Number of thermal power plants, wind power plants, and solar power plants |

: | The cost of wind power |

: | Number of simulations Monte-Carlo |

: | The cost of solar power |

: | Power at the i^{th} generator node |

: | Total emission |

: | Power at the j^{th} load node |

: | Power losses |

VG, VL: | Voltage magnitude at the generator and load buses |

: | Direct cost coefficients of the j^{th} wind power and the k^{th} solar power plant |

: | Minimum output power for power i^{th} plant |

: | Reserve cost coefficients of j^{th} wind power plant and k^{th} solar power plant |

: | Maximum output power for power i^{th} plant |

: | Penalty cost coefficients of j^{th} wind power plant and k^{th} solar power plant |

: | Voltage magnitude limit at generator nodes |

: | Real power and reactive power of the ith thermal power plant |

: | Voltage magnitude limit at load nodes |

: | Scheduled power of the j^{th} wind power plant |

: | Transmission capacity and line limit |

: | Scheduled power of the k^{th} solar power plant wind speed (the wind speed random variable) (m/s) |

: | Number of power plants, load nodes, and transmission lines cut-in wind speed (m/s) |

: | Number of transformer tap and switchable capacitor capacity rated wind speed (m/s) |

: | Voltage phase angle of node i^{th} and node j^{th} cut-out wind speed (m/s) |

L: | The angle of the bus conduction matrix Y with two nodes i^{th} and j^{th}. |

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

We acknowledge Ho Chi Minh City University of Technology (HCMUT), VNU-HCM, for the support of time and facilities for this study.