Solving the Optimal Reactive Power Dispatch Using Marine Predators Algorithm Considering the Uncertainties in Load and Wind-Solar Generation Systems

: The optimal reactive power dispatch (ORPD) problem is an important issue to assign the most e ﬃ cient and secure operating point of the electrical system. The ORPD became a strenuous task, especially with the high penetration of renewable energy resources due to the intermittent and stochastic nature of wind speed and solar irradiance. In this paper, the ORPD is solved using a new natural inspired algorithm called the marine predators’ algorithm (MPA) considering the uncertainties of the load demand and the output powers of wind and solar generation systems. The scenario-based method is applied to handle the uncertainties of the system by generating deterministic scenarios from the probability density functions of the system parameters. The proposed algorithm is applied to solve the ORPD of the IEEE-30 bus system to minimize the power loss and the system voltage devotions. The result veriﬁes that the proposed method is an e ﬃ cient method for solving the ORPD compared with the state-of-the-art techniques.


Introduction
Solving the optimal reactive power dispatch (ORPD) problem plays a significant role in the efficient operation and planning of the power system. The aim of solving the ORPD is to determine the best operating point of system for maximizing the voltage stability, minimizing the system loss and the voltage deviations [1,2]. The best operating point includes the terminal voltages of the generators, taps of transformers and the injected reactive powers of the shunt compensators [3].
The ORPD is a nonconvex and nonlinear optimization problem. Several conventional techniques such as the quadratic programming method [4], the linear programming [5], the nonlinear programming [6], and the interior point [7] have solved the ORPD problem. However, although these techniques succeed in solving the problem for some cases, they suffer from stagnation at the local optimal for other cases, especially in solving optimal active and reactive power flows of the large-scale systems [8]. The shortages of the aforementioned techniques can be avoided by applying numerous optimization techniques, which can be classified based on their inspiration as follows: • Swarm-based algorithms such as particle swarm optimization (PSO) [9], ant lion optimizer [10], whale optimization algorithm [11], improved social spider optimization algorithm [12], improved antlion optimization algorithm [13], and moth swarm algorithm [14]; • Evolutionary-based algorithms such as differential evolution [15], specialized genetic algorithm (SGA) [16], evolutionary programming [17], modified differential evolution [18], pareto evolutionary algorithm [19], comprehensive learning particle swarm optimization [20], and enhanced grey wolf optimizer (EGWO) [21]; • Human-based algorithms such as harmony search algorithm [22], teaching learning-based optimization [23], and biogeography-based optimization [24]; • Physical-based algorithms such as gravitational search algorithm [25], improved gravitational search algorithm [26], lightning attachment procedure optimization (LAPO) [27], modified sine cosine algorithm [28], and water cycle algorithm [29]; • Hybrid-based algorithms such as hybridization of the particle swarm optimization method and the tabu-search technique [30].
Renewable energy resources (RERs) widely incorporated in electrical systems provide an elegant solution from an economic, environmental and technical perspective. In other words, the inclusion of the RERs can effectively minimize the dependence of the generation obtained by fossil fuel, reduce the greenhouse gases and the harmful emissions, minimize the generation cost and enhance the system operation. The wind and solar power generation systems are the most applied technologies for RERs. However, some technical issues related to incorporating the wind and solar generation units, such as the wind speed and the solar irradiance, are continuously varied and intermittent in nature, which leads to increasing the uncertainties of the electrical power system, especially fluctuations in load demand [31,32]. Considering the uncertainties of the electrical power system is a strenuous task for the decision makers for efficient planning [33]. Thus, several approaches have been presented for modeling the uncertainty of the system including probabilistic methods [34], possibilistic methods [35], hybrid possibilistic-probabilistic methods [36], information gap theory [37], robust optimization [38], and interval analysis [39]. Moreover, the authors in [40] presented a comprehensive review of the stochastic techniques that have been implemented for optimization of solar-based renewable energy systems. In [41], a distributed operation strategy using a double deep Q-learning method is employed to manage the operation of the battery energy storage system considering the uncertainty of the renewable distributed generators. The authors in [42] proposed a robust optimization model for analyzing the interdependency of natural gas, coal and electricity infrastructures considering their operation constraints and wind power uncertainties.
A marine predators algorithm (MPA) is one of the newest optimization techniques that simulates the foraging behavior and movement of the marine predators proposed by A. Faramarzi et al. in 2020 [43]. The authors in [44] have applied the MPA to determine the optimal parameters for the adaptive neuro-fuzzy inference system to predict the spread of coronavirus . In [45], an improved version of the MPA was presented to assign the people infected by COVID-19 with X-ray image segmentation.
The aim of our study is to apply the newest algorithm, called the marine predators algorithm (MPA), in order to solve the ORPD problem for the first time to best of our knowledge. Then, the validity of the proposed MPA for minimizing the power loss and the total voltage deviations is investigated and compared with state-of-the-art techniques. The contributions of this paper can be summarized as follows: 1.
Solving the ORPD problem by utilizing one of the newest algorithms, called the marine predators algorithm (MPA); 2.
The validity of the proposed algorithm for minimizing the power loss and the total voltage deviations is investigated and compared with the state-of-the-art techniques; 3.
The ORPD problem is solved by incorporating the renewable energy resources, including wind turbine and solar PV systems; 4.
The ORPD problem is solved to minimize the expected power loss with considering the uncertainties of the load demand and the output powers of wind and solar generation systems; 5.
Weibull PDF, Beta PDF and normal PDF are utilized for modeling the uncertainties of wind speed, the solar irradiance and the load demand, respectively. In addition, the scenario-based method is used to model the combination of load-generation uncertainty.
The organization of paper is listed as follows: Section 2, illustrates the problem formulations of the ORPD problem. Section 3 models the uncertainty of the renewable sources and the load demand. Section 4 illustrates the step procedure of the MPA. Section 5 shows the obtained results and the corresponding discussion. Finally, the paper's conclusion is presented in Section 6.

Problem Formulation
The solution of the ORPD problem is formulated as an optimization problem applied for assigning set control parameters for a certain objective function, satisfying the operating constraints of the system. Generally, the ORPD problem is represented as follows where g k and h n represent the equality and inequality constraints. u is a vector of the control parameters which includes the generator voltages, the compensators of injected reactive powers and the transformer taps while x denotes the vector of the dependent variables which includes the slack bus power, the voltages of the load buses and the apparent power flow in transmission lines. u and x are vectors represented as follows

Objective Function
Power Loss Voltage Deviations

Constraints
Inequality Constraints P min Gk ≤ P Gk ≤ P max Gk k = 1, 2, . . . , N G (8) T min n ≤ T n ≤ T max n n = 1, 2, . . . , N Q Q min Cn ≤ Q Cn ≤ Q max Cn n = 1, 2, . . . , N Q Quality Constraints The system constraints should be considered to ensure that the obtained solution is at the suitable solution. This can be accomplished by using the concept of the weighted square variables as follows

Uncertainty Modeling
To consider the uncertainties of the load demand and output powers of wind and solar photovoltaic (PV) generation systems, the continuous probability density function (PDF) is used for modeling the uncertainties of the system where the PDF is divided into subsections to obtain a number of scenarios from the load demand, wind speed and solar irradiance.

Modeling of Load Demand
Uncertainty of load demand is modeled using normal PDF [46] which can be described as where P d is the probability density of normal distribution of the load while σ d and µ d are the standard and mean deviation values, respectively. The portability of load demand and its corresponding expected load scenario can be calculated as follows [47] where P max d,i and P min d,i are the minimum and the maximum limit of the selected interval i. In this paper, three scenarios of load demand are generated from the previous equations as depicted in Figure 1 and their corresponding mean values are µ d − 1.5259σ d and µ d + 1.5259σ d , while their probabilities are 0.1587, 0.6826 and 0.1587, respectively. In this paper, the selected value of the σ d = 0.02µ d . The load scenarios and their corresponding probability are depicted in Table 1.

Modeling of Wind Speed
The uncertainty of wind speed is modeled by using the Weibull probability density function ( ) which can be represented as follows [48] where and are the scale and the shape parameters of the Weibull PDF. The output power of wind turbine is defined as follows [49] where is the rated output power of the wind turbine while , and are the rated, cutin and the cut-out speeds of the wind turbine. In this paper, the wind farm consists of 25 turbines and the rated power of each turbine is 3 MW, while its rated speeds of , and are 16, 3 and 25 m/s, respectively [50]. The portability of wind speed for each scenario can be calculated as follows [32] , = ( )

Modeling of Wind Speed
The uncertainty of wind speed is modeled by using the Weibull probability density function f v (v) which can be represented as follows [48] f where α and β are the scale and the shape parameters of the Weibull PDF. The output power of wind turbine is defined as follows [49] where P wr is the rated output power of the wind turbine while v ωr , v ωi and v ωo are the rated, cut-in and the cut-out speeds of the wind turbine. In this paper, the wind farm consists of 25 turbines and the rated power of each turbine is 3 MW, while its rated speeds of v ωr , v ωi and v ωo are 16, 3 and 25 m/s, respectively [50]. The portability of wind speed for each scenario can be calculated as follows [32]  where τ wind,k denotes the probability of the wind speed being in scenario k. In addition, v min k and v max k are the starting and ending points of wind speed's interval at k-th scenario. In this paper, three scenarios are generated from the previous equations. The wind speed scenarios and their corresponding probability are depicted in Table 1.

Modeling of Solar Irradiance
The uncertainty of solar irradiance called (G) is modeled by the Beta PDF, which is used to describe the solar irradiance [51]. Thus, the Beta PDF is formulated as follows where α and β are parameters of the beta probability function which can be calculated in terms of the mean (µ s ) and standard deviation (σ s ) of the random variables as follows The output power of the PV unit is calculated as a function of solar irradiance as follows [52,53] where P sr is the rated power of the solar PV. G std denotes the standard solar irradiance which is set as 1000 W/m 2 . Moreover, X c represents a certain irradiance point which is set as 120 W/m 2 [50]. In this paper, the rated power of the PV system is 50 MW. The portability of solar irradiance for each scenario can be calculated as follows [32] τ Solar,m = G min m G min m f G (G)dG (28) where τ Solar,i denotes the probability of the solar irradiance being in scenario m. G min m and G min m are the starting and ending points of solar irradiance's interval at m-th scenario. In this paper, three scenarios are generated from the previous equations. The solar irradiance scenarios and their corresponding probability are depicted in Table 1. It should be highlighted here that the three most occurring scenarios and the corresponding probabilities (τ d,i , τ wind,k , τ Solar,m ) are selected similarly to [32] and [47].

The Combined Load-Generation Model
Combined scenarios of the load, wind speed and solar irradiance model are captured by multiplying the probabilities of (19), (23) and (28) together which result as

Optimization Algorithm
The marine predators algorithm (MPA) is a new efficient algorithm that is conceptualized from the foraging behavior of the marine predators like marlines, tunas, sunfish, sharks, and swordfish with their prey in oceans. The foraging technique of the marine predators depends upon two random movements processes including the Lévy flight walk and the Brownian movements which are depicted in Figure  predators when searching for food; however, when it comes to foraging, the pattern is prevalently switched to Brownian type [54]. It should be highlighted here that the Lévy flight walk is a random process of transition of an object from one position to another based on the probability distribution factor [43].
Energies 2020, 13, x FOR PEER REVIEW 7 of 19 depicted in Figure 2. Humphries et al. indicated that Lévy motion is a widespread pattern among marine predators when searching for food; however, when it comes to foraging, the pattern is prevalently switched to Brownian type [54]. It should be highlighted here that the Lévy flight walk is a random process of transition of an object from one position to another based on the probability distribution factor [43]. The predator moves based on the Lévy flight at the food in a prey-sparse environment. On the other hands, the predator moves in a Brownian pattern when being located in a prey-abundant area. Another behavior related to the motion of the predators like sharks is recorded by the fish aggregating device (FAD), which reveals that the sharks move in a sudden vertical jump. The following steps describe the MPA: Step 1: initialization The initialization of the populations of the MPA randomly as follows where and are the maximum and the minimum boundaries of the control variable; rand is a random number where 0 ≤ ≤ 1. The objective functions of the initial population are represented as follows Step 2: Detecting top predator First, the populations are arranged in a matrix called the prey matrix which can be represented as follows where is the number of population, while denotes the number of control variables. In this step, the top predator is determined by the arrangement of the solutions based on their objective function values. A matrix is constructed that includes the top predator called the Elite matrix, which can be listed as follows where denotes the top predator vector.
Step 3: Lévy flight and the Brownian movements The predator moves based on the Lévy flight at the food in a prey-sparse environment. On the other hands, the predator moves in a Brownian pattern when being located in a prey-abundant area. Another behavior related to the motion of the predators like sharks is recorded by the fish aggregating device (FAD), which reveals that the sharks move in a sudden vertical jump. The following steps describe the MPA: Step 1: initialization The initialization of the populations of the MPA randomly as follows where X max i and X min i are the maximum and the minimum boundaries of the control variable; rand is a random number where 0 ≤ rand ≤ 1. The objective functions of the initial population are represented as follows Step 2: Detecting top predator First, the populations are arranged in a matrix called the prey matrix which can be represented as follows where n is the number of population, while d denotes the number of control variables. In this step, the top predator is determined by the arrangement of the solutions based on their objective function values. A matrix is constructed that includes the top predator called the Elite matrix, which can be listed as follows where E denotes the top predator vector.
Step 3: Lévy flight and the Brownian movements In this step, the positions of the prey and the predators are updated based on three phases which depend upon the velocity ratio between the prey and predator velocity, which can be stated as follows: Phase 1: This phase, representing the exploration phase of the algorithm, is applied at a high-velocity ratio. In other words, the velocity of the predator is higher than the velocity of the prey. In this phase, the prey and the predator are updated based on Brownian movement which can be mathematically represented as follows where T and T max are the current iteration number and the maximum number of iterations, respectively; ⊕ denotes the entry wise multiplication; R Br is a vector containing random numbers based on normal distribution representing the Brownian motion. In addition, SZ i is the step size vector, and P is a constant number equal to 0.5. Phase 2: This phase is an intermittent phase between the exploration and the exploitation phase that is applied when the velocity of predator equals the velocity of prey. In this section, the populations are divided into two groups. The first group is employed from exploitation and the other group for exploration, which can be represented as follows The second group or exploration group where R Levy is a vector containing random numbers based on Lévy distribution which mimics the movement of prey in levy manner. CF is an adaptive operator utilized to control the step size of predator's movement. Phase 3: This phase is a fully exploitation phase which is applied when the velocity of the predator is more than the velocity of the prey at the final iterations of the optimization technique.
In this phase, the predator is moving in levy strategy and the previous equations describe the prey movement based on the elite vector.
Step 4: In eddy formation and FADs' effect, the predator changes its movement behavior due to environmental issues, as mentioned before. These predators move in the eddy formation or fish aggregating devices which represent the local optima, while they take a longer jump to find a new environment that has abundant regions. This step can be represented as follows where r is a random value within the range 0 to 1. r1 and r2 represent random indices from prey matrix. FADS represents FADs probability, which equals 0.2. U is a binary vector.
Step 5: Marine memory The marine predators can efficiently remember the best location of foraging. Thus, in the MPA technique, the updated solution is compared with those in the previous iteration to capture the optimal solution.

Simulation Results
In this section, the MPA technique is utilized to solve the ORPD in IEEE 30-bus system with and without considering the uncertainty of system. The program code of the proposed algorithm of ORPD was written for the MATLAB R2018b programming software and run on a PC with Core I5 @ 1.7 GHz with 4GB RAM. The IEEE 30-bus system data are given in [55]. The boundaries of the control variables are listed in Table 2. The studied cases of system are listed as follows. Table 2. The optimal control variables for P Loss and voltage deviations(VD).

Case 2: Solution of ORPD Problem Considering the Uncertainty
In this case, the ORPD is solved by considering the uncertainties of the wind and solar generation systems and load demand. To solve the ORPD with the stochastic natural of wind and solar generation resources, the IEEE 30-bus system is modified as depicted in Figure 5. A wind farm and solar PV system are incorporated with bus 5 and bus 8, respectively. The wind farm consists of 25 turbines and the rated power of each turbine is 3 MW, while the rated speeds of each turbine v ωr , v ωi and v ωo are 16, 3 and 25 m/s, respectively. The rated power of the PV system is 50 MW and the standard solar irradiance (G std ) is 1000 W/m 2 [50].

Case 2: Solution of ORPD Problem Considering the Uncertainty
In this case, the ORPD is solved by considering the uncertainties of the wind and solar generation systems and load demand. To solve the ORPD with the stochastic natural of wind and solar generation resources, the IEEE 30-bus system is modified as depicted in Figure 5. A wind farm and solar PV system are incorporated with bus 5 and bus 8, respectively. The wind farm consists of 25 turbines and the rated power of each turbine is 3 MW, while the rated speeds of each turbine , and are 16, 3 and 25 m/s, respectively. The rated power of the PV system is 50 MW and the standard solar irradiance ( ) is 1000 W/m 2 [50]. As depicted in Table 1, this shows three individual scenarios and their corresponding probabilities for modeling the uncertainties of the load, the solar irradiance and the wind speed. The number of the generated scenarios from combining the uncertainties of load, solar irradiance and the wind speed is 27 according to Section 3.4. Table 5 lists the combined scenarios and their probabilities  As depicted in Table 1, this shows three individual scenarios and their corresponding probabilities for modeling the uncertainties of the load, the solar irradiance and the wind speed. The number of the generated scenarios from combining the uncertainties of load, solar irradiance and the wind speed is 27 according to Section 3.4. Table 5 lists the combined scenarios and their probabilities along with the load, solar irradiance and the wind speed. It should be highlighted here that the wind speeds in Table 5 are a percentage of the rated wind speed (16 m/s) while the solar irradiances are a percentage of the standard solar irradiance (1000 W/m 2 ). The aim of solving the ORPD is to minimize the expected the power loss over the whole scenarios which can be expressed as follows where N S is a number of the generated scenarios; EPL h is the expected power loss for the h-th scenario; TEPL describes the total expected power losses. The total expected power loss without incorporating the renewable energy resources in system is 10.746 MW, while it equals 7.1223 MW when incorporating the renewable energy resources and considering their probabilities. In other words, the TEPL is reduced to 33.72% with the inclusion of the renewable energy resources. Table 6 shows the output power of wind turbine, solar PV unit and the corresponding probabilities as well as the expected power loss for each scenario. The IEEE 30-bus voltage profile for each scenario is depicted in Figure 6. It is clear that the system voltages are within limits.

Conclusions
This paper solved the ORPD problem efficiently using a new efficient algorithm called the marine predators algorithm (MPA) for losses and the voltage deviations minimization. The ORPD was solved with and without the uncertainties of the system. The proposed technique was validated in an IEEE 30-bus system and the obtained results were compared with well-known techniques. Three uncertainty parameters were considered, which are the uncertainty of load demands and uncertainties associated with renewable energy resources including the wind speed and solar irradiance. The uncertainties of load demands, wind speed and the solar irradiance were examined using normal PDF, Weibull PDF and Beta PDF. The scenario-based model has been applied to generate various scenarios using the individual scenarios of the uncertainty parameters of the system. The expected power loss is evaluated with numerous scenarios of load demands, wind and solar power. The obtained results confirmed that the proposed MPA is an efficient technique for solving the ORPD compared with other reported techniques. Furthermore, in cases solving the ORPD considering the uncertainty, the expected power loss reduced considerably with the inclusion of renewable energy resources, from 10.746 to 7.1223 MW.
Author Contributions: M.E. and A.A. proposed the idea, wrote the simulation program, obtained results and wrote the important parts of the article. S.K. and F.J. contributed by drafting and critical revisions and summarized results in tables. F.J. has edited the whole paper. All authors together organized and refined the

Conclusions
This paper solved the ORPD problem efficiently using a new efficient algorithm called the marine predators algorithm (MPA) for losses and the voltage deviations minimization. The ORPD was solved with and without the uncertainties of the system. The proposed technique was validated in an IEEE 30-bus system and the obtained results were compared with well-known techniques. Three uncertainty parameters were considered, which are the uncertainty of load demands and uncertainties associated with renewable energy resources including the wind speed and solar irradiance. The uncertainties of load demands, wind speed and the solar irradiance were examined using normal PDF, Weibull PDF and Beta PDF. The scenario-based model has been applied to generate various scenarios using the individual scenarios of the uncertainty parameters of the system. The expected power loss is evaluated with numerous scenarios of load demands, wind and solar power. The obtained results confirmed that the proposed MPA is an efficient technique for solving the ORPD compared with other reported techniques. Furthermore, in cases solving the ORPD considering the uncertainty, the expected power loss reduced considerably with the inclusion of renewable energy resources, from 10.746 to 7.1223 MW.
Author Contributions: M.E. and A.A. proposed the idea, wrote the simulation program, obtained results and wrote the important parts of the article. S.K. and F.J. contributed by drafting and critical revisions and summarized results in tables. F.J. has edited the whole paper. All authors together organized and refined the manuscript in the present form. All authors have read and agreed to the published version of the manuscript.
Funding: This research receives no external funding.

Acknowledgments:
In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

Conflicts of Interest:
The authors declare no conflict of interest. Penalty factors max, min : Superscript of maximum and minimum limit lim

Abbreviations
Superscript of limit boundary f d Probability density function of load f v Weibull probability density function f G Beta probability density function of solar irradiance σ d , µ d The standard and mean deviation values of the load demand EPL The expected power loss TEPL The total expected power losses