Application of Equilibrium Optimizer Algorithm for Optimal Power Flow with High Penetration of Renewable Energy

: In recent decades, the energy market around the world has been reshaped to accommodate the high penetration of renewable energy resources. Although renewable energy sources have brought various beneﬁts, including low operation cost of wind and solar PV power plants, and reducing the environmental risks associated with the conventional power resources, they have imposed a wide range of difﬁculties in power system planning and operation. Naturally, classical optimal power ﬂow (OPF) is a nonlinear problem. Integrating renewable energy resources with conventional thermal power generators escalates the difﬁculty of the OPF problem due to the uncertain and intermittent nature of these resources. To address the complexity associated with the process of the integration of renewable energy resources into the classical electric power systems, two probability distribution functions (Weibull and lognormal) are used to forecast the voltaic power output of wind and solar photovoltaic, respectively. Optimal power ﬂow, including renewable energy, is formulated as a single-objective and multi-objective problem in which many objective functions are considered, such as minimizing the fuel cost, emission, real power loss, and voltage deviation. Real power generation, bus voltage, load tap changers ratios, and shunt compensators values are optimized under various power systems’ constraints. This paper aims to solve the OPF problem and examines the effect of renewable energy resources on the above-mentioned objective functions. A combined model of wind integrated IEEE 30-bus system, solar PV integrated IEEE 30-bus system, and hybrid wind and solar PV integrated IEEE 30-bus system is performed using the equilibrium optimizer technique (EO) and other ﬁve heuristic search methods. A comparison of simulation and statistical results of EO with other optimization techniques showed that EO is more effective and superior and provides the lowest optimization value in term of electric power generation, real power loss, emission index and voltage deviation.


Background
The urgent need for reducing the fuel cost of the conventional power generation units and minimizing the greenhouse gases emitted from the thermal power generators have led various electric power companies to go toward utilizing renewable energy resources. Furthermore, the advanced technologies of renewable energy resources have contributed significantly to them becoming the most 1. 3

. Contribution and Paper Organization
In the present work, an equilibrium optimizer [34], which is a novel optimization method inspired by controlling the volume mass balance model for estimating both equilibrium and dynamic states, is used to prove its performance in solving the OPF problem. It is implemented on (i) IEEE 30-bus system, (ii) wind integrated IEEE 30-bus system, (iii) solar PV integrated IEEE 30-bus system, and (iv) hybrid wind and solar PV integrated IEEE 30-bus system. Real power loss minimizations, total cost minimization of generating units and emission index minimization are considered to be the objective functions of the OPF problem. Weibull and lognormal probability distribution functions are used to model the wind speed and solar irradiance to forecast the output power of wind and solar PV systems. Furthermore, aiming to fill the gap in the literature, this paper investigates the impact of the presence of only wind or only solar PV or both of them on enhancing the objective functions of the OPF problem. In addition, a comprehensive statistical analysis for the equilibrium optimizer technique (EO) and other optimization techniques are used.
The rest of this paper is organized as follows: the formulation of OPF problem is described in Section 2. Then, a mathematical model of wind and solar PV plants is introduced in Section 3. Section 4 presents the equilibrium optimizer technique (EO) and its implementation to solve the OPF problem. Section 5 presents the test systems and the input parameters of the test systems and the optimization methods. Simulation results are explained in Section 6. Finally, Section 7 draws the conclusion of this work.

General Structure of OPF
Generally, OPF aims to minimizes some objective functions. f o is the objective function to be minimized, and h and g are the equality and inequality constraints in the power system network; OPF can be expressed as [14,35]: x is a state vector of dependent variables including the real power of swing generator (P G 1 ), (V L i ) is the voltage magnitude of load buses, (Q G i ) is the reactive power of generator at i th bus and (S l i ) is the loading of the i th transmission line. x can be expressed as follows [14,35]: x = [P G 1 , V L 1 , . . . , V L npq , Q G1 , . . . , Q G NG , S l 1 , . . . , S l n l ] T where npq and n l are the number of PQ buses and transmission lines. S l and n l are the loading of transmission lines and the number of transmission lines, respectively. u is a vector consisting of control variables, (P G i ) is the real power of all generators excluding swing generator, (V G i ) is the voltage magnitude of generators, (TS) is the branch transformer tap, and (Q C ) is the shunt capacitor. u can be expressed as follows [14,35]: u = [P G 2 , . . . , P G NG , V G 1 , . . . , V G NG , Q C1 , . . . , Q CN c , TS 1 , . . . , TS N T ] T where, NG, N c and N T are the number of generators, shunt VAR compensator and transformers, respectively.

Objective Functions of OPF
Here, the first four cases dealt with solving single objective OPF and the last one addressed the multi-objective OPF.

•
Case 1: real power loss minimization Due to the presence of the inherent resistance for the transmission lines, the aim of this function is to minimize the active power losses and it is expressed as [14,35]: where G q(ij) is the conductance of q th transmission line, and V i and V j are the voltage magnitude of terminal buses of transmission line. • Case 2: emission index minimization In the present case, the target is to reduce the harmful gases emission from the thermal generation units, and the coefficients of the gas emission of the thermal power generators are given in Table 1. Emission in tons per hour (t/h) can be described by [14,35]: where α, β, γ, ω and µ are the emission coefficients and G1,G2,G3,G4,G5, and G6 represent thermal power generators at buses 1, 2,5,8,11, and 13, respectively, as given in Table 1. • Case 3: Basic fuel cost minimization The relationship between fuel cost ($/h) and the power generated from the thermal generating units can be approximately given by the quadratic relationship and it is expressed as [14,35]: where a i , b i , c i are the cost coefficient of the thermal generators and these coefficients are provided in Table 2. It also plays a significant role in keeping the voltage quality and security of the electrical power network.This case is expressed as [14,35]: • Case 5: Minimization of basic the fuel cost, emission index, voltage deviation and the real power losses.
The aim of this case is to reduce quadratic fuel cost, active power transmission losses, environmental emission index and voltage deviation index simultaneously. It can be defined as follows [14,35]: where λ p , λ VD and λ E are weight factors and they are assumed to be 22, 21 and 19, respectively as in [14].

Constraints
The constraints of OPF are usually categorized into [14,35]: 1. Equality constraints The equality constraints of OPF are usually represented by the load flow equations: where P D i , Q D i , N B , and θ ik are the active and reactive load demand, the number of buses and the angle difference between bus i and k, respectively. G ik and B ik are the transfer and susceptance conductance.

Inequality constraints
It can be defined by operating limits on the equipment of the power system, transmission loading and voltage of load buses.
(a) Constraints of thermal and renewable energy generating units (b) Constraints of the transformer tap setting (c) constraints of the shunt compensator Q C,j,min ≤ Q C ≤ Q C,j,max j = 1, . . . , N C

Constraint Handling
In order to decline the infeasible solutions of OPF and keep the dependent variables within the allowable ranges, a penalty function was modeled and added to the objective functions defined in Section 2.2 [14,35].
where K Q , K p , K V and K S are the values of penalty factors associated with generation reactive power, generation real power of the swing generator, load bus voltages and line flow of transmission lines. They are assumed to be 100, 100, 100, and 100,000, respectively [14,36], and x Lim is the value of the violated limit of dependent variables(x). It is equal to x max in case of x > x max or x min in case of x < x min .

Uncertain and Power Model of Wind Turbines
The wind speed of the wind turbines follows the Weibull probability distribution function. The characteristic of the output power generated by the wind turbine is a random variable depending on wind speed. The Weibull probability distribution function with dimensionless shape factor (k) and scale factor (c) is used to model the wind speed f v (v). The wind speed ( f v (v)) can be expressed mathematically as [1,2,37,38]: The electrical energy generated by a wind turbine is the result of converting the kinetic energy of wind. The actual output power of wind turbines (P w (v)) can be presented as [1,2,37,38]: where (P wr ), (v in ), (v out ) and (v r ) are the rated power of the wind turbine, the cut-in wind speed of the wind turbine, the cut-out wind speed and the rated wind speed, respectively.

Calculation of Direct, Underestimation and Overestimation Cost of Wind Power
The direct cost of wind power plant can be defined as [4,[39][40][41]: C w,j (P ws,j ) = g j P ws,j where g j is the direct cost coefficient of wind plant. The cost function is overestimated because the actual generated power from the wind turbine is less than the estimated power by mathematical equations. The overestimation cost is used for reverse the requirements when the estimated output power of the wind turbine is more than actual output power. Reserve cost for the j th wind turbine can be defined as [4,[39][40][41]: C Rw,j (P ws,j − P wav,j ) = K Rw,j (P ws,j − P wav,j ) where K Rw,j , P wav,j , P ws,j and f w(P w,j ) are the reserve cost coefficient pertaining to the j th wind turbine, the actual available power from the same plant, the estimated power from the j th wind turbine and the wind power probability density function for j th wind turbine. Underestimation cost function of the wind turbine is due to not using the whole power which is generated from the wind turbine. In other words, when the generated power from the wind turbine is more than the estimated power, the underestimation cost function is applied as a penalty due to wasting the surplus power. The penalty cost for the j th wind turbine can be defined as [4,[39][40][41]: where K Pw,j is a coefficient representing the penalty cost of the j th wind turbine and P wr,j is the rated output power which is generated from the j th wind turbine. As shown in Section 3.1.2, the total cost of wind power turbines (C W T ) can be described as follows: C w,j (P ws,j ) + C Rw,j (P ws,j − P wav,j ) + C Pw,j (P wav,j − P ws,j ) (24) where N w is the number of wind power turbines.

Uncertain and Power Model of Solar PV Plants
Solar irradiance can be modelled by lognormal probability distribution function due to its uncertain and stochastic nature. The lognormal probability distribution is a function of solar irradiance (G) with mean µ and standard deviation σ, it can be expressed mathematically as [3,4]: The main role of PV systems is to convert the solar irradiance to electrical energy. The output power of PV system (P s (G)) as a function of irradiance can be estimated as [4,39]: where G std represents the solar irradiance in standard environment,R c is a certain irradiance point, and P sr is the rated output power which is generated from the solar PV system.

Calculation of Direct, Underestimation, and Overestimation Cost of Solar PV Power
The direct cost of solar power plant can be defined as [4,39]: where h k is a coefficient which represents the direct cost of the solar photovoltaic plant. The same case as in the wind energy system, the solar energy system involves overestimation and underestimation cost due to its uncertain output power. Reserve cost for the overestimation of k th solar PV system is [4,39]: C Rs,k (P ss,k − P sav,k ) = K Rs,k (P ss,k − P sav,k ) = K Rs,k * f s (P sav,k < P ss,k ) * [P ss,k − E(P sav,k < P ss,k )] (28) where K Rs,k is a coefficient which represents the reserve cost pertaining to k th solar PV system, P sav,k is the actual available power from the same plant, f s (P sav , k < P ss,k ) is the probability of solar power shortage occurrence than the scheduled power (P ss,k ) and E(P sav,k < P ss,k ) is the expectation of solar PV power below P ss,k .
In the case of the underestimation of k th solar PV system, the penalty cost is given as [4,39]: where K Ps,k is a coefficient represents the penalty cost pertaining to k th solar PV system, f s (P sav,k < P ss,k ) is the probability of solar power surplus than the scheduled power (P ss,k ) and E(P sav,k < P ss,k ) is the expectation of solar PV power above P ss,k . As explained in Section 3.2.2, the total cost of solar PV plants (C PV T ) consists of three terms(direct, underestimation and overestimation cost) and it can be given as follows [4,39]: C s,k (P ss,k ) + C Ps,k (P sav,k − P ss,k ) + C Rs,k (P ss,k − P sav,k ) (30) where N PV is the number of the solar PV plants.

Inspiration and Mathematical Model
The main inspiration for this algorithm is the dynamic mass balance equation which describes the conservation of mass that enters, leaves or generates in a control volume. This equation is a first-order ordinary differential equation and it is described as the following [34]: where V dC dt is the rate of change of mass in volume, (V), C is the concentration inside the volume(V), V is the control volume, Q is the volumetric flow rate into and out of the control volume, C eq is the concentration at an equilibrium state, and G is the mass generator rate inside the control volume.
After reaching the steady equilibrium state of Equation (31) that is reformulated as a function of Q V , which is called turnover rate λ = Q V . The following equations are derived from Equation (31) to solve for (C) as a function of time (t) [34]: where F is an exponential term to assist EO having a balance between exploitation and exploration, t 0 is the initial start time, and C 0 is the initial concentration. The Equation (35) introduces three rules for updating the concentration of each particle. The equilibrium concentration is the first term which is described as one of the best-so-far solutions randomly chosen from the equilibrium pool. The difference between a concentration of a particle and the equilibrium state is the second term which helps particles to globally explore the domain. The final term is called the generation rate which mainly acts as an exploiter or solution refiner [34].

Initialization and Function Evaluation
Firstly, the optimization process starts with the initial population. The Equation (36) describes the initial concentration process which depends on the number of particles and dimensions that are initialized in the search space in a uniform random manner [34].
where C initial i is the initial concentration vector of the ith particle, C min is the minimum value for the dimensions, C max is the maximum value for the dimensions and rand i is a random vector ranging between zero and one. After that, the fitness function of the particles are evaluated and then solved to determine the equilibrium conditions.

Equilibrium Pool and Candidates (C Eq )
The global optimum of EO is represented by the equilibrium state. At the beginning, no information about the equilibrium state exists, but equilibrium candidates are identified to provide a search domain for the particles. There are five equilibrium candidates as given in Equation (37). Four of them are the best-so-far particles determined during the optimization process and the last one is the arithmetic mean of the previous-mentioned four particles. The main goal of the first four candidates is to improve the exploration capability, whereas the fifth candidate enhances the exploitation [34]

Exponential Term (F)
The exponential term (F) helps EO to have an acceptable balance between exploration and exploitation. Referring back to Equation (34), the time (t) in Equation (34) depends on the iteration (Iter) and it is described as follows [34]:

Iter
Max iter (39) For the purpose of convergence, t 0 in Equation (10) is proposed to slow down the search speed as well as enhancing the exploration and exploitation ability of EO [34].
where a 1 and a 2 are constant values for controlling exploration and exploitation ability, sign( − → r − 0.5) is a factor that determines the direction of exploration and exploitation and r is a random vector ranges between zero and one.

The Generation Rate (G)
The generation rate aims to provide the exact solution by enhancing the exploitation ability of EO and can be described as [34]: After assumption that k = λ, the equation of generation rate was updated as follows [34]: where r 1 and r 2 are a random number between zero and one, GCP is the generation rate control parameter. The generation rate control parameter (GCP) mainly depends on generation probability (GP) which defines the number of particles of the generation term to update their states.
State of the art state that EO at GP = 0.5; EO can achieve a good balance between exploration and exploitation. The updating rule of EO is given as: The second and third terms of Equation (45) can increase variation and thus help EO to better explore in case of they have same signs or to decrease the variation and aiding EO in local searches in case of having opposite signs [34].

Particle's Memory Saving
This can help each particle track with its coordinates in the space. It aids EO in exploitation capability and avoids getting trapped in local minima [34].

Implementation of EO to Solve the OPF Problem
The proposed EO is applied to solve OPF problem including wind and solar PV generation units. The following pseudo code and flowchart shown in Figure 1 explain the steps of the application of EO for OPF problem.
1. Define the control and dependent variables and their limits, as well as the target objective function defined in Section 2.2 [34]. 2. Collect and read the input data of the power system under test, such as data of transmission lines, transformers, shunt VAR compensator, loads and generation units. 3. Calculate the estimated output power of solar PV and wind power generation units, as defined and explained in Section 3 [34]. 4. Initialize the particle's populations [34]. 5. Assign a large number to the fitness of equilibrium candidates and let a1 = 2; a2 = 1; GP = 0.5 [34]. 6. Do the main while loop as the following [34]:

Description of the Test Power Systems
• Test system 1: IEEE 30-bus system The IEEE 30-bus system consists of six thermal power generators, as presented in Figure 2.
The data about transmission lines, tap changing transformers, AVR compensators, limitations on generators and load voltages, active and reactive power demand are given in [42][43][44]. The general specifications of this system are described in Table 3.  Figure 3. The IEEE 30-bus system has been modified by replacing the thermal power generating units at buses 5, 11, and 13 with wind generator at bus 11 and solar PV at buses 5 and 13. In addition, the new solar PV and wind power generators are constructed at bus 24, and 30, respectively. The objective functions defined in Section 2.2 are modified by adding the output power of solar PV plants (P s (G)) and the output power of wind plants (P w (v)) given in Section 3. Case 3 and case 5 described in Section 2.2 are modified by adding the total cost of solar PV plants (C PV T ) and the total cost of wind plants (C w T ) defined in Section 3. The specification of this hybrid power system is given in Table 3. The data of wind and solar PV plants are described in Tables 4 and 5, respectively.   Figure 4. The objective funtions defined in Section 2.2 are modified by adding the output power of solar PV plants (P s (G)) given in Section 3. Case 3 and case 5 described in Section 2.2 are modified by adding the total cost of solar PV plants (C PV T ) defined in Section 3. The general data of this system and solar PV plants are presented in Tables 3 and 6, respectively.  4. Solar PV integrated IEEE 30-bus system.
• Test system 4: wind integrated IEEE 30-bus system In this system, the IEEE 30-bus system is modified by replacing the thermal power generating units at buses 5, 11, and 13 with wind power generators. Moreover, two new wind generators have been added at buses 24, and 30, as seen in Figure 5. The objective functions defined in Section 2.2 is modified by adding the output power of wind plants (P w (v)) given in Section 3. Case 3 and case 5 described in Section 2.2 are modified by adding the total cost of wind plants (C w T ) defined in Section 3. The general specifications of this system and the data of wind power plants are given in Tables 3 and 7, respectively.

Control Parameters of Optimization Methods
The number of iterations, population size, testing ranges and other parameters of the optimization methods are given in Table 8.

Results and Discussion
The performance and effectiveness of the EO are verified for solving OPF problem by carrying out 20 independent test trial runs for all cases discussed in Section 2.2. The EO [34] and other five metaheuristic optimization techniques: MFO [47], TACPSO [45], AGPSO1 [45], TLBO [46] and MPSO [45] have been tested on four power test systems given in Section 5.1. All these optimization techniques are implemented on 2.8-GHz i7 PC with 16 GB of RAM using MATLAB 2017. The EO [34], TLBO [46], MPSO [45], MFO [47], AGPSO1 [45], and TACPSO [45] algorithms are implemented on the test system 1, test system 2, test system 3, and test system 4 for the minimization of the real power loss as defined in Section 2.2. Figure 6a shows the convergence characteristics of real power loss yielded by the best solution of the EO and other optimization methods for test system 1. It observed that the better convergence characteristic is yielded by the EO. Furthermore, Figure 6b,c display voltage and loading profiles of test system 1 for case 1. It is clear that the EO and other optimization methods obey the voltage limits of buses and loading limits of transmission lines. The results of EO and other techniques for test system 1 are displayed in Table 9. It can be observed that EO achieves the minimum real power loss, but other optimization techniques such as TLBO and MFO have less fuel cost at 967.24 $/h and 967.44 $/h, respectively. Furthermore, other techniques have less voltage deviations at minimization this objective function. However, it is clear that the loading of transmission lines for EO is healthy and less than other methods. (c) Loading profile. The loss and loading profiles using EO for all test systems are given in Figure 7. The optimal (best) results yielded by the EO method for the test system 1, test system 2, test system 3, and test system 4 are tabulated in Table 10. From Figure 7 and Table 10, it is seen that the losses of test system 2, test system 3, and test system 4 reduced by 23.6%, 31.52%, and 33.32%, respectively, compared to test system 1. Additionally, it can be seen that the contribution of power generation of Test system 2 from wind, solar PV and thermal power generation are 15.05%, 33.01%, and 51.93%, respectively. With respect to Test system 3, the contribution of power generation from solar PV and thermal power are 53.33% and 46.66%, respectively. Besides, the wind power plants of Test system 4 contributes 52.35% of the total power generation.  The statistical results (the best, the worst, the mean, and the standard deviation) of the real power loss for the EO and other optimization techniques are given in Table 11. As shown in Table 11, the minimum best, standard deviation, and mean are resulted from the EO. As expected, the addition and location of the renewable energy resources in the power system have a significant impact on reducing the real power loss.

Emission Index Minimization
In this case, the emission index defined in Section 2.2 was minimized for all test systems. Figure 8 demonstrates the convergence characteristics, loss profiles, and loading profiles for emission minimization using EO and other methods. It can be noticed from Figure 8a that the EO has the smoothest and speediest convergence curves in comparing with other techniques, as well as Figure 8b,c showing that there is no violation in the voltage limits of buses and loading limits of transmission lines. As we can see from Table 13 and Figure 8c that EO can achieve the lowest real power loss and the lowest loading of the transmission lines while minimizing this objective function, but other optimization methods can obtain less voltage deviations in comparison to EO.
The best (optimal) results obtained using the EO for all test systems for case 2 are shown in Table 12. As we can see from Figure 9 and Table 12 that emission index reduced by 55.54% for test system 2, test system 3, and test system 4 compared to test system 1. In this case, the contribution of power generation from wind power plants for Test system 2 and Test system 4 are 12.25% and 54.12% of the total generation, respectively. In addition, the contribution of power generation from solar PV for Test system 2 and Test system 3 are 41.69% and 53.93% of the total generation, respectively.     Table 13 presents the results of the EO and other methods for test system 1 with the minimization of emission index. For example, the objective function of case 2 for EO was 0.204819 ton/h compared to 0.204862 ton/h and 0.204885 ton/h for MFO [47] and TLBO [46] algorithms, respectively.  Table 14 summarizes the statistical results for the present case. It can be found from Table 14 that the EO provides the smallest best, standard deviation, and median than other methods. Table 14. Summary of the statistical analysis of case 2 for Test system 1.

Best
Worst The comparative convergence characteristics, loading profiles, and loss profiles for test system 1 for the EO and other optimization techniques are presented in Figure 10. As observed in Figure 10, the voltage and loading profiles are kept within the acceptable ranges and the EO gives the best convergence characteristics compared to other methods. The optimal results of the EO and other techniques for test system 1 are summarized in Table 15. From Table 15, the EO leads to 800.4486 $/h total cost of generators which is better than the total cost obtained by the other compared methods. From Figure 10b,c, it can be found that even though EO can obtain the minimum value of the total cost of power generation, the loading of the transmission lines is more than other methods and voltage deviation of EO is higher than other optimization techniques.   The statistical results yielded by the EO and other optimization techniques are given in Table 16. From Table 17 and Figure 11, it can be observed that the total cost of generating units for test system 2, test system 3, and test system 4 declined by 3.54%, 3.47%, and 2.91%, respectively, compared to test system 1. In this case, as shown in Table 17; the contribution of power generation from wind power plants for Test system 2 and Test system 4 are 9.58% and 35.68% of the total generation, respectively. Moreover, the contribution of power generation from solar PV for Test system 2 and Test system 3 are 21.33% and 41.17% of the total generation, respectively.  Figure 12 demonstrates the voltage profiles for all test systems for this case using EO. The optimal solution obtained by EO for test system 1, test system 2, test system 3, and test system 4 are tabulated in Table 18. As shown in Figure 12 and Table 18, the presence of the renewable energy resources improves the voltage profiles and reduced the voltage deviation for test system 2, test system 3, and test system 4 by 22.46%, 37.39%, and 29.61%, respectively, compared to test system 1. Besides, it can be observed that the power generation contribution of Test system 2 from wind, solar PV and thermal power generation are 17.57%, 17.87%, and 64.54%, respectively. With respect to Test system 3, the contribution of power generation from solar PV and thermal power are 34.69% and 65.30%, respectively. Moreover, the wind power plants of Test system 4 contribute 55.34% of the total power generation.  It is clear from Table 19; the minimum best, standard deviation, and median are obtained by the EO.
From Figure 13, the voltage and loading profiles for this case for all optimization methods obey the constraints of voltages at load buses and transmission line loading. It can also be observed that the EO convergence characteristic outperforms the convergence characteristics of other methods. The results of EO and other methods for test system 1 are given in Table 20. From Figure 13b,c and Table 20, it can be seen that EO achieves the minimum emission index while minimizing the objective function of voltage deviation. In addition, EO and MFO can obtain the lowest real power loss at 6.528 MW and 5.965 MW, respectively. Nevertheless, MPSO, TLBO, TACPSO, and ACPSO1 obtain lower fuel cost in comparison to EO.   It is clear from Figure 14 that the EO has the best convergence characteristics compared to the other optimization algorithms and the voltage and loading profiles for all algorithm ranges within the allowable limits. The results of EO and other methods for test system 1 of this case are shown in Table 21. It is clear from Figure 14b,c and Table 21 that while EO achieves the minimum value of voltage deviation and fuel cost in comparison to other methods, other optimization techniques obtain lower values of emission index and real power loss than EO.   The statistical analysis of the EO and other methods for test system 1 is given in Table 22. As shown in the table, the EO gives the minimum best, median and standard deviation. It is clear from Figure 15 and Table 23 that the objective function for this case for test system 2, test system 3, and test system 4 dropped by 3.90%, 7.77%, and 7.84%, respectively compared to test system 1. It is found from Table 23 that the real power loss for test system 2, test system 3, and test system 4 dropped by 30.94%, 20.75%, and 46.06%, respectively compared to test system 1. It can be noted in Figure 15 that the contribution of power generation from wind for Test system 2 and Test system 4 are 15.67% and 49.37% of the total power generation, respectively. While the solar PV contributes 33.36% and 44.69% of the total power for test system 2 and test system 3, respectively. Moreover, it is observed from Table 23 that emission index for test system 2, test system 3, and test system 4 dropped by 61.24%, 54.91%, and 58.58%, respectively compared to test system 1.

Conclusions
In this study, a novel proposed EO method has been successfully applied to solve single and multi-objective OPF with integrated wind turbines and solar PV generators. Its performance and effectiveness were evaluated on four power system, namely: IEEE 30-bus system, wind integrated IEEE 30-bus system, solar PV integrated IEEE 30-bus system, and hybrid wind and solar PV integrated IEEE 30-bus system. Realistic models for the wind turbines and solar PV systems have been proposed and thus real power outputs of wind turbines and solar PV power plants have been accurately forecasted. Therefore, a correct and efficient decision can be taken for inclusion the wind turbines and solar PV power plants in the proper locations. The simulation and statistical results indicate and approve that the EO [34] method outperforms other optimization techniques, namely: TLBO [46], MPSO [45], MFO [47], AGPSO1 [45], and TACPSO [45]. Our research has highlighted the importance of the proper locations of the renewable energy resources on improving the objective functions of OPF problem. Furthermore, adding wind turbines and solar PV play an integral role in enhancing the performance of the standard IEEE 30-bus system. For example, they significantly reduce the fuel cost and emission of the conventional power generators, as well as minimize real power loss and voltage deviation.
Author Contributions: Both authors have made equal contributions to this work. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.