Optimal Operation of Conventional Power Generation with High Penetration of Renewable Energy using Equilibrium Optimizer Technique

21 1 Protection and Metering Department, National Electric Power Company, Amman 11181,Jordan; khalednusair2016@yahoo.com 2 Department of Power Engineering, Hijjawi Faculty for Engineering Technology, Yarmouk University, Irbid 21163,Jordan;lina.hmoud@yu.edu.jo * Correspondence: lina.hmoud@yu.edu.jo Abstract: Over the last decades, the energy market around the world has reshaped due to accommodating the high penetration of renewable energy resources. Although renewable energy sources have brought various benefits, 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 difficulties in power system planning and operation. Naturally, classical optimal power flow (OPF) is a nonlinear problem. Integrating renewable energy resources with conventional thermal power generators escalates the difficulty 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 flow, 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 30bus system, solar PV integrated IEEE 30-bus system, and hybrid wind and solar PV integrated IEEE 30-bus system are performed using the equilibrium optimizer technique (EO) and other five heuristic search methods. A comparison of simulation and statistical results of EO with other optimization techniques showed that EO is more effective and superior.


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, advanced technologies of renewable energy resources have contributed significantly to be the most inexpensive and environmentally friendly. Integrating wind and solar PV in proper locations and appropriate settings of the variables of the conventional power networks may have a significant impact on the performance of power system control and operation.
To make the modeling of wind and solar PV more accurate and realistic, the Weibull probability distribution function was used to forecast the wind speed [1]- [2]. Whereas lognormal probability distribution function is used to mimic the intermittent nature of solar irradiance in [3] [4].

Literature review
Numerous publications in the literature studied the optimal power flow (OPF) problem for systems consisting of conventional power generation and renewable-energy power plants. Deterministic, stochastic or hybrid optimization methods are used extensively to address the issues associated with increased penetration of non-dispatchable renewable energy, advanced controls such as FACTs devices and deregulated electricity markets.
Various conventional optimization techniques are used to solve the OPF problem. For instance, continuous nonlinear programming was proposed [5]. An extended conic quadratic format [6] is presented to solve the economic dispatch and decrease real power loss. Besides, the predictor-corrector interior point algorithm is proposed to fit the OPF for solving nonlinear programming problems [7]. Quadratic programming is used to derive a loss formula based on the incremental power flow [8]. Sequential quadratic programming is used to address large scale OPF; it also depends on transforming the original problem to a sequence of linearly constrained sub-problem by applying an augmented lagrangian [9]. Mixed-integer linear programming to minimize transmission losses and reactive generator outputs are adapted [10]. Although these methods have excellent convergence characteristics, they have various drawbacks, including failing to find the global solution because of non-convexity and facing difficulty while handling the problems with non-differentiable and discontinuous objective functions.
Recently, metaheuristic optimization algorithms have been gaining much attention due to flexibility, free of derivation, and local optima avoidance. Thus, single and multi-objective optimization methods overcome the shortcomings attributed to classical techniques. A gravitational search algorithm to find the optimal solution for OPF and IEEE 30-bus and 57-bus systems are examined [11]. The basic fuel cost, voltage profile, voltage stability, and non-smooth quadratic cost are minimized and optimized using a differential evolution algorithm [12]. The Black hole-based optimization method is used to address the OPF problem for IEEE 30-bus and Algerian 59-bus power systems [13]. Constrained OPF problem for IEEE 30-bus, 57-bus, and 118-bus is optimized using a moth swarm algorithm [14]. A multi-objective OPF to minimize the generation cost and environmental pollution using a fuzzy membership function to choose a compromise solution from the Pareto optimal solutions is discussed [15]. The fuel cost, voltage deviation, and real power loss are minimized as a multi-objective OPF problem using a gravitational search algorithm [16]. A modified teaching learning-based optimization algorithm added a self-adapting wavelet mutation strategy and a fuzzy clustering [17]. A hybrid of fuzzy evolutionary and swarm optimization is proposed to minimize the cost of active power generation and real power losses [18].
A fuzzy-based modified bee colony is presented to solve discrete OPF using multi-objective mixed integer nonlinear [19]. Emission, real power losses, and voltage deviation are all minimized as a multi-objective OPF using a multi-objective modified imperialist competitive algorithm [20]. The particle swarm optimization and the shuffled frog leaping algorithm are hybridized to solve OPF using the generator's constraints such as prohibited zones and valve point effect [21]. A chaotic invasive weed optimization algorithm is proposed to solve the OPF problem with non-smooth and non-convex fuel cost curves [22]. Brainstorming optimization and teaching-learning optimization are hybridized to minimize the fuel cost of thermal generation units [23]. A hybrid optimization algorithm is based on sequential quadratic programming to generate an initial population. Then a differential evolution took that population to find the optimal solution more effectively and it was used to minimize the fuel cost with valve point and the transmission line real losses [24].
A growing and considerable effort have been made in recent years to solve and model the OPF problem, including renewable energy sources. The OPF problem with taking into account uncertainties in the wind, solar, and load forecast and optimized using a genetic algorithm and two-point estimate method [25]. A hybrid method called moth swarm algorithm and gravitational search algorithm is used to solve the problem of OPF, including wind power [26]. A modified two-point estimation method is used to solve probabilistic OPF incorporating wind and solar photovoltaic [27]. Hybrid wind photovoltaic power systems are optimized using the unscented transformation method, which can carry out probabilistic OPF with high accuracy and less computational time 3 of 34 [28]. The OPF, including wind is optimized using a fuzzy-based particle swarm optimization. A fuzzy set modeled the forecast load demand and wind speed [29].
Besides, OPF incorporating wind power energy is optimized by a hybrid algorithm called a hybrid dragonfly with aging particle swarm optimization [30]. Adaptive differential evolution with proper constraint handling method is addressed OPF, including wind and solar. The forecast wind and solar photovoltaic are modeled using Weibull and lognormal probability distribution functions [31]. An optimal reactive power dispatch with solar photovoltaic power and its impact on minimizing real power losses is addressed using the Jaya algorithm to solve this issue [32]. A constrained multi-objective population external optimization method in [33] is presented to minimize the fuel cost and emission in the presence of renewable energy sources. A grey wolf optimization algorithm in [34] was proposed to tune the parameters of a thyristor controlled series compensator and address OPF, including wind and solar power. A gbest guided artificial bee colony optimization in [1] was to find the optimal setting of conventional and renewable power generation.

Contribution and paper organization
In the present work, an equilibrium optimizer [35], 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 analysed.
The rest of this paper is organized as follows: the formulation of OPF problem is described in Section 2. Then, a mathematical models of wind and solar PV plants are introduced in Section 3. Section 4 presents the equilibrium optimizer technique (EO) and its implementation to solve the OPF problem. 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,36] : x is a state vector of dependent variables including 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,36]: where npq, and n l are the number of PQ buses and transmission lines. S l and n l are 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, (T S) is the branch transformer tap, and (Q C ) is the shunt capacitors. u can be expressed as follows [14,36]: 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,36]: Where G q(i j) 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,36]: where α, β , γ, ω and µ are the emission coefficient and they are 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,36]: where a i , b i , c i are the cost coefficient of the thermal generators and these coefficients are provided in Table 2.  • Case 4: Voltage deviation minimization The voltage deviation index is the cumulative deviation of all load buses from nominal value of unity. It also play a significant role in keeping the voltage quality and security of the electrical power network.This case is expressed as [14,36]: • 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,36]: where λ p , λ V D 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,36]: 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 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 T S k,min ≤ T S k ≤ T S k,max k = 1, . . . , N T (c) constraints of the shunt compensator Q C, j,min ≤ Q C ≤ Q C, j,max j = 1, . . . , N C (15) (d) Constraints of the voltages at load buses (e) Constraints of the transmission line loading

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,36].
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,37], 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 . 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,38,39]: Mean of Weibull distribution (M wbl ) can be expressed as [31,[40][41][42]: The electrical energy generated by a wind turbine (P w (v)) is the result of converting of the kinetic energy of wind and it can be estimated as [1,2,38,39]: where (P wt ), (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 [31,[40][41][42]: C w, j (P ws, j ) = g j P ws, j (22) 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 [31,[40][41][42]: where K Rw, j , P wav, j , P ws, j and f w(pw, j) are the reserve cost coefficient pertaining to 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, underestimation cost function is applied as a penalty due to waste the surplus power. The Penalty cost for the j th wind turbine can be defined as [31,[40][41][42]: where K Pw, j is a coefficient represent the penalty cost for 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 ) (25) 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],and [4]: Mean of lognormal distribution M lgn can be expressed as: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 September 2020 doi:10.20944/preprints202009.0344.v1 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 [31,40]: 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 [31,40] : where h k is a coefficient represents the direct cost of solar photovoltaic plant. The same case as in wind energy system, solar energy system in involves overestimation and underestimation cost due its uncertain output power. Reserve cost for the overestimation ofk th solar PV system is [31,40]: where K Rs,k is a coefficient 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 case of the underestimation of k th solar PV system, the penalty cost is given as [31,40]: 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 follow [31,40]: 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 which enters, leaves or generates in a control volume. This equation is a first-order ordinary differential equation and it is described as following [35]: 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.
After reaching the steady equilibrium state of equation (33) that is reformulated as a function of Q V which is called turnover rate λ = Q V . The following equations are derived from equation (33) to solve for (C) as a function of time (t) [35]: dC where F ,t 0 is the initial start time and C 0 is the initial concentration.
The equation (37) 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 [35].

Initialization and function evaluation
Firstly, the optimization process starts with the initial population. The equation (38) describes the initial concentration process which depends on the number of particles and dimensions that initialized in the search space in a uniform random manner [35].
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 ranges 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 is existed, but equilibrium candidates are identified to provide a search domain for the particles. There are five equilibrium candidates as given in equation 39. 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 [35] Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 September 2020 doi:10.20944/preprints202009.0344.v1

Exponential term (F)
The exponential term (F) helps EO to have an acceptable balance between exploration and exploitation. Referring back to equation (36), the time (t) in equation (36) depends on the iteration (Iter) and it is described as follows [35] : Iter Max iter (41) 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 [35].
where a 1 and a 2 are constant values for controlling exploration and exploitation ability, sign( − → r − 0.5) is a factor that determine 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 [35]: After assumption that k = λ ,the equation of generation rate was updated as follows [35]: where r 1 and r 2 are 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 which uses 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 (47) can increase variation and thus helps 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 [35].

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 [35].

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 explains 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. 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. 4. Initialize the particle's populations. 5. Assign a large number to the fitness of equilibrium candidates and let a1=2; a2=1; GP=0.5. 6. Do the main while loop as following [35]: (a) While (current iteration (Iter) <maximum number of iteration (Max-iter)) (b) For i=1: particles' number (n) (c) Find the fitness value of the i th particle i. If fitness (C i ) <fitness (C eq1 ) then Substitute (C eq1 ) with (C i ) and fitness (C eq1 ) with fitness then Substitute (C eq4 ) with (C i ) and fitness (C eq4 ) with fitness (C i ) (d) End (if ) (e) End(for) 7. Find the − − → C avg according to equation (39). 8. Construct the equilibrium pool according to equation (40) [35]. 9. In case of the current iteration >1, accomplish memory saving [35]. 10. Assign t according to equation (41). 11. Do the second for loop as following: For i=1: particles' number (a) Select one candidate from the equilibrium pool randomly. (b) Create the two random vector (λ and r).
(c) Construct F, GCP, G 0 and G according to the equation (36), equation (46), equation (45) and equation (44), respectively [35]. (d) Update the concentration C according to equation (47)/ End the second for loop. 12. Increase the current iteration by one. 13. End the main while loop. 14. Extract and analysis of the results.

Results and Discussion
The performance and effectiveness of the EO is verified for solving OPF problem by carrying out 20 independent test trial runs for all cases discussed in Section 2.2. The EO [35] and other five metaheuristic optimization techniques: MFO [43], TACPSO [44], AGPSO1 [44], TLBO [45] and MPSO [44] have been tested on four power test systems given in Section 6.1 . All these optimization techniques are implemented on 2.8-GHz i7 PC with 16 GB of RAM using MATLAB 2017. The number of iteration, population size ,testing ranges and other parameters of the optimization methods are given in Table 3.

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 1.The data about transmission lines, tap changing transformers, AVR compensators, limitations on generators and load voltages, active and reactive power demand are given in [46][47][48]. The general specifications of this system are described in Table 4.  • Test system 2: 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, and13 with wind power generators. Moreover, new two wind generators have been added at buses 24, and 30, as seen in Figure 2 .The objective functions defined in Section2.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 Table 4 and Table 5, respectively.  Figure 2. Wind integrated IEEE 30-bus system.
• Test system 3: Solar PV integrated IEEE 30-bus system This system is modified by locating solar PV generators at buses 5, 11, and13 instead of the thermal power generators. Furthermore, two new solar power generation units are installed at buses 24, and 30, as shown in Figure 3.The objective funtions defined in Section2.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 Table4 and Table 6, respectively.  Figure 3. Solar PV integrated IEEE 30-bus system.
• Test system 4: Hybrid wind and solar PV integrated IEEE 30-bus system This modified test power system is simulated to show its behavior in the presence of both wind and solar PV power generating units, as depicted in Figure 4. 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 Section2.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 4. The data of wind and solar PV plants are described in Table 7, and Table 8, respectively.   The EO [35], TLBO [45], MPSO [44], MFO [43], AGPSO1 [44], and TACPSO [44] 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. Fig. 5a shows the convergence characteristics of real power loss yielded by the best solution of the EO and other optimization methods for test system1. It observed that the better convergence characteristic is yielded by the EO. Furthermore, Fig. 5b and Fig.5c   The loss and loading profiles using EO for all test systems are given in Fig. 6.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 Fig.6 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.   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. Fig. 7 demonstrates the convergence characteristics, loss profiles, and loading profiles for emission minimization using EO and other methods. It can be noticed from Fig. 7a that the EO has the smoothest and speediest convergence curves in comparing with other techniques, as well as Fig.7b and Fig.7c show that there is no violation in the voltage limits of buses and loading limits of transmission lines.  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 Fig. 8 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.   Table 13 that the EO provides the smallest best, standard deviation, and median than other methods.  Table 14 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 [43] and TLBO [45] algorithms, respectively. The comparative convergence characteristics, loading profiles, and loss profiles for test system 1 for the EO and other optimization techniques are presented in Fig.9. As observed in Fig.9, 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$/hr total cost of generators which is better than the total cost obtained by the other compared methods.
(c) Loading Profile Figure 9. Comparative convergence, voltage and loading Profiles for case 3 for all test systems. From Table 17 and Fig.10, 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.   11 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 Fig.11 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 system1.   Fig.12, 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 be also 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.    Fig.13, the EO has the best convergence characteristics compared to the other optimization algorithms and the voltage and loading profiles for all algorithms ranges within the allowable limits. The results of EO and other methods for test system 1 of this case are shown in Table 21.   The statistical analysis of the EO and other methods for test system 1 are given in Table 22. As shown in table, the EO gives the minimum best, median and standard deviation. It is clear from Fig.14 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 system1. 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 system1.   and approve that the EO [35] method outperforms other optimization techniques ,namely: TLBO [45], MPSO [44], MFO [43], AGPSO1 [44], and TACPSO [44] .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 .