An Improved Marine Predators Algorithm for Optimal Reactive Power Dispatch With Load and Wind-Solar Power Uncertainties

With the growing usage of renewable energy resources, the fluctuation and randomness of solar and wind power output have a significant influence on the safe and stable operation of a power system. Deciding how to deal with these uncertainties has become a hot topic for researchers and practitioners. In this paper, a hybrid improved marine predators algorithm is proposed to solve the optimal reactive power dispatch problem with load demand, and wind-solar power uncertainties. The $\varepsilon $ -constraint handling technique is adopted to deal with the constraints in the optimization problem. The uncertainty factors in the system are modeled using appropriate probability density functions. The backward reduction method is used to determine a representative set of scenarios from a large number of scenarios produced by a Monte Carlo simulation method. IEEE 30-bus test system is utilized to verify the proposed algorithm. The simulation results demonstrate that the algorithm can produce better results than other advanced metaheuristic algorithms.


I. INTRODUCTION
The optimal reactive power dispatch (ORPD) problem is one of the significant research directions in the field of power system optimization. The primary objective of the ORPD problem is to maximize the voltage stability by optimizing the generator terminal voltages, transformer taps, shunt VAR compensators, and other equipment to minimize active power loss and voltage deviation [1], [2]. From the mathematical viewpoint, the control variables of the power system include both continuous variables (i.e., generator voltages) and discrete variables (i.e., tap of transformers and shunt VAR compensators), which make the ORPD problem difficult to The associate editor coordinating the review of this manuscript and approving it for publication was Sotirios Goudos . solve. Therefore, seeking a better solution has received much attention.
In the early stages, many traditional mathematical methods, including the newton method [3] and linear programming [4], have been widely employed to address the ORPD problem. Unfortunately, these techniques have weaknesses that cannot be ignored, in particular poor convergence, complex computation, and inappropriate handling of discrete variables. To better solve the non-linear ORPD problem, scholars have developed numerous intelligent optimization techniques, like differential evolution (DE) [5], grey wolf optimizer (GWO) [6], cuckoo search algorithm (CS) [7], bat algorithm (BA) [8], and social spider optimization (SSO) [9].
Nowadays, with the increasing use of renewable energy resources (RERs) in the power system, the uncertainties of wind power and photovoltaic power output has brought serious challenges to safe and stable operation. To realize efficient management, the uncertainty of the system must be taken into account. Therefore, the solutions of ORPD incorporating RERs are studied in [10], [11], [12], and [13]. The marine predators algorithm (MPA) is a recently proposed algorithm by Faramarzi et al. [14]. The research demonstrates that while the MPA performs favorably with different algorithms, it still has poor convergence and local optima traps when solving complex problems [15], [16], [17].
In this paper, considering the uncertainties of load demand and RERs, a hybrid improved marine predators algorithm (HIMPA) is applied to solve the ORPD problem. The HIMPA has been modified in the following three aspects: i) to improve the quality of the initial solution, an oppositionbased learning method is adopted; ii) to balance the exploration and development capabilities, the mid-term stage of traditional MPA is replaced by the modified DE/best − to-rand/1 strategy; iii) SS operator is used to enhance the development stage of the algorithm and accelerate the convergence speed. Since the constraints are difficult to deal with, we use the ε-constraint handling [18] and HIMPA integration technology to optimize the decision variables in the system. Moreover, the uncertainty in the power system is also modeled using various probability density functions (PDFs) [11]. In addition, 1000 scenarios are generated employing the Monte Carlo Simulation method [19] and then chosen by applying a backward reduction technique. Finally, the feasibility and effectiveness of the HIMPA algorithm are verified on IEEE 30-bus test system.
The remainder of the paper is structured as follows: Section II describes the ORPD problem formulation and ε-constraint handling technique. Section III describes the uncertainties model of RERs and load demand. In Section IV, the background knowledge is briefly recalled. The improved MPA is detailed in Section V. The related experimental analysis and discussion are provided in Section VI. Section VII concludes the work.

II. PROBLEM FORMULATION
The mathematical model of the ORPD problem can generally be defined as follows: where F (x, u) is the objective function, g i (x, u) and h j (x, u) indicate the ith inequality constraint and the jth equality constraint, respectively. The vectors x and u represent control and state variables, respectively, i.e., where V G i represents the voltage at the ith generator bus, Q C j denotes the shunt compensation at the jth bus and T k signifies the kth branch transformer tap; NG, NC, and NT are the number of generators buses, shunt compensators, and transformers, respectively. In this paper, Q C and T k are discrete variables. V L m is the voltage at the mth load bus, Q G i denotes the reactive power at the ith generator bus and S l q represents the line loading at the qth line. NL and nl are the number of load buses and the transmission lines, respectively.

A. OBJECTIVE FUNCTIONS 1) MINIMIZATION OF THE REAL POWER LOSS
The mathematical model with the objective of reducing the real power loss (F Loss ) can be expressed as [20]: where V i is the voltage magnitude at the ith bus. G ij denotes the transfer conductance between bus i and j, and δ i is the voltage angle at the ith bus.

2) MINIMIZATION OF THE VOLTAGE DEVIATION
The voltage deviation is calculated as the sum of all load bus voltages deviating in the network from the ideal value by 1.0 p.u. Mathematically, the voltage deviation objective function (F VD ) can be defined as [20]: B. CONSTRAINTS

1) EQUALITY CONSTRAINTS
The power balance equations are stated as follows: where P G and Q G are the active and reactive power of the generator, respectively. P D and Q D represent the active and reactive load demands, respectively. B ij denotes the susceptance between bus i and j, and NB is the number of buses.

2) INEQUALITY CONSTRAINTS
To ensure a safe and stable operation, control variables and state variables necessarily satisfy the following inequality constraints.
∀k ∈ NT (10) VOLUME 10, 2022 The equality constraints in the ORPD problem can be converted into inequality constraints, and all constraints can be written as: where δ signifies the tolerance parameter. Therefore, the individual's total constraint violation can be formulated as: where w i is a weight parameter and G i,max represents the maximum violation of the constraint G i (x, u) obtained so far. The ε-constraint handling method is an effective constraint handling technique, which is first introduced by Takahama and Sakai [18]. The core idea is to find feasible solutions by controlling the ε-parameter [11]. The traditional ε setting method is stated as follows: where x θ is the θth individual with the highest overall constraint violation in the initial population and θ = 0.05 * N . ε(1) is the initial epsilon parameter and t indicates the current iteration. T c represents the number of generations that the ε-parameter is managed before being set to zero. According to [18], the relevant parameters are set as cp ∈ [2,10], where t max is the maximum number of iterations. We can follow the feasibility rule proposed by Deb to find an optimal feasible solution [21].

A. UNCERTAINTY MODELING OF SOLAR IRRADIANCE
The lognormal PDF is applied to represent the solar irradiance [11] and expressed as follows: where µ s and σ s denote the mean and standard deviation, respectively. In this paper, the 1000 Monte Carlo scenarios of solar irradiance are created using lognormal PDF (µ s = 5.5, σ s = 0.5).
The output power of the photovoltaic unit as a function of irradiance can be calculated as [22]: where P sr is the photovoltaic unit rated power which equals 50 MW; G sd and X c denote the standard solar irradiance and a certain irradiation point which are equal to 1000W /m 2 and 120W /m 2 , respectively [11].

B. UNCERTAINTY MODELING OF LOAD
The uncertainty of load is modeled using the normal PDF defined as follows [23]: where P d represents probability density of normal distribution of load; σ d and µ d are the standard deviation and mean values which are equal to 10 and 70, respectively [11]. Here, the Monte Carlo method is used to establish load demand scenarios (sample size 1000).

C. UNCERTAINTY MODELING OF WIND SPEED
The Weibull PDF is utilized to signify the wind speed [11] and is expressed as follows: where ν w represents the wind speed, and its unit is meters per second. α and β are the scale and shape parameters of the Weibull PDF, respectively. Weibull PDF is employed in the Monte Carlo method to construct the wind speed distribution scenarios (sample size 1000, α = 9, β = 2). The output power of a wind turbine can be formulated as follows [24]: where P ωr is the wind turbine rated output power; ν wr = 16 m/s, ν wi = 3 m/s, and ν wo = 25 m/s are the rated, cutin, and cut-out speeds of the wind turbine, respectively. Note that a wind farm is made up of 25 turbines with a total output of 75 MW.

D. BACKWARD REDUCTION TECHNIQUE
In this paper, the backward reduction algorithm [25] is applied to simplify 1000 Monte Carlo scenarios into 25 representative scenarios sets. Detailed steps of scenario reduction are described in [11]. The selected scenarios, along with their associated parameters and probabilities, are provided in Table 3.

IV. BACKGROUND A. MARINE PREDATORS ALGORITHM
MPA is a novel nature-inspired metaheuristic algorithm that simulates the biological behavior of marine predators foraging for prey [14]. The main operation idea is that predators adopt two random walks calledĹevy flight and Brownian motion to improve the encounter rate with prey and update their position based on the optimal solution. Meanwhile, the prey also acts as a predator. Because the prey is hunting for its own food while the predator is looking for prey. Similar to other population-based metaheuristic algorithms, the MPA also randomly initializes the location of prey through (22): where − → ub and − → lb are the upper and lower boundaries of the search space, respectively. Also, − − → rand represents a random number within the range from 0 to 1. The prey matrix is described as follows: where n is the population size and d signifies problem dimensions. In this step, according to the fitness value of prey, we determine the best solution as the top predator, namely the Elite matrix, which is defined as (24): The entire iterative process is divided into three stages based on different speed ratios while simulating the lifespan of predators and prey. The description of these three stages are as follows [14]: The first stage is carried out in the first third of the iteration when the prey is faster than the predator. The updated process for the prey is as follows: where − − → Step i is the step size vector, and − → R B is a normally distributed random vector, representing the Brownian motion. − → R ∈ [0, 1] denotes a random vector, P = 0.5, and ⊗ indicates the entry-wise multiplication. In addition, t and t max are the present and maximum iterations, respectively.
The second stage is implemented in the middle of the iteration when the predator and prey move at the same speed. The population is split into two groups, one is used for exploration and the other is utilized for development. The mathematical descriptions of both are as follows: where − → R L represents the random vector based on the Lévy distribution.

− − →
Step is an adaptive parameter to control the predator's moving step.
The last stage is executed in the last third of the iteration when the predator outruns the prey. The updated formula for the prey is as follows: However, some studies have shown that eddy formation or the effect of fish aggregating devices (FADs) would affect the foraging behavior of predators, and easily lead to local optimum. This process can be mathematically described as in (33), shown at the bottom of the page, where − → U is a binary vector, FADs = 0.2, and r ∈ [0, 1]. r 1 and r 2 are random indices from the prey matrix.

Spherical search algorithm (SS) is proposed by Abhishek
Kumar to solve global optimization problems with bounded constraints [26]. In the SS algorithm, each solution creates a spherical boundary according to the target direction and makes the trial solution lies on the surface of the spherical boundary. Each spherical boundary is constructed using the axis determined by the individual position and the target position, with the target direction here referring to the primary axis of the spherical boundary.

VOLUME 10, 2022
Similar to other metaheuristic algorithms, the SS algorithm initializes the search process by randomly generating a set of populations, and the expression is shown in (22). Then the ith trial solution − → Y i corresponding to the ith solution − → X i can be calculated as: where c i represents the step-size control vector of the ith trial solution, whose value is randomly calculated within the range of [0.5,0.7] obtained from the experiment. Pr i is the projection matrix, which determines the value of − → Y i on the spherical boundary, and − → Z i denotes the search direction. Calculate the rank of projection using (38) 9.
Update the half population of better solutions using (35) 11. else 12.
Update the other half of the population using (36) 13. endif 14. end for 15. Calculate the trial solution 16. Apply the greedy selection 17. end while The SS algorithm uses two ways to calculate the search direction, namely toward-rand and toward-best. The population is divided into two groups on average according to the fitness value. When determining the search direction, the group with superior fitness uses toward-rand, while the group with poor fitness uses toward-best. The mathematical descriptions of both are as follows: where pi, qi and ri are randomly selected indices from among 1 to n such that pi = qi = ri = i.pbest denotes the randomly selected individuals from the top p solutions in the current iteration.
The projection matrix Pr i in (38) is a symmetric matrix, and its expression is as follows: where A is an orthogonal matrix and − → b i represents a binary vector. The calculation method of binary diagonal matrix diag − → b i is as follows: Followed by using greedy selection to update the ith solution of the population, the mathematical expression can be expressed as follows: The basic steps of the SS algorithm are summarized by the pseudo-code shown in Algorithm 1.

V. HYBRID IMPROVED MARINE PREDATORS ALGORITHM A. OPPOSITION-BASED LEARNING
Opposition-based learning (OBL) is an effective method of population initialization proposed by Tizhoosh [27]. In the metaheuristic algorithm, some randomly generated search agents are frequently distributed far away from the optimal solution or in the opposite position of the optimal solution, which reduces the search efficiency of the population. Here, the OBL strategy improves the quality of the initial population by constructing a solution closer to the optimal solution search region. Suppose − → X i = X i,1 , X i,2 , · · · , X i,d is a solution in d-dimensional space, then the opposite solution − → X i is calculated as follows: Therefore, the fittest search agents are chosen from the population generated by random and its opposite population to construct an initial population.

B. TRANSITIONAL STAGE OF INTEGRATING DIFFERENTIAL EVOLUTION
In the preliminary and intermediate stages of the MPA, predators and prey may skip the most promising area in the search area because they search for food quickly [28]. Inspired by the DE algorithm, the improved DE/best − to-rand/1 strategy is used to replace the middle stage of the traditional MPA, which balances the development and exploration capabilities and improves the convergence speed of the algorithm. The updated mode for the prey position is as follows: where rand represents a random value within [0, 1] . The subscripts r 4 , r 5 and r 6 indicate that integers are randomly selected from 1 to n with r 4 = r 5 = r 6 = i. Note that this phase is executed between t > 1 20 t max and t < 18 20 t max .

C. ADAPTIVE MUTATION OPERATION
In the third stage of iterative optimization, the SS operator is used to strengthen the development stage of the traditional MPA, to improve the search efficiency of the improved algorithm. Here, the probability value of p = 0.7 is utilized to perform the adaptive mutation operation, and the prey updates their position by performing MPA or SS operators to exploit the finest solutions. The corresponding mathematical expression of this strategy is: The basic steps of the proposed HIMPA are summarized by the pseudo-code shown in Algorithm 2.

VI. SIMULATION RESULTS AND DISCUSSION
In this section, the HIMPA is studied for solving the ORPD problem on IEEE 30-bus systems with and without RERs uncertainty. The parameters and generator data of the studied system are taken from [11]. For all given cases, the population size, the number of trial runs, the maximum number of function evaluations, and the maximum number of iterations are specified as 40, 5, 30000, and 374, respectively. The proposed HIMPA runs on MATLAB R2020a with an Intel Core i3 @4 GHz CPU with 64GB RAM.

A. CASE 1: STANDARD IEEE 30-BUS SYSTEM
Cases (1) and (2) aim to minimize real power loss (P loss ) and aggregate voltage deviation (VD), respectively. In both cases, the shunt compensator and transformer taps are set to discrete variables in accordance with the recommendations in [11]. In addition, the active power value of generators should be reasonably selected within the specified range of the generators. To effectively compare with previous studies, the objective functions of Cases (1a) and (2a) are to reduce P loss and VD, respectively. And these cases take into account the active power value of the generator. Table 1 provides the optimization results of all variables of the four study cases. In Case (1), the P loss obtained by the HIMPA is 4.4125 MW, while the VD in Case (2) is 0.08846 p.u. The voltage profiles of load buses are shown in Fig. 1. Voltage of few load buses operates near the maximum limit to achieve the minimum P loss . While Cases (2) (1) and (2), respectively. It can be seen that the MPA-EC is liable to fall into local optimization in the early and middle stages of optimization, and the optimal solution cannot be found quickly. However, the HIMPA can converge to the optimal value quickly in the early stage of optimization. This shows that the improved marine predator algorithm effectively balances the global and local search ability. Table 2 presents a comparison between the proposed algorithm in this study and the reported algorithm. None of the other references in the table specifies the active power settings for the particular generators, except [29] and [30]. Assume that the data in these references are set according to the data in [29]. The study in [11] analyzes in detail the infeasibility of some algorithm solutions in Table 2. In the comparison between Case (2) and (2a), many reported algorithms have better VD values than the proposed HIMPA, but there is no documented study to examine the corresponding actual reactive power generation status or voltage distribution of load buses [11].

B. CASE 2: MODIFIED IEEE 30-BUS SYSTEM
This subsection adopts the stochastic ORPD approach based on wind power uncertainty proposed in [31] and [32] to calculate the expected power loss (EPL) and expected voltage deviation (EVD). When RERs are in short supply, swing generators necessarily provide additional power to meet system demand. The details of stochastic ORPD are explained in [16]. The modified IEEE 30-bus system is to switch the VOLUME 10, 2022        The EPL value of Case (3) is 1.8316 MW, and the EVD value of Case (4) is 0.05856 p.u. According to the data analysis in Table 3, under the condition of minimum network load, P loss is also small (scenario 4). Minimum loading represents the lowest current in the network and therefore the lowest power loss. On the other hand, when the network loading  is the largest, the system active power loss is the largest (scenario 20).
The experimental results of the EVD show that under the minimum load (scenario 4), the bus voltage levels of the network can be maintained close to the desired 1.0 p.u. Similarly, the low current means the minimum aggregate VD. On the contrary, under the maximum load (scenario 20), the aggregate VD of the system is the highest. Fig. 4 depicts the distribution of load bus voltage in these two extreme scenarios. Figs. 5 and 6 show the optimal values of generator bus voltages for each scenario in Cases (3) and (4), respectively. Generator bus voltage values in Case (3) are higher than those in Case (4); particularly in Case (4), bus 11 has a greater voltage setting, since it is connected to the nearby load buses without any shunt compensators.

VII. CONCLUSION
In this paper, a hybrid improved marine predators algorithm (HIMPA) is proposed, which can improve the performance of the original MPA when solving the ORPD problem. The three improvements of the HIMPA are based on opposition-based learning, the transitional stage of integrating differential evolution, and the adaptive mutation strategy. In the first section, the proposed HIMPA is applied to solve the deterministic ORPD problem, and the basic IEEE 30-bus system is selected for testing. The optimization objectives are real power loss and voltage deviation. The test results show that the optimal solution obtained by HIMPA is better than other competitive algorithms. Subsequently, stochastic ORPD solutions with uncertainties of solar power, wind, and load demand have been investigated by a scenario-based approach for the modified IEEE 30-bus system. Various scenarios are created by Monte Carlo simulation, and then a group of representative scenarios is selected by using the backward reduction technique. The expected power loss and expected voltage deviation values in different scenarios are calculated by optimizing these parameters. The results show that the combination of the HIMPA and ε-constraint technique can effectively solve the ORPD problem.