An Optimization-Based Strategy for Solving Optimal Power Flow Problems in a Power System Integrated with Stochastic Solar and Wind Power Energy

: In an effort to reduce greenhouse gas emissions, experts are looking to substitute fossil fuel energy with renewable energy for environmentally sustainable and emission free societies. This paper presents the hybridization of particle swarm optimization (PSO) with grey wolf optimization (GWO), namely a hybrid PSO-GWO algorithm for the solution of optimal power ﬂow (OPF) problems integrated with stochastic solar photovoltaics (SPV) and wind turbines (WT) to enhance global search capabilities towards an optimal solution. A solution approach is used in which SPV and WT output powers are estimated using lognormal and Weibull probability distribution functions respectively, after simulation of 8000 Monte Carlo scenarios. The control variables include the forecast real power generation of SPV and WT, real power of thermal generators except slack-bus, and voltages of all voltage generation buses. The total generation cost of the system is considered the main objective function to be optimized, including the penalty and reserve cost for underestimation and overestimation of SPV and WT, respectively. The proposed solution approach for OPF problems is veriﬁed on the modiﬁed IEEE 30 bus test system. The performance and robustness of the proposed hybrid PSO-GWO algorithm in solving the OPF problem is assessed by comparing the results with ﬁve other metaheuristic optimization algorithms for the same test system, under the same control variables and system constraints. Simulation results conﬁrm that the hybrid PSO-GWO algorithm performs well compared to other algorithms and shows that it can be an efﬁcient choice for the solution of OPF problems.


Introduction
The rapid increase in energy demand and higher prices of conventional thermal generators led the world to a challenging problem of greenhouse gas emissions and an unstable economy. In this regard, energy production sectors are trying to integrate environmentally friendly and cost-effective renewable energy sources (RES) into the grid [1]. The integration of RES provides economic, social and environmental benefits for a sustainable society. Nowadays solar photovoltaics (SPV) and wind turbines (WT) are the two most integrated RES into the grid [2]. The incorporation of RES into the system is a tough and complicated task due to the sporadic nature of solar photovoltaics (SPV) and WT [3]. Due to the intermittent nature of RES, such as SPV and WT sources, it is hard to transmit power over a large distance between generation and consumers. This hurdle can be overcome by integrating enormous storage systems into the transmission network that will also help to reduces fluctuations in the network [4]. Energy storage devices, such as those that are lithium and sodium based [5,6], can be integrated along with conventional thermal and RES. Lithium-ion batteries have gained a lot of significance due to their high energy density, long life span and low self-discharge property, which can be added in the future expansion of microgrids [7]. OPF solutions provide the best operating values of control parameters, which are generators of real power and voltages, to reduce the total generation cost of the system [8]. Optimizing the generation cost system constraints, which are power flow balance, bus voltages, line capacity, and generator capability, is to be satisfied. OPF solutions must incorporate the uncertainties of these RES, which makes OPF a more complex and nonlinear problem.
For OPF problem formulation and solutions, many optimization algorithms and techniques have been developed and used in the past couple of years. These methods can be classified into classical methods, evolutionary based methods, and advanced metaheuristic algorithm-based methods. The classical method consists of linear programming [9], quadratic programming [10], non-linear programming [11], interior point method [12], and dynamic programming etc. For large scale problems, these classical methods are inadequate due to their complex and extensive calculation and their high dependency on initial guess values.
Evolutionary based algorithms copy the evolutionary idea of generating next generations from the parentages by mating to obtain the optimal solution. One of the most popular evolutionary based algorithms is the Genetic Algorithm, which has been used to solve OPF as mentioned in Dashtdar et al. [13], Moeini-Aghtaie et al. [14], and Bakirtzis et al. [15]. Another algorithm is the Differential Evolution (DE), that falls under this category and has been successfully implemented for constraint handling while optimizing the different objectives of OPF [16].
To overcome the discrepancies of the classical and evolutionary based methods, various metaheuristic algorithms have been effectively utilized for OPF solutions. Examples of various metaheuristic algorithms are particle swarm optimization (PSO), the moth-swarm algorithm (MSA), whale optimization algorithm (WOA), grey wolf optimization (GWO), and elephant herding optimization (EHO) etc. Physics based algorithms that are stimulated from different nature/physics laws and swarm-based algorithms which are inspired from the social manner of a cluster of animals have been included in the group of metaheuristic algorithms. The Gravitational Search Algorithm (GSA) has been suggested by Rashedi et al. in [17], that imitates the concept of Newton's gravity law. Mirjalili et al.,in [18], developed a Multi-Verse Optimization (MVO) algorithm which mimics the idea of cosmology that is comprised of wormholes, white holes and black holes. Physics based algorithms that have been utilized for the solution of OPF problems are GSA, presented by researchers in [19,20]. The Electromagnetism-like Mechanism Method (EMM) presented by the authors in [21] and Pinheiro et al. in [22] used the Doubly Modified logarithmic Barrier Function (M2BF) algorithm. Other swarm-based algorithms that have been utilized for OPF solutions exist such as PSO, presented by Kennedy and Eberhart in [23], Moth-flame Optimization (MFO) [24], Ant Lion Optimization (ALO) [25], GWO [26], and the Grasshopper Optimization Algorithm (GOA) [27]. The authors in [28][29][30] used PSO and its enhanced variants to provide an optimal solution for OPF problems. A novel meta-heuristic moth swarm algorithm (MSA) has been used for OPF solutions, which is inspired by the path finding processes of moths [31]. The performance of the moth swarm algorithm has been improved with arithmetic crossover (MSA-AC) and has been successfully utilized in the optimization of OPF problems for HVDC systems [32]. The results obtained show the better exploration and speedy convergence capability of MSA-AC for a global optimum. A modified sine cosine algorithm (MSCA) has been used for OPF solutions by the researcher in [33]. The effectiveness of modified MSCA has been verified in terms of a decrease in computational time and improvement of finding global OPF solutions with different objective functions. In reference [34], a novel improved artificial bee colony (IABC) optimization algorithm has been used for OPF solutions. The poor exploitation capability of the original ABC technique has been improved by an orthogonal learning approach, which guesses the optimum solution among the solution vectors created on restricted trials rather than a complete search. The dominance of IABC has been verified on IEEE 30 and 118 bus test systems, while optimizing the OPF solution, where IABC achieves a lower cost in comparison to other heuristics techniques for all test cases.
To improve the global searching ability, hybridization of two metaheuristic algorithms has been conducted to obtain the global optimum results; a hybrid of dragon fly and PSO (DA-PSO) have been proposed in reference [35] in which the exploitation and exploration capability of both algorithms are enhanced. The competence of the proposed algorithm (DA-PSO) has been confirmed on IEEE 30 and 57 bus test systems, while taking into account multi-objective and single-objective cases of power losses, cost and emissions. The dominance of the suggested algorithm has been proven by simulation results. However, computation time is increased due to the sequential operation of the hybrid technique. The authors in [36] present a hybrid DE along with a harmony search algorithm (Hybrid DE-HS) for OPF solutions. The usefulness of the applied hybrid DE-HS algorithm has been verified on IEEE 30, 118 and 300 bus test systems with reactive power compensation and transformer tap settings as a control variable. Other hybrid algorithms utilized for OPF solutions are hybrid PSO-GSA [37], imperialist competitive algorithm (ICA) hybridized with sequential quadratic programming by Ben Hmida et al. in [38], hybrid of WOA with Pattern Search Algorithm (PS) presented by authors in [39], hybrid of PSO and Artificial Physics Optimization (HPSOAPO) presented by Teeparthi et al. in [40], hybrid of PSO, and Shuffle Frog Leaping algorithm (HPSO-SFLA) presented by Narimani et al. in [41]. Other techniques that have been used with different multi-objective functions for OPF solutions are fuzzy logic-based techniques. A fuzzy logic-based approach has been used for the contingency constrained OPF problem formulation and solution in a decomposed form [42]. The authors in [43] used fuzzy based membership functions to form fuzzy optimal power problem formulations with an artificial bee colony (ABC) algorithm. In [44], the authors proposed an OPF solution using fuzzy logic with integrated GA and PSO. A hybrid of fuzzy based PSO techniques has also been used in [45], while considering the uncertainties in load demand and wind speed. The uncertainties in load demand and wind speed have been forecasted using fuzzy set notation.
Generally, there is not a single algorithm that can be considered better than others to solve all of the objectives of OPF problems, due to the inconsistent nature of OPF solutions. Thus, the performance of an algorithm or technique depends upon the nature and space of the problem that needs to be considered. Therefore, there must exist a space for new algorithms which can efficiently solve various objectives of OPF solutions. The objective of this paper is to accomplish a renewable integrated energy system with the use of WT and SPV sources to fulfill the energy requirements in the best possible way with a reduction in emissions for a sustainable environment. In this paper, utilization of the HPSO-GWO algorithms for OPF solutions incorporated with uncertain RES is presented. The results achieved with the suggested hybrid PSO-GWO algorithm (here after will be called HPSO-GWO) will be compared with PSO, GWO, success history based adaptive DE combined with superiority of feasible technique (SHADE-SF) [8], and with those results of GOA and Barnacles Mating Optimizer (BMO) reported in the literature [46]. The rest of the sections in this paper are set out as follows: Section 2 provides OPF problem mathematical formulations. Section 3 formulates frequency distributions of stochastic SPV and WT. Section 4 presents details of selected metaheuristic algorithms. Section 5 provides a detail discussion on the simulation results of each case. Section 6 provides a conclusion and future work.

Mathematical Formulation of OPF Problem
The standard IEEE 30 bus system consists of twenty-four load buses and six generator buses to which thermal generators are connected. In the proposed modified system, three thermal generators are replaced by two WT and a SPV source to formulate OPF, which is an optimization problem with inequality and equality constraints. The generation at a bus should not violate the line capacity constraints while replacing RES; in this study, the same method of replacing thermal generators with RES has been adopted as in in reference [8,46]. To forecast the output power of SPV and WT sources, lognormal and Weibull probability density functions (PDF) have been used, respectively. The application of the proposed HPSO-GWO algorithm is tested on the modified IEEE 30-bus system (hereafter test system) while considering different objective functions such as total generation cost, power losses, emissions, and cumulative voltage deviation in the system. The control parameters and other details of the test system under study are summarized in Table 1. Along with thermal generators, the adopted test system also includes SPV and WT, which are presented in Figure 1. The uncertain output power of wind and SPV sources must be compensated by the output power of thermal generators or reserve sources. The total generation cost of the entire system represents a combination of costs of all thermal generators, SPV and WT, plus reserve and penalty costs for overestimation and underestimation of wind and SPV sources, respectively. The subsequent sections explain cost function formulations for each type of source.

Cost Function Formulation for Thermal Generators
For the operation of thermal generators, they require fossil fuels and their output power depends on the consumption of fossil fuels. The fuel cost ($/h) and output power (MW) of the thermal generators have a relationship which can be expressed as below: where C TG represents the total cost of n thermal generators producing P TG (MW) real power, while a i , b i and c i represent fuel cost coefficients for ith thermal generators whose values are mentioned in Table 2. There is a greater divergence in the cost function of thermal generators with multi valve steam turbines. Therefore, for modeling a precise and more accurate cost function for thermal generators, valve point effects must be considered. Valve point effect is a sinusoidal function and must be included in the basic cost of thermal generators as an absolute value. The equation for the cost function, after counting the valve point effect, now becomes: where e i and d i represent the valve point effect coefficients of steam turbines. All cost and valve point coefficients used in Equations (1) and (2) are taken from reference [8] and are summarized in Table 2 along with emission coefficients.

Cost Function Formulation for Emission and Carbon Tax
During the operation of thermal generators, harmful emissions are released into the atmosphere due to consumption of fossil fuels. The release of these emissions increases with the escalation of output power from thermal generators. The expression for calculating the emission in ton/h is given below: Emission coefficient values for thermal generators are given in Table 2 and are the same as in reference [8]. To control greenhouse gas emission, various countries are now putting immense pressure on power production sectors due the greenhouse effect [48]. As a result, the world has seen a shift towards solar photovoltaics, due to its tremendous potential across the globe and green renewable nature [49]. To reduce harmful emissions and encourage the integration of environmentally friendly and cost-effective RES such as SPV and wind, carbon tax is applied on a per unit quantity of produced emissions. The emission cost ($/h) for thermal generators is formulated as below:

Cost Function Formulation for SPV and Wind Turbine
RES involves no running costs, unlike fossil fuel based thermal generation plants. RES per unit price is determined by either reimbursement cost to the initial expenditure or renovation and maintenance cost in case of possession by independent systems operators (ISO). RES may be operated and owned by private sectors; in such cases, ISO must pay a per unit price proportionate to the scheduled output energy of these plants. Direct cost of jth wind energy plants is mathematically modeled as a function of scheduled output power: where P w−sch,j represents the schedule output power for jth wind turbine and g w,j represents the direct cost coefficient, linked with the same wind turbine.
Similarly, a direct cost assigned to kth SPV sources producing P s−sch,k scheduled power is: where h k represents the direct cost coefficient of the kth SPV source.

Mathematical Modeling of Uncertainties in RES
There might be a situation in which the RES real power is below the forecasted power, which is referred to as over estimation. The system must be capable of having a reserve source to deliver the continuous power supply to consumers. The cost assigned to generating units of the reserve source, supplied to fulfill the deficit of overestimated units, is named the reserve cost.
Reserve cost for the overestimated units can be expressed as: where P wsch,j represents the scheduled power output from jth WT power and P wac,j represents the actual power output from jth WT. K Rw,j represents the reserve cost coefficient of jth WT. f (p w,j ) represents PDF of WT, where details about probabilities calculation for different wind speeds can be obtained in Biswas et al. [8].
The reserve cost for the kth SPV source may be expressed as below: C Rs,k = K Rs,k (P ssch,k − P sac,k ) = K Rs,k * f s (P sac,k < P ssch,k ) * [P ssch,k − E(P sac,k < P ssch,k )] where P ssch,k represents the scheduled output power from kth SPV source, P sac,k represents actual power output from kth SPV source and K Rs,k represents the reserve cost coefficient of the kth SPV source. f s (P sac,k < P ssch,k ) represents the SPV source deficient amount of power than the scheduled SPV source power P ssch,k . E(P sac,k < P ssch,k ) represents the probability of SPV power below P ssch,k . On other hand, there might be a situation in which the estimated output power from SPV and WT is more than the estimated power. Under such conditions, the output power from RES is underestimated. The excess amount of energy is expected to be useless if not stored or used properly. ISO must pay a penalty cost for the additional amount of energy.
Penalty cost for the additional amount of energy for jth WT may be expressed in the form: C Pw,j = K Pw,j (P wac,j − P wsch,j ) P wr,j represents the rated output power of jth WT, and K Pw,j represents penalty cost coefficient of jth WT.
Similarly, ISO should pay a penalty cost for the additional amount of solar energy of kth SPV sources, and may be expressed below as: where K Ps,k represents the penalty cost coefficient of the kth SPV source, and f s (P sac,k > P ssch,k ) represents the possibility of the left-over amount SPV power. E(P sac,k > P ssch,k ) − P ssch,k represents the probability of P sac,k above P ssch,k . The cost function evaluation of SPV and wind turbines are based upon the distribution of the lognormal PDF and Weibull PDF, respectively.

Objective Function Formulation
For the OPF solution we formulate two main objective functions, one without carbon tax and the other with the carbon tax. Both objective functions incorporate all of the cost functions of all generators.
First objective function: where N s and N w represents the total number of SPV sources and WT in the network, respectively. Second objective function: While optimizing these OPF objective functions, system inequality, equality, and line capacity constraints must be satisfied.

OPF Equality Constraints
OPF has two equality constraints, that are real (MW) and reactive power (MVAr) balances in the network. The total real and reactive power produced must compensate for the total power loss and demand in the network.
Here, N B represents the total number of buses while B ij and G ij represent the susceptance and transfer conductance between bus i and j, respectively.

Inequality Constraints in OPF
For finding a solution to the OPF problem, there are several inequality constraints that should be satisfied. OPF inequality constraints are the generator voltages, real and reactive power limits, load bus voltages limits, and line capacity limits.
Generator inequality constraints are expressed as follow: Line and load bus inequality constraints:

Power Losses and Voltage Deviation
The most significant factor besides the reduction of total generation cost is the reduction of power loss in the network. Due to inherent line resistance, power losses are inevitable in transmission systems. The total network loss can be calculated as: The difference between the voltage angles of two buses is represented by δ ij = δ i − δ j ; G ij represents transfer conductance and nl shows the total number of lines in the network.
The voltage profile of the system has been expressed in terms of voltage deviation (VD), which is a gauge of the voltage quality within the system. Voltage deviation is computed by the following expression:

Mathematical Formulation for Stochastic SPV and WT Power
The adopted IEEE 30 bus system is changed at bus number 5 and 11 to WT, and at bus number 13 a SPV power source is integrated. Weibull PDF and wind speed distribution possess a similar shape [50]. A Weibull PDF with shape factor (k) and scale factor (c) can be used to obtain the probability of wind speed (m/s), as follows: Appl. Sci. 2021, 11, 6883 9 of 27 The output power of WT, which is a function of wind speed, can be modeled as follows: The values for the shape factor (k) and scale factor (c) are taken from reference [8] and are mentioned in Table 3. A Weibull fitting curve and wind speed frequency distribution curve are illustrated in Figures 2 and 3, which are obtained after simulation of 8000 Monte Carlo scenarios.
The values for the shape factor ( k ) and scale factor ( c ) are taken from reference [8] and are mentioned in Table 3. A Weibull fitting curve and wind speed frequency distribution curve are illustrated in Figures 2 and 3, which are obtained after simulation of 8000 Monte Carlo scenarios.
Output power from SPV sources relies on solar irradiance that follows the lognormal PDF [51]. The lognormal PDF, having a mean ( μ ) and standard deviation (σ ), may be used to obtain the probability for solar irradiance ( G ) by utilizing the following expression: The output power of the SPV source, which is a function of solar irradiance ( G ), can be mathematically modeled as follows: The parameters for the lognormal PDF are taken from reference [8] and are presented in Table 4. After simulation of 8000 Monte Carlo scenarios, solar irradiance frequency distribution and a lognormal fitting curve is obtained and is illustrated in Figure 4, while Figure 5 shows the real power distribution of the SPV source at bus 13.    Output power from SPV sources relies on solar irradiance that follows the lognormal PDF [51]. The lognormal PDF, having a mean (µ) and standard deviation (σ), may be used to obtain the probability for solar irradiance (G) by utilizing the following expression: The output power of the SPV source, which is a function of solar irradiance (G), can be mathematically modeled as follows: The parameters for the lognormal PDF are taken from reference [8] and are presented in Table 4. After simulation of 8000 Monte Carlo scenarios, solar irradiance frequency distribution and a lognormal fitting curve is obtained and is illustrated in Figure 4, while Figure 5 shows the real power distribution of the SPV source at bus 13.

Optimization Algorithm for OPF Solutions
Several optimization algorithms and techniques are currently being developed and have been applied for OPF solutions. The common objective of all of these nature inspired algorithms is trying to achieve the best global optimum with a fast convergence speed. In this paper a HPSO-GWO algorithm is utilized to provide the best OPF solution for the adopted test system.

Particle Swarm Optimization
The PSO is a probabilistic, population-based algorithm that is stimulated from the social behavior of bird flocking. While searching for food, they either fly together or in a scattered way before reaching the food. During the search for food, there must be a bird that can sense the location of the food precisely, where it can be found. The one having the correct information of the food position communicates the message to the group. During this search process, they change their positions continuously until they have found the food.
This behavior is inspired from the social behavior of birds, to find the best optimum value of the function/problem, where every bird of the flock is named a particle. The position of every particle in the PSO algorithm can be updated with the following equation:

Grey Wolf Optimization
Mirjalili et al. [26] in 2014 developed a GWO algorithm, which copies the leadership hierarchy and hunting mechanism of grey wolves. In the GWO algorithm, the hunting mechanism is implemented in three main steps to perform optimization, which are: searching, encircling, and attacking the prey. In the GWO algorithm, the leadership hierarchy is simulated by employing the grey wolves into four groups that are: alpha (α), beta (β), delta (α) and omega (γ) wolves. The alpha (α) wolves, which are a male and female, lead the pack and are liable for making the decisions of sleeping, hunting and many other aspects. The second upper level of wolves in the pack are known as beta (β) wolves, which are the secondary wolves that can help alphas (α) in decision making and other pack actions. The beta (β) wolves reinforce the alphas' (α) orders and provide feedback to alphas (α) during the pack actions. The third level of wolves in the pack are known as omegas (γ), which perform the character of a scapegoat. The omega (γ) wolves must always comply to all other leading alpha (α) and beta (β) wolves. The omegas (γ) are the last wolves that are permitted to eat, and do not have an important role in the hierarchy; however, if they lose the omegas (γ), the entire pack faces internal struggle and trouble. If a wolf does not belong to an alpha (α), beta (β) or omega (γ), she/he may be named as a delta (α) or subordinate. Delta (α) wolves dominate omegas, but they must comply to alpha (α) and beta (β) wolves in the command. To perform optimization, three main steps are to be executed, which are searching, encircling and hunting the prey.
The encircling movement of each wolf in the group is calculated through the following equation: where c and a are calculated as: Hunting behavior is mathematically expressed as: If A represents an arbitrary value between −2a and 2a, then an arbitrary value of |A| < 1 means that the wolves are enforced to encircle and attack the prey. While |A| > 1 means the wolves are forced to diverge from the prey. Exploration is the searching ability of GWO, whereas exploitation is the ability to attack the prey.

HPSO-GWO Optimizer
To find the global optimum solution with an excellent convergence capability is the common aim of all nature inspired metaheuristics algorithms. To find an optimal solution with better convergence performance, an algorithm should be furnished with the exploitation and exploration capabilities. Exploration is the capability of an algorithm to search the complete space of the solution, where exploitation is the capability to converge to the utmost result of the solution. The aim of all nature inspired hybrid algorithms is to establish a balance between these two capabilities to find a global optimal solution.
PSO has the best exploitation capability, but may become stuck at a local optimum, whereas the GWO has a best exploration capability of searching the whole search space, but has a slow convergence and poor local searching ability. A HPSO-GWO algorithm is developed in [52], in which the exploitation capability in PSO is improved with the exploration capability of the GWO to produce a combined algorithm strength.
In the developed HPSO-GWO algorithm, the position of delta, beta and alpha in the search space are updated with the following consecutive equations: In the above equations, both capabilities of grey wolves are balanced in the search area by an inertia constant. The flow chart of proposed HPSO-GWO algorithm is presented in Figure 6.
In the above equations, both capabilities of grey wolves are balanced in the search area by an inertia constant. The flow chart of proposed HPSO-GWO algorithm is presented in Figure 6.
In order to combine both variants (GWO and PSO) the equations for velocity and position are modified as follows:  In order to combine both variants (GWO and PSO) the equations for velocity and position are modified as follows:

Simulation Results and Discussion
The suggested HPSO-GWO algorithm is utilized for OPF solutions using MATLAB software. The simulation results are executed on a system with Intel (R) Core (TM) i5-4210U CPU @ 1.70 GHz 2.4 GHz with 4.0 GB RAM.
Simulations are carried out for several cases of adopted test systems, in which an SPV and two WT are integrated. Simulation results for each case obtained with utilization of the HPSO-GWO algorithm are presented with necessary description in this section. In the first two cases, variation in generation costs for SPV and WT are examined with variation in respective schedule powers and PDF parameters. In the third case, reduction of total generation cost is achieved without considering carbon tax and valve point effects for thermal generators. In the fourth case, total generation costs with valve point effects for thermal generators are considered. The fifth case optimizes the total generation cost with valve point effects of thermal generators, and a carbon tax of 20$/ton applied on emission. Each optimization algorithm is run five times, and a total of 500 iterations are executed in a whole run of the algorithm. During the complete run of each optimization algorithm for five times, the minimum cost of generation thus found and consequent control variable values are recorded.

Case 1: Total Cost, Reserve Cost, Direct Cost and Penalty Cost vs. PDF Parameters
In this case, total cost, penalty cost, direct cost and reserve cost of two WTs are plotted against Weibull PDF parameters, keeping the shape parameter (k = 2) constant while the scale parameter (c) is changed from 2 to 16 and respective changes in costs are plotted at a fixed scheduled WT power. A suitable schedule WT power of 25 MW and 20 MW are selected for WT 1 and WT 2, respectively. The schedule power selected is about one-third of the total installed capability, which is reasonable because practical WTs have a capacity between 30% and 45% [53]. The curves for WT 1 and WT 2 are plotted in Figures 7 and 8, respectively. The minimum value of total cost occurs around a value of 8 of the scale parameter (c); as the value of c increases, penalty cost increases while the reserve cost becomes constant after a certain value of c.
The suggested HPSO-GWO algorithm is utilized for OPF solutions using MATLA software. The simulation results are executed on a system with Intel (R) Core (TM) i 4210U CPU @ 1.70 GHz 2.4 GHz with 4.0 GB RAM.
Simulations are carried out for several cases of adopted test systems, in which an SP and two WT are integrated. Simulation results for each case obtained with utilization the HPSO-GWO algorithm are presented with necessary description in this section. In th first two cases, variation in generation costs for SPV and WT are examined with variatio in respective schedule powers and PDF parameters. In the third case, reduction of tot generation cost is achieved without considering carbon tax and valve point effects fo thermal generators. In the fourth case, total generation costs with valve point effects fo thermal generators are considered. The fifth case optimizes the total generation cost wi valve point effects of thermal generators, and a carbon tax of 20$/ton applied on emissio Each optimization algorithm is run five times, and a total of 500 iterations are executed a whole run of the algorithm. During the complete run of each optimization algorithm f five times, the minimum cost of generation thus found and consequent control variab values are recorded.

Case 1: Total Cost, Reserve Cost, Direct Cost and Penalty Cost vs. PDF Parameters
In this case, total cost, penalty cost, direct cost and reserve cost of two WTs are plotte against Weibull PDF parameters, keeping the shape parameter ( k = 2) constant while th scale parameter ( c ) is changed from 2 to 16 and respective changes in costs are plotted a fixed scheduled WT power. A suitable schedule WT power of 25 MW and 20 MW a selected for WT 1 and WT 2, respectively. The schedule power selected is about one-thir of the total installed capability, which is reasonable because practical WTs have a capaci between 30% and 45% [53]. The curves for WT 1 and WT 2 are plotted in Figures 7 and respectively. The minimum value of total cost occurs around a value of 8 of the sca parameter ( c ); as the value of c increases, penalty cost increases while the reserve co becomes constant after a certain value of c .
In a similar way, the SPV penalty cost, reserve cost and direct cost are plotted again lognormal PDF parameters from 2 to 7 at a fixed value of standard deviation (σ = 0.6). suitable value of scheduled SPV power is chosen at 20 MW. The minimum value of tot cost occurs at a lognormal mean value of 5.5; the penalty and reserve cost have the sam values at μ = 5.7. After this, the point penalty cost increases, due to which a shar increase occurs in total cost. SPV output power is highly dependent on lognormal mea μ i.e., at a low value of μ , a high reserve is necessary and vice versa. The variation SPV cost vs. lognormal mean ( μ ) is plotted in Figure 9.

Case 2: Total Cost, Penalty Cost, Reserve Cost and Direct Cost vs. Schedule Power
In this case, the scheduled wind power of two WT connected at bus 5 and 1 varied from zero to WT rated power; variation of total cost, direct cost, reserve cos penalty cost are plotted in Figures 8 and 9. The coefficients of direct cost, reserve cos penalty cost are taken from [8] and provided in Table 3. Simulation results show reserve cost decreases as scheduled WT power increases, while penalty costs increa scheduled power increases. There is a linear relationship between direct cost and schedule power. The total cost increased as schedule power increased, as it is the su direct cost, reserve cost and penalty cost.
Similarly, SPV total cost, direct cost, reserve cost and penalty costs are plotted ag SPV scheduled power in Figure 10. The SPV different parameters and cost coefficien provided in Table 4. The total minimum cost occurs somewhere between 15 MW an MW in the case of solar power, while in the case of WT sources there is a linear rel between cost and wind schedule power, which can be observed from Figures 10-12. In a similar way, the SPV penalty cost, reserve cost and direct cost are plotted against lognormal PDF parameters from 2 to 7 at a fixed value of standard deviation (σ = 0.6). A suitable value of scheduled SPV power is chosen at 20 MW. The minimum value of total cost occurs at a lognormal mean value of 5.5; the penalty and reserve cost have the same values at µ = 5.7. After this, the point penalty cost increases, due to which a sharp increase occurs in total cost. SPV output power is highly dependent on lognormal mean µ i.e., at a low value of µ, a high reserve is necessary and vice versa. The variation in SPV cost vs. lognormal mean (µ) is plotted in Figure 9.

Case 2: Total Cost, Penalty Cost, Reserve Cost and Direct Cost vs. Schedule Power
In this case, the scheduled wind power of two WT connected at bus 5 and 1 varied from zero to WT rated power; variation of total cost, direct cost, reserve cos penalty cost are plotted in Figures 8 and 9. The coefficients of direct cost, reserve cos penalty cost are taken from [8] and provided in Table 3. Simulation results show reserve cost decreases as scheduled WT power increases, while penalty costs increa scheduled power increases. There is a linear relationship between direct cost and w schedule power. The total cost increased as schedule power increased, as it is the su

Case 2: Total Cost, Penalty Cost, Reserve Cost and Direct Cost vs. Schedule Power
In this case, the scheduled wind power of two WT connected at bus 5 and 11 are varied from zero to WT rated power; variation of total cost, direct cost, reserve cost and penalty cost are plotted in Figures 8 and 9. The coefficients of direct cost, reserve cost and penalty cost are taken from [8] and provided in Table 3. Simulation results show that reserve cost decreases as scheduled WT power increases, while penalty costs increase as scheduled power increases. There is a linear relationship between direct cost and wind schedule power. The total cost increased as schedule power increased, as it is the sum of direct cost, reserve cost and penalty cost.
Similarly, SPV total cost, direct cost, reserve cost and penalty costs are plotted against SPV scheduled power in Figure 10. The SPV different parameters and cost coefficients are provided in Table 4. The total minimum cost occurs somewhere between 15 MW and 20 MW in the case of solar power, while in the case of WT sources there is a linear relation between cost and wind schedule power, which can be observed from Figures 10-12.

Case 3: Minimization of Total Generation Cost without Valve Point Effect
In this case, total generation costs without considering valve point effects of thermal generators is optimized. Cost coefficients and all other parameters for wind and SPV sources remain the same, as in Case 2. The optimized values obtained for all thermal generators, RES and other computed parameters are given in Table 5. The obtained results from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 5 for comparison. The convergence curve of four selected algorithms is shown in Figure 13. The voltage deviation VD, power losses, emissions and total generation costs are also provided in the table, and the best value obtained is highlighted in bold. The minimum generation cost obtained by HPSO-GWO is 772.736 $/h, and the worst result obtained by SHADE-SF is 773.979 $/h. The difference between the total generation cost obtained by HPSO-GWO and SHADE-SF is 1.243 $/h, which becomes a significant amount of (1.243 × 24 × 365) 10,889 $/Y. The minimum power losses and emissions are achieved by HPSO-GWO, while a minimum voltage deviation is achieved by SHADE-SF. Generator bus voltages, real powers and reactive powers for all algorithms are combined in Figures 14-16, respectively. All of the values obtained by these algorithms are within specified limits, as shown clearly. The load bus voltages profile obtained through all algorithms is combined in Figure 17. The five run results of each algorithm are combined in Figure 18, and the minimum value achieved by each algorithm in the respective run of the simulation is indicated by bold.  generators is optimized. Cost coefficients and all other parameters for wind and sources remain the same, as in Case 2. The optimized values obtained for all th generators, RES and other computed parameters are given in Table 5. The obtained r from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 5 for compa The convergence curve of four selected algorithms is shown in Figure 13. The vo deviation VD, power losses, emissions and total generation costs are also provided  Figure 17. The five run results of each algorithm are combined in Figure 18, an minimum value achieved by each algorithm in the respective run of the simulat indicated by bold.  In this case, total generation costs without considering valve point effects of th generators is optimized. Cost coefficients and all other parameters for wind and sources remain the same, as in Case 2. The optimized values obtained for all th generators, RES and other computed parameters are given in Table 5. The obtained r from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 5 for compa The convergence curve of four selected algorithms is shown in Figure 13. The vo deviation VD, power losses, emissions and total generation costs are also provided table, and the best value obtained is highlighted in bold. The minimum generation obtained by HPSO-GWO is 772.736 $/h, and the worst result obtained by SHADE 773.979 $/h. The difference between the total generation cost obtained by HPSO-GWO SHADE-SF is 1.243 $/h, which becomes a significant amount of (1.243 × 24 × 365) 1 $/Y. The minimum power losses and emissions are achieved by HPSO-GWO, w minimum voltage deviation is achieved by SHADE-SF. Generator bus voltages powers and reactive powers for all algorithms are combined in Figures 1  respectively. All of the values obtained by these algorithms are within specified lim shown clearly. The load bus voltages profile obtained through all algorithms is com in Figure 17. The five run results of each algorithm are combined in Figure 18, an minimum value achieved by each algorithm in the respective run of the simulat indicated by bold.

Case 4: Minimization of Total Generation Cost with Valve Point Effect
In this case, total generation cost with valve point effect of thermal generators is considered as the objection function for optimization. Cost coefficients and all other parameters for wind and SPV sources remain the same as in Case 2. The obtained results from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 6 for comparison with the results obtained by BMO and GOA [46]. Comparing the results of Table 6 with  Table 5, it is clear that penetration from thermal generators decreases as we add up valve point effects in the cost function. The convergence curve of four selected algorithms is shown in Figure 19. The voltage deviation VD, power losses, emissions and total generation costs are also provided in the table, and the best value obtained is highlighted

Case 4: Minimization of Total Generation Cost with Valve Point Effect
In this case, total generation cost with valve point effect of thermal generators is considered as the objection function for optimization. Cost coefficients and all other parameters for wind and SPV sources remain the same as in Case 2. The obtained results from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 6 for comparison with the results obtained by BMO and GOA [46]. Comparing the results of Table 6 with Table 5, it is clear that penetration from thermal generators decreases as we add up valve point effects in the cost function. The convergence curve of four selected algorithms is shown in Figure 19. The voltage deviation VD, power losses, emissions and total generation costs are also provided in the table, and the best value obtained is highlighted in bold. The minimum generation cost obtained by HPSO-GWO is 780.587 $/h, and the worst result obtained by GOA is 785.711 $/h. The difference between total generation cost obtained by HPSO-GWO and GOA is 5.124 $/h, which become a significant amount of (5.124 × 24 × 365) 44,886.24 $/Y. The minimum power losses and emissions are achieved by PSO, while a minimum voltage deviation is achieved by SHADE-SF. Generator bus voltages, real powers and reactive powers for all algorithms are combined in Figures 20-22, respectively. All of the values obtained by these algorithms are within the defined limits, as shown clearly. The load bus voltages obtained by all algorithms are combined in Figure 23. The five run results of each algorithm are combined in Figure 24, and the minimum value achieved by each algorithm in the respective run of each simulation is indicated by bold.   Figure 23. The five run results of each algorithm are combined in Figure 24, and the minimum value achieved by each algorithm in the respective run of each simulation is indicated by bold.

Case 5: Minimization of Total Generation Cost with Valve Point Effect and Carbon Tax
In this case, total generation cost with valve point effects and carbon tax imposed on emission of thermal generators is considered as the main objective function for optimization. The rate of carbon tax is kept the same, as in [8] and [46], that is 20$/ton. Cost coefficients and all other parameters for wind and SPV sources remained the same, as in Case 2. The obtained results from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 7 for a comparison with the results achieved by BMO and GOA [46]. Due to carbon tax imposed on emissions of thermal generators, their penetration decreases whereas RES penetration increases in the system; this can be noticed by

Case 5: Minimization of Total Generation Cost with Valve Point Effect and Carbon Tax
In this case, total generation cost with valve point effects and carbon tax imposed on emission of thermal generators is considered as the main objective function for optimization. The rate of carbon tax is kept the same, as in [8] and [46], that is 20$/ton. Cost coefficients and all other parameters for wind and SPV sources remained the same, as in Case 2. The obtained results from PSO, GWO, HPSO-GWO and SHADE-SF are summarized in Table 7 for a comparison with the results achieved by BMO and GOA [46]. Due to carbon tax imposed on emissions of thermal generators, their penetration decreases whereas RES penetration increases in the system; this can be noticed by comparing the results of Tables 5-7. The convergence curve of four selected algorithms is shown in Figure 25. The voltage deviation VD, power losses, emissions and total generation costs are also provided in the table, and the best value obtained is highlighted in bold. The minimum generation cost obtained by HPSO-GWO is 809.277 $/h, and the worst result obtained by GOA is 822.307 $/h. The difference between total generation cost obtained by HPSO-GWO and GOA is 13.03 $/h, which become a significant amount of (13.03 × 24 × 365) 114,142.8 $/Y. The minimum power losses and emissions are achieved by GWO while a minimum voltage deviation is achieved by SHADE-SF. Generator bus voltages, real powers and reactive powers for all algorithms are combined in Figures 26-28, respectively. All of the values obtained by these algorithms are within specified limits, as shown clearly. The load bus voltages obtained by all algorithms are illustrated in Figure 29. The five run results of each algorithm are combined in Figure 30, and the minimum value achieved by each algorithm in the respective run of each simulation is indicated by bold.  comparing the results of Tables 5-7. The convergence curve of four selected algorithms is shown in Figure 25. The voltage deviation VD, power losses, emissions and total generation costs are also provided in the  Figure  29. The five run results of each algorithm are combined in Figure 30, and the minimum value achieved by each algorithm in the respective run of each simulation is indicated by bold.    Figure  29. The five run results of each algorithm are combined in Figure 30, and the minimum value achieved by each algorithm in the respective run of each simulation is indicated by bold.

Conclusions
This paper presents a solution to OPF problems, integrated with RES using a HPSO-GWO algorithm. The stochastic nature of SPV and WT is modeled with Lognormal and Weibull probability density functions, respectively. Important system constraints such as equality constraints, inequality constraints and line capacity constraints have been taken into consideration. The usefulness and robustness of the proposed HPSO-GWO algorithm is tested on adopted test systems integrated with RES. Based upon the simulation results obtained, it has been concluded that:

Conclusions
This paper presents a solution to OPF problems, integrated with RES using a HPSO-GWO algorithm. The stochastic nature of SPV and WT is modeled with Lognormal and Weibull probability density functions, respectively. Important system constraints such as equality constraints, inequality constraints and line capacity constraints have been taken into consideration. The usefulness and robustness of the proposed HPSO-GWO algorithm is tested on adopted test systems integrated with RES. Based upon the simulation results obtained, it has been concluded that:

•
The proposed HPSO-GWO algorithm has been effectively applied to adopted test systems to find the optimal values of the control settings. • Minimization of total system generation costs without considering valve point effect, with valve effect, and valve effect with carbon tax imposed on emission) is an objective function which has been achieved using the HPSO-GWO algorithm.

•
The dominance of the HPSO-GWO algorithm in comparison with the PSO, GWO, SHADE-SF, GOA and BMO algorithms has been confirmed.

•
The HPSO-GWO has a speedy and steady convergence characteristic curve in comparison with PSO, GWO and SHADE-SF. • Penetration from renewable energy sources further increased, and emissions have been reduced by 0.864 ton/h after imposing a carbon tax on emissions.
In future, the HPSO-GWO algorithm will be effectively utilized in other power system expansions and planning studies, with the possible addition of hydro power generation in the system.