Optimal Power Flow Solution Based on Jellyfish Search Optimization Considering Uncertainty of Renewable Energy Sources

Today’s electrical power system became more complex interconnected network that is expanding every day. The transmission lines of the power system are more severely loaded than ever before. Hence, the power system is facing many problems such as power losses increasing, voltage instability, line overloads, etc. The optimization of real and reactive powers due to the installation of energy resources at appropriate buses can minimize the losses and improve the voltage profile especially, for congested networks. As a result, the optimal power flow problem (OPF) is considered more important tool for the processes of planning and operation of power systems. OPF is a very significant tool for power system operators to meet the electricity demand of the consumers efficiently, and for the reliable operation of the power system. However, the incorporation of renewable energy sources (RESs) into the electrical grid is a very challenging problem due to their intermittent nature. In this paper, the proposed power flow model contains three different types of energy sources: thermal power generators representing the conventional energy sources, wind power generators (WPGs), and solar photovoltaic generators (SPGs) representing RESs. Uncertain output powers from WPGs and SPGs are forecasted with the aid of Weibull and lognormal probability distribution functions (PDF), respectively. The under and overestimation output powers of RESs are taken into consideration while formulating the objective function through adding a penalty and reserve cost, respectively. Moreover, carbon tax is imposed to the main objective function to help in reducing carbon emissions. A jellyfish search optimizer (JS) is employed to reach optimization in the modified IEEE 30-bus test system to validate its feasibility. To examine the effectiveness of the proposed JS algorithm, its simulation results are compared with the results of four other nature-inspired global optimization algorithms. The developed OPF algorithm considers several practical cases such as generation uncertainty of renewable energy sources, time-varying load and the ramp rate limits of thermal generators. The simulation results show the effectiveness of the JS algorithm in solving the OPF problem in terms of minimization of total generation cost and solution convergence.


PDF
probability density function ISO Independent system operator SF Superiority of feasible solution TG Thermal generator WPG Wind power generator SPG Solar photovoltaic generator RES Renewable energy sources a i , b i , c i Cost coefficients of the i -th thermal power generator VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ P TGi Output power of the i-th thermal generator l i and m i Valve-point loading effect coefficients P Ws,j Scheduled wind power of the j-th wind power generator d j Direct cost coefficient of the j-th wind power generator P Ss,k Scheduled solar PV power of the k-th solar PV power plant K RW .j Reserve cost coefficient of the j-th wind power plant P Wa.j Available power from the j-th wind power plant K PW .j Penalty cost coefficient of the j-th wind power plant K RS,k Reserve cost coefficient of the k-th solar PV plant P Sa,k Available power from the k-th solar PV power plant P 0

TGi
Power of i th thermal generator at previous hour K PS,k Penalty cost coefficient of the k-th solar PV power plant C tax Carbon emission tax in $/ton Probability of wind speed following Weibull PDF (m/s) k, c Weibull PDF shape and scale factors f T (T ) probability of solar irradiance (T) following lognormal PDF (W/m 2 ) σ and µ Standard deviation and mean of lognormal PDF P loss Active power losses in the network V d Voltage deviation T Solar irradiance (W/m 2 ) DR i , UR i down and up ramp-rate limits of i-th thermal power generator

I. INTRODUCTION A. LITERATURE REVIEW
In 1962, optimal power flow (OPF) was first formulated by Carpentier [1]. Then several techniques for solving OPF have been proposed. OPF is employed for optimizing generation cost, reducing gas emissions, minimization of power losses, and keeping voltage stability. This optimization is usually forced by the physical limitations of the power system, where power generator capability, transmission line capacity, bus voltage, power cable flows, and any other technical constraints have to be satisfied. This can represent a fancy problem, especially in large power systems. Thus, distinctive care must be taken to guarantee that these technical constraints are not violated. Traditional OPF includes just conventional generation sources operating by firing fossilfuel, and this already leads to an extremely non-linear, mixed integer, and non-convex optimization problem [2]- [5]. With high penetration levels of renewable energy sources (RES) in the power systems, OPF study becomes essential incorporating the uncertainty nature of these sources due to the associated challenges at the planning and operational phases. Since the initiation of OPF, many classical optimization methods have been proposed for solving OPF. These methods include interior-point methods, mixed-integer linear programming, non-linear programming, and quadratic programming [6]- [8]. Some of these methods have been successfully implemented by the industry sector because they have a fast convergence in addition to the strength in obtaining an optimum solution. However, such optimization methods require to linearizing the optimization function firstly. For this issue, some properties are often approximated for the optimization function like the non-convex, non-differentiable, and non-smooth properties.
To provide a solution for this problem, heuristic optimization techniques have also been suggested [9], [10]. In this regard, numerous heuristic techniques are used to solve the OPF.
Ref. [11], has presented a reliable and efficient Tabu search best method, which has been tested on the standard IEEE 30-bus power system, for achieving different objectives functions within several operating constraints.
Simulated annealing that is a single solution-based algorithm can provide a global or near global optimal solution but with high computational time [12].
In [13], a simple genetic algorithm (SGA) has been applied for solving OPF problem, where a sequential GA solution scheme has been employed to achieve suitable control variable resolution without violation of system constraints.
Ref. [14] has provided a hybrid genetic algorithm combining GA with the linear programming (LP) and sequential quadratic programming (SQP) algorithms available in the MATPOWER software.
Ref. [15] has presented a simple refined genetic algorithm (RGA) with the ability to code large number of control variables in a practical power system, within a reasonable length of chromosome.
In [16], a multi-parent crossover based genetic algorithm (GA-MPC) has been presented to solve OPF problem of power system consisting of thermal generators considering various objective functions such as piecewise quadratic cost, total fuel cost, valve-point loading effects and emission.
In [17], an enhanced genetic algorithm has been presented for solving OPF problem. This algorithm has been validated using IEEE 30-bus power system and the three area IEEE RTS-96.
Differential evolution (DE) has been applied to solve OPF problem in many previous works. It is suitable for OPF problems including complex variables or transient stability constraints and it has fast convergence characteristics but the probability of converging to a local, rather than global, optima is high [18]- [22].
Particle swarm optimization (PSO) has been used in several difficult OPF problems. However, as with many heuristic approaches, a primary drawback of traditional PSO is premature convergence when the parameters are not chosen correctly [23]- [28].
In [36], authors provided a modified bacteria foraging algorithm (MBFA) to solve OPF incorporating with a model for doubly fed induction generator (DFIG), thermal power generators and static synchronous compensator (STATCOM).
The authors in Ref. [40] have proposed a new modification on JAYA algorithm (MJAYA) to provide a solution for the premature convergence problem of the original JAYA.
Ref. [41] proposed a global best artificial bee colony algorithm, not only to improve the initialization phase of ABC algorithm, but to enhance the search technique to reach an optimum solution as well.
Ref. [37] proposed a modification on DE called SHADE algorithm where the selection process of future control parameters is guided through the successful control parameters settings to guarantee an appropriate balance concerning the exploration and exploitation phases. Moreover, this helped in achieving comparatively fast convergence rate for OPF problems.
Ref. [42] proposed a combination between SHADE and a technique for handling constraints, called the superiority of feasible solution (SF) which helped in increasing the efficiency of SHADE.

B. CONTRIBUTION AND PAPER ORGANIZATION
This paper comes in the context of completing the exerted efforts to find an optimal solution for the OPF problem. In this paper, a new recent metaheuristic optimization algorithm called jellyfish search (JS) which is developed in [43] is used to efficiently solve the OPF problem. The proposed algorithm is applied to IEEE 30 -bus system incorporating with two wind generators and one solar PV generator to verify its validity in obtaining the optimal solution for OPF problem with renewable energy sources during theoretical and practical conditions. Four other optimization techniques: Giza pyramids construction (GPC) [44], chaos game optimization (CGO) [45], flower pollination algorithm (FPA) [46] and artificial bee colony (ABC) [29] are applied to the proposed system to compare their results with the results of JS. The results of JS are also compared with the results obtained by SHADE-SF in [42].
The structure of the rest of this paper is prepared as follows. Section II shows the mathematical model and associated applicable constraints applied for the OPF problem. Section III introduces the uncertainty output models of WPG and SPG. In section IV, the new proposal for applying JS to OPF including the uncertain RES is developed. Section V provides simulation results of numerous realistic case studies for the five algorithms under study. Finally, section VI presents conclusions for this paper. Table 1 summarizes the parameters of IEEE 30-bus system under work. The modified network involves three different types of power generation resources i.e. thermal power generators (TGs) with constant outputs, solar PV generator (SPG) with variable output, and wind generators (WPG) with variable outputs. This variation in PV and wind outputs must be balanced through the mixture of all generators and reserve power, so the overall generation cost comprises costs of operation for all generators, reserve cost and penalty cost.

A. COST MODEL FOR THERMAL POWER GENERATORS
Thermal power generators use the fossil fuel to operate. The relation between the fossil fuel cost in $/hr and thermal generator output power in MW is given by (1).
where a i , b i , c i represent the cost coefficients of the i -th thermal power generator with power output P TGi . N TG is the total number of thermal power generators.
However, taking the valve point loading effect in our consideration gives more precise and realistic cost, so the modified quadratic relationship with a valve point loading effect will be as following: where, l i and m i represent coefficients of the valve-point loading effect. P min TGi is the minimum power of i-th thermal unit generator. Table 2 provides all cost and emission coefficients used in calculations related to the thermal generating units.

B. DIRECT COST OF POWER PRODUCED FROM WIND AND SOLAR PHOTOVOLTAIC SOURCES
WPGs and SPG do not need fossil fuel to operate. They depend on wind and the sun to generate power. When an independent system operator (ISO) owns the solar PV/wind generators, the cost function might not be as the solar PV/wind generators need no fuel, unless ISO desires to allocate some reimbursement cost to the initial expenditure  [42].

TABLE 2.
Coefficients of emission and cost for thermal power generators in the system under study [48].
for these plants or to allocate this in a form of the cost required for maintenance and renewal [42]. On another side, when private parties own solar PV or wind generators, ISO is obliged to pay a price in proportion to the scheduled power contracted.
The direct cost of the j-th wind power generator as a function of scheduled power is molded as following, where P Ws,j and d j are the scheduled wind power and the direct cost coefficient related to the j-th wind power plant respectively. In a same way, the direct cost of the k-th solar PV power generator is: where P Ss,k and e k are the scheduled solar PV power and the direct cost coefficient related to the k-th solar PV power plant respectively.

C. EVALUATION OF COST DUE TO UNCERTAINTIES IN WIND POWER
There are two scenarios may occur due to the intermittent nature of wind energy. The first scenario takes place when the output power of the wind plant is less than the expected output power. This situation is termed as overestimation of output power. In this situation, the system operator uses spinning reserve to provide reliable power supply to its consumers. The cost required for committing the reserve generating units to overcome the overestimation situation is referred to reserve cost.
Reserve cost for the j-th wind power plant is determined by: where K RW ,j and P Wa,j are the reserve cost coefficient and the available power relating to the j-th wind power plant respectively, and f W (P W ,j ) is referred as the PDF for the output power of the j-th wind power farm. The second scenario occurs when the output power of the wind power plant is higher than the anticipated value of the output power. This situation is expressed as underestimation of output power. Consequently, the remaining power will be lost if it is not possible to consume through decreasing the output power of conventional generators. In this situation ISO has to compensate a penalty cost associated to the remaining power.
Penalty cost of wind power farm j is determined by: where K PW ,j and P Wr,j are the penalty cost coefficient and the rated power relating to the j-th wind power plant respectively.

D. EVALUATION OF COST DUE TO UNCERTAINTIES IN SOLAR PHOTOVOLTAIC POWER
Similar to wind energy, solar energy has also uncertain and intermittent output. The approach used to solve under and overestimation of solar output power should be in principle the same as the wind output power.
However, for more simplicity in calculations, penalty and reserve cost models are formed following the proposed concept in Ref. [49].
This concept is used because solar radiation follows lognormal PDF [50], which represents a meaningful difference comparing to the Weibull PDF used for wind speed distribution.
More details are provided in Section III for calculating stochastic solar PV and wind power.
The reserve cost of solar PV power plant k is given by: where K RS,k represents the coefficient of reserve cost and P Sa,k is the available power related to solar PV power plant k. f s (P Sa,k < P Ss,k ) represents the solar power shortage probability, while E(P Sa,k < P Ss,k ) referrers to the expected output of the solar PV power generator below P Ss,k . The penalty cost of solar PV power generator k is given by: where K PS,k is the penalty cost and f s (P Sa,k > P Ss,k ) represents the probability of remaining power generated by the solar PV power plant k comparing to P Ss,k , while E(P Sa,k < P Ss,k ) referrers to the expected remaining output power.

E. CARBON EMISSIONS AND TAX
Using conventional energy sources to produce power is well known as a source of emitting greenhouse gases into the environment. Two of these harmful gases are SO x and NO x which increase in emitting with growth in produced power from conventional power generators. The relationship between emission in tonnes per hour (t/h) and produced power (in p.u. MW) is represented through Eq. (9).
where α i , β i , γ i , ω i and µ i represent emission coefficients associated to the thermal power generator i. Values of emission coefficients related to the thermal power generator are illustrated in Table 2. Which are similar to values mentioned in [48] with some slight change in the value of µ for the thermal power generator linked to bus 1.
In current years, owing to the harmful impact of climate change, energy regulators in several countries are setting huge regulations on the whole energy industry sector to decrease the carbon emission [51]. To encourage investors to invest in clean and renewable energy sources like solar and wind, carbon tax (C tax ) is enforced per ton of greenhouse gasses. The cost of carbon emission in $/h is calculated as: The objective optimization for OPF is formed including all models of cost functions as illustrated in Eqs. 2-8.
In objective function (F1), emission cost is neglected. To clarify the variation in the scheduling of generation when considering a carbon tax, the objective function (F2) is formulated including emission cost as illustrated in Eqs.9-10.
Consequently, the first objective function is to minimize: (11) where, N WG and N SG represent the number of wind power generators and solar PV generators existing in the system respectively.
The second objective function is to minimize: The above-mentioned OPF objective functions are exposed to both inequality and equality constraints of the system as illustrated below.

1) EQUALITY CONSTRAINTS
Equality constraints are necessary load flow equations for providing power balance of both active and reactive powers generated in the system. This mean that these powers must be equal to the total demand and total losses in the system. The equality constraints are represented as follows [42]: where δ ij = (δ i− δ j ) is the voltage angles difference between buses i and j, NB refers to the total number of network buses, P Di and Q Di represent active and reactive components of load demand connected to bus i respectively whilst the active and reactive power components of generation at bus i are represented by P Gi and Q Gi respectively from any of the energy resources either conventional power plants or renewable energy plants. G ij and B ij are the transfer conductance and susceptance is between buses i and j respectively.

2) INEQUALITY CONSTRAINTS
The inequality constraints represent the operating limits of components and equipment in power system. These constraints include also security constraints on transmission lines and load buses.
(a) Generator's constraints: (b) Security constraints: The generation limits of active power for thermal generators, wind generators, and solar PV generators are represented by Eqs. (15) - (17) respectively. With the same arrangement, Eqs. From (18) to (20) represent reactive power limits of all generators. While N G represents the total number of generator buses or generators. Eq. (21) defines constraints on generator buses voltage, and Eq. (22) provides voltage limits subjected to PQ buses (load buses) and N L refers to the number of PQ buses. Eq. (23) defines the constraints on the transmission Lines capacity for total N L numbers of transmission lines within the system.
It is important to note that power flow convergence to an optimal solution guarantees automatically satisfying the imposed equality constraints through equations of power balance. Amongst inequality constraints, voltages of generator buses and active powers of generators (excluding swing or slack generator assumed to be linked to bus 1) are characterized as self-limiting control variables.
For each variable from the control variables, the applied optimization technique chooses a reasonable value lied in the accepted range of this variable. Thus, inequality constraints related to the active and reactive power outputs of the slack generator, reactive power outputs of other generators, PQ buses' voltage limits, and transmission line capabilities require special attention.

3) PRACTICAL CONSTRAINTS
Thermal power generators are regulated to operate always between two adjacent limits.
Generally, thermal power generators are connected with steam valves to decrease or increase their power generation. In real-time, due to physical limitations, the power generated from a thermal generator tracks ramp-functions i.e. to decrease the produced power, the generator follows down ramp-rate, and to increase the produced power, the generator follows up ramp-rate. The ramp rate limits can be mathematically formed as: (24) where, P 0 TGi represents the i th generator power generation at previous hour. DR i and UR i represent the respective down and up ramp-rate limits of i-th thermal power generator, respectively.
The ramp-rate constraints can be given in OPF problem as follows:

generation increases
In OPF problem, due to the existence of ramp-rate limits, the generation fuel cost increases because the operating points of the thermal generating units are changed from its operating limits. The new upper and lower limits of a thermal generator considering ramp-rate limits must not exceed the following limits: P high TGi and P low TGi are the new limits of thermal unit-i.

G. LOAD BUS MODELLING
The reactive power capability of the generator is a vital issue in the OPF study. In this work, narrower ranges are applied for thermal power generators than what have been applied in [47]. In recent years, the reactive power capability of wind power plants has significantly developed. Commercially, wind turbines (WTs) including full reactive power capability and other important features are already available [52].
With the assistance of the Enercon FACTS -WT reactive capability curve, it is clear that a wind turbine can deliver reactive power from −0.4 p.u. to 0.5 p.u. throughout its range of output active power. Delivering negative reactive power indicates the capability of the generator to absorb reactive power.
A rooftop solar PV system can be characterized as PQ bus (load bus) with reactive power equal to zero (Q = 0). While large scale solar PV systems are connected with inverters. As a result, a full capability model of solar PV generator is necessary due to the dynamic performance of the inverters [53].
In this aspect, ref. [54] provides converter and controller models when implementing a detailed study of reactive power capabilities of solar PV generator. In ref. [55], the PV converter capability study has been extended taking into consideration the impact of change in solar irradiance and surrounding temperature. In this work, the reactive power capability chart of solar PV generator is nearly considered from −0.4 p.u. to 0.5 p.u. All limits of active and reactive powers of the power generators in the system under study are illustrated in Table 5 in addition to simulation results. Parameters of the system like voltage deviation and active power loss in the transmission lines are also important in OPF problem study. The active power loss in transmission lines is unavoidable due to the inherent resistance of the lines. The real power loss of the network is calculated by: where δ ij represents the difference of voltage angles between buses i and j and G ij refers to the transfer conductance. The voltage deviation indicator represents the cumulative deviation of voltages of all PQ buses (load buses) in the system from a nominal value (1 p.u.). It provides an indicator for the quality of power system voltage and mathematically expressed as:

III. STOCHASTIC WIND POWER, SOLAR POWER AND UNCERTAINTY MODELS
Refs. [36] and [49] show that wind speed distribution follows Weibull probability density function (PDF) for mean power calculations of wind turbines. The probability of wind speed according Weibull PDF is given by: where v is wind speed in m/s, k is the shape factor and c is scale factor. Eq. 30 defines the mean of Weibull distribution: where is gamma function that is defined as: The IEEE-30 bus system is adjusted in this study by replacing conventional generators connected to buses 5 and 11 with two wind power plants. Table 3 provides the selected values of Weibull scale (c) and shape (k) parameters. These values have been followed in this study except other values used for a specific case study. Figs. 1 and 2 provide Weibull fitting and distributions of wind frequency. They are got after running Monte Carlo simulation with 8000 iterations [42]. Ref. [56] regulates the requirement of wind turbine design and states the maximum turbulent class IA of the wind turbine which is suitable for operation at the maximum annual average wind speed (10 m/s at hub height). Shape (k) and scale (c) parameters for the wind farms are carefully selected to remain the maximum value of Weibull mean around 10 to guarantee both realistic and diverse geographic locations for the two wind farms.
In the same way, we replaced the conventional generator connected to bus 13 of the IEEE-30 bus system with solar PV generator. The output of solar PV generator depends on the solar irradiance (T) which follows lognormal PDF [50]. The probability of solar irradiance (T) following lognormal PDF is given by: where σ and µ are the standard deviation and mean of lognormal PDF respectively. Fig. 3 shows lognormal PDF of solar irradiance after Monte Carlo simulation with 8000 iterations.    Table 3 provides the chosen values for lognormal PDF parameters. These values have been followed in this study except other values used for a specific case study.

A. WIND POWER AND SOLAR PHOTOVOLTAIC POWER MODELS
As mentioned above, the modified IEEE-30 bus system has two wind power farms. Wind farm 1 is connected to bus 5, it consists of 25 wind turbines with a rated power of 3 MW for each turbine, and thus the cumulative active power of wind farm 1 is 75 MW. Similarly, wind farm 2 is connected to bus 11, it contains 20 wind turbines, and each turbine has a rated power of 3 MW, thus the cumulative active power of wind farm 2 is 60 MW.
The actual output power of wind power generator depends on the wind speed. A wind turbine output power is given by [42]: where v r , v out , and v in are the rated, cut-out, and cut-in wind speeds of the wind power turbine respectively. P wr represents the rated output power of the wind turbine. Based on the product datasheet of Enercon E82-E4 wind turbine, various speeds of a 3-MW wind turbine are v r s = 16, m/s v out = 25 m/s, and v in = 3 m/s. Similarly, the energy conversion for solar PV in terms of solar irradiance (T) is expressed by [57]: where T std is the solar irradiance set as 800 W/m 2 . in standard environment, R c represents a certain irradiance point fixed as 120 W/m 2 . and P Sr is the rated output power of the solar PV unit.

B. WIND POWER PROBABILITY MODEL
Based on (34), the variable output power of wind generators is discrete in some regions of wind speed. The output power is zero when wind speed (v) is lower than cut-in speed (v in ) and higher than cut-out speed (v out ), while the output power will be at its rated value P wr when the wind speed lies between rated wind speed (v r ) and cut-out speed (v out ). The wind output power probabilities for these discrete zones are described by [58]: The output power of wind turbine is continuous between cut-in speed (v in ) and rated speed (v r ) of wind. The probability of wind output power in the continuous zone are described by [13]:

C. SOLAR PV POWER OVER/UNDERESTIMATION COST
The stochastic output power from the solar PV plant is represented by the histogram in Fig. 4. The dotted line illustrates the scheduled power supposed to be delivered to the network by the solar PV generator. It is necessary to observe that scheduled solar PV power is a variable value, so there is a jointly power agreement between the ISO and the owner of the solar PV generator. Eqs. (39) and (40) are used in the model to calculate over and underestimation cost of the solar PV generator respectively: (40) where P Sn− and P Sn+ are the shortage power and surplus power, as lying on the right and left half plane of schedule power P Ss in the histogram of Fig. 4. Likewise, f Sn− and f Sn+ are relative frequencies for the occurrence of P Sn− and P Sn+ . N − and N + represent number of discrete bins on the left and right planes of P Ss for PDF generation.

IV. OPTIMISATION TECHNIQUE
In 2020, a novel metaheuristic optimizer titled artificial jellyfish search optimizer has been proposed by J.-S. Chou and D.-N. Truong [43]. The behavior of jellyfish for finding food in the ocean inspired the JS algorithm. At the beginning, they follow ocean current, then they move inside swarms as time passes on, in the presence of a mechanism to control the time of switching amongst these motions.

A. MATHEMATICAL MODEL
The proposed optimization technique depending on three main rules. 1. Jellyfish either pursue the ocean currents or proceed inside swarms, and a time control mechanism regulates the switching process between these motions. 2. Jellyfish proceed in the ocean seeking food. They are highly attracted to places that have a greater amount of food. 3. The amount of food found by jellyfish is determined through the position and its correspondent objective function.

1) CURRENT OF OCEAN
The ocean current includes huge quantities of nutrients that attract jellyfish. Eq. (41) determines the ocean current direction through taking the average of the entirely vectors from every jellyfish in the ocean to that is presently in the best position. Set Therefore trend → is determined by: where N Pop represents jellyfish number, X * represents the jellyfish presently in the best position in the jellyfish swarm, while e c refers to the control factor that manages the attraction process, the mean position of whole jellyfish in the ocean is represented by µ s , df refers to the difference between the current best position of the jellyfish and the mean position of all jellyfish. Assume that jellyfish distribute with a normal spatial distribution in all dimensions, thus there will be a distance (±β s σ s ) nearby the mean position includes a certain probability of all jellyfish.
where σ s refers to the standard deviation of the normal spatial distribution. Therefore, Set Hence, df = β s * rand f (0, 1) * rand α s (0, 1) * µ s To make it simpler, Eq. 43 is formed as follows: where e c = β s * rand(0, 1) Thus, So, we can define the new position of each jellyfish by: Eq. 47 can be written as: where β s > 0 refers to the distribution coefficient, associated with the Trend → length. In this case, β s =3.

2) SWARM OF JELLYFISH
In a swarm, there are two types of motions: type A (passive motion) and type B (active motion). In the beginning, when the jellyfish swarm has just been shaped, most of jellyfish show passive motion. As time passes on, they increasingly show active motion. In passive motion, jellyfish move around their original locations. Consequently, each jellyfish updates its location by (52).
where U B and L B represent the upper and lower bounds of the search space, respectively. γ s > 0 represents the motion coefficient, associated with the motion length around locations of jellyfish. In this paper, γ s = 0.1.
To simulate the active motion (type B), we will randomly select a jellyfish (Jf j ) different to the jellyfish of interest (Jf i ). Consequently, to determine the direction of motion, a vector from Jf i to the Jf j is utilized to define the direction of motion. When the amount of food founded at the position of Jf j is greater than the amount at the position of Jf i , Jf i moves in the direction of Jf j . In contrast, if the available amount of food to Jf j is lesser than the available amount to Jf i , Jf i moves directly in opposite direction of Jf j . Thus, each jellyfish in a swarm moves directly to the better trend to find food. The motion direction and the updated location of a jellyfish are simulated by (55) and (56), respectively. This movement represents an effective way to exploit the search space. Step where Step → = rand(0, 1) * Direction → (54) where, f (X ) represents the objective function and X represents the location of jellyfish. Therefore, A mechanism for time control is set to define the motion type over time. Its function is to control not only passive and active motions in a jellyfish swarm but also the movements in the direction of the ocean current.

3) MECHANISM OF TIME CONTROL
Ocean currents contain large quantities of food that attract jellyfish. As time passes on, a swarm of jellyfish is formed through gathering more jellyfish together. When the ocean current changes by wind or temperature, the jellyfish existing in the formed swarm move in the direction of another suitable ocean current and consequently another swarm of jellyfish is formed.
As mentioned above, there are two types of jellyfish motion inside a swarm: passive motions (type A) and active motion (type B), the motion of jellyfish is switched between them. In the beginning, jellyfish prefer the motion of type A; as time passes on, the motion of type B is favored.
The mechanism of time control is presented to provide a simulation for this situation. The time control mechanism regulates the switching process between the two types of motion by determining the time control function C(t) represented by (57). The value of C(t) fluctuates from 0 to 1 overtime. This value is compared with a constant C o ., when it is greater than C o , the jellyfish follow the ocean current. When its value is lower than C o , jellyfish tend to move inside the swarm. The value of C o is assumed to 0.5.
where t refers to the iteration number, while Max iter represents the maximum number of iterations.

4) INITIALIZATION OF JELLYFISH POPULATION
In normal, Jellyfish population is randomly initialized. This population approach has some disadvantages such as trapping at local optimal solution and the slow convergence. There are many chaotic maps have been advanced to enhance the initial population diversity including tent map, Liebovitch map, and logistic map. In this study, the logistic map is used because it helps in providing more diversity in initial populations than the diversity of random initialization. It also reduces the premature convergence probability.
where, X i refers to the logistic chaotic value of i th jellyfish location, η is a parameter set to 4.0, and X o is used to generate the initial population of jellyfish, where X o ∈ (0, 1), X o / ∈ {0.0, 0.25, 0.75, 0.5, 1.0}.

5) BOUNDARY CONDITIONS
Depending on the fact of spherically of the earth and oceans are found around it, when a jellyfish leaves the bounded of search area, it will automatically go back to the opposite bound of the area. This re-entering procedure is presented by (59).
where X i,d represents the i th jellyfish location in a dimension (d th ), X i,d refers to the updated location of the i th jellyfish after checking the boundary constraints. While L b,d and U b,d are lower and upper bounds of the d th dimension in food search spaces, respectively.

6) SCHEMATIC REPRESENTATION OF JS ALGORITHM
Exploration and exploitation phases represent the two key phases of a metaheuristic optimization algorithm. In the JS optimization algorithm, the motion of jellyfish in the direction of an ocean current represents the exploration phase and the motion inside a jellyfish swarm represents the exploitation phase, and the switching process between them is implemented using a time control mechanism. In the beginning, the exploration probability is higher than the exploitation probability because in this time jellyfish try to discover areas that includes promising optimal locations; as time passes on, the exploitation probability exceeds the exploration probability, and hence the jellyfish recognize the best position inside the discovered areas. Table 4 and Fig.5 illustrate the algorithm and flowchart of the JS optimization algorithm, respectively.

V. CASE STUDIES AND SIMULATION RESULTS
In this section, different case studies are executed for the modified IEEE-30 bus system. Results of the various case studies using the proposed Jellyfish search (JS) optimization algorithm are presented and explained. To verify the validity of JS optimization algorithm, four different optimization algorithms: ABC, CGO, FPA, and GPC are applied. The results obtained by JS optimization algorithm are compared with the results obtained by the four above mentioned optimization algorithms and the results obtained by SHADE-SF in [42].
Case studies one and two are applied to study the impact of change of schedule wind and solar powers and PDF parameters on the generation costs. In the remaining case studies, the aim is to optimize the schedule generation from all sources under theoretical and practical conditions. MATLAB software is used to perform the simulation of optimization techniques with computer properties: Intel (R) Core (TM) i3 CPU M380 @2.53 GHz and installed memory (RAM) 4 GB.
After running the algorithm five times, the optimal value of the objective function of each case study is found and control variables settings are listed, a maximum of 1000 iterations are implemented as the ending condition.

A. CASE-1: SCHEDULED OUTPUT POWER OF WIND AND SOLAR POWER PLANTS AGAINST DIFFERENT COSTS
In this case, the parameters of Weibull PDF are the same as provided in Table 3. Relevant parameters of wind turbine are presented in section III.A. The coefficients of direct cost for the wind power are g 1 =1.6 and g 2 = 1.75. The coefficients of penalty cost for not fully using wind power are K PW ,1 = K PW ,2 = 1.5 and the coefficients of reserve cost corresponding to overestimation of wind power are assumed K RW ,1 = K RW ,2 = 3. Noteworthy, the direct cost of wind and solar powers is lower than the average cost value of thermal power and the penalty cost for not utilizing wind power is lower than the direct cost [59]. For wind farm, the scheduled power varies from 0 to rated power and Figs. 6 and 7 show variations of total, direct, reserve, and penalty costs for the two wind farms. The total cost refers to the summation of direct, penalty and reserve costs in proportional to the scheduled wind power. In a similar way, Fig. 8 shows variations of total, direct, reserve, and penalty costs for solar PV power plant corresponding to the schedule solar power.
Based on Ref. [60], the yearly cost needed for operating and maintenance of a solar PV power plant lies almost in a range similar to that of the onshore type of wind power plants. Thus, in this work, the direct, reserve and penalty cost coefficients related to the solar PV power plant are assumed to be h = 1.6, K RS = 3 and K PS = 1.5 respectively. Any other parameters related to the solar PV power plant are provided in section III.A. It is worth to be mentioned, the total cost for solar power does not increase uniformly with increasing of the scheduled power of the solar PV power plant with the selected lognormal PDF parameters. Actually, the minimum total cost of solar PV power plant is achieved at around 20 MW of the scheduled power.

B. CASE 2: PROBABILITY DENSITY FUNCTION PARAMETER OF WIND AND SOLAR POWER PLANTS AGAINST DIFFERENT COSTS
This case is presented to study the change in wind power costs for a fixed arbitrarily schedule power with the variation of Weibull distribution scale parameter (c) (with fixed shape parameter for the two wind farms, k = 2). The schedule powers for WG1 and WG2 are fixed at one third of their installed capacity, i.e. 25MW and 20 MW, respectively. This assumption is reasonable since a practical wind power plant's capacity factor lies in a range from 30% to 45% [37]. Cost coefficients for the two wind farms are the same as used in case 1. Figs. 9 and 10 show variations of costs with the variation of the Weibull scale parameter for wind power plant 1 and wind power plant 2 respectively.
The minimum total cost for wind power is achieved at the middle value of the selected range for the scale parameter. With the increase of the scale parameter value, there will VOLUME 9, 2021  be a domination of wind speeds with specific probabilities. With keeping the scheduled power the same, the penalty cost increases leading to raise the total cost. However, the level of reduction in the reserve cost is insignificant above a certain value of the Weibull scale parameter.
For evaluating the variation in the cost of solar PV power with the change in the mean of the lognormal PDF (µ), the value of µ is changed with an increase value of 0.5 in a range from 2 to 7. The scheduled power of solar PV plant is constant at 20 MW with standard deviation (σ ) = 0.6. Cost coefficients for the solar PV plant are the same as used in case 1. Fig. 11 illustrates cost curves for solar power. From Fig.11, it is clear that the total cost of solar power is gradually reducing to its minimum value at lognormal mean  (µ) = 5.5. It is also clear that the reserve cost and penalty cost are the same at about lognormal mean = 5.8. After this value, the penalty cost increase sharply, leading to suddenly increasing in the total cost of solar PV power.
It is necessary to select an appropriate value of the scheduled power of solar PV due to the high sensitivity of solar irradiance for the value of µ and consequently the solar power. When the value of µ is low, solar irradiance is low and consequently the output power is also low, thus approximately the full reserve power is necessary to compensate the output power.
However, if the value of µ is high, the solar irradiance will be high and consequently the solar PV output power will be high.

C. CASE-3: TOTAL GENERATION COST MINIMIZATION
Case-3 is implemented for optimization the scheduled power of both conventional and renewable power generators to minimize the total generation cost depending on Eq.11. A comparison between the convergences of different OPF optimization methods is provided by Fig. 12. Table 5 provides the optimal results of total generation cost, reactive power (Q), control variables, and other important calculated parameters. For clarification, V i refers to the i-th bus voltage, P Ws,1 and P Ws,2 represent the scheduled output power from wind farm 1 and wind farm 2, respectively, and P Ss represents the solar PV schedule power. P loss and V d refer to power loss and voltage deviation that are calculated by Eqs. 24-25, respectively. Simulation results of case 3 show the effectiveness of JS algorithm, fast convergence, and high solution quality comparing to the other OPF optimization algorithms. The minimum total generation cost reached by JS is 781.64. Consequently, for this case, JS exceeds SHADE-SF [42] in addition to all other applied optimization algorithms regarding minimization of total generation cost and solution convergence.

D. CASE-4: MINIMIZATION OF TOTAL GENERATION COST WITH CARBON EMISSION TAX
This case study aims at minimizing the total cost of generation considering the carbon tax (Ct) imposed on the conventional power sources due to their emissions of CO 2 based on Eq. 12. For this case, the imposed carbon tax is assumed to be $20/ton [42]. It is expected with the existence of the carbon tax, increasing the level of penetration from renewable energy sources, and this concept is clear from simulation results. The level of penetration of solar PV and wind powers in the schedule of optimum generation is solely depending on the volume of carbon emissions and the value of imposed carbon tax. For this case, Fig. 13 provides a comparison between the convergence of JS optimization technique and other applied techniques.
For this case study, Table 6 provides the optimal results of total generation cost, reactive power (Q), control variables, and other important calculated parameters.  The minimum generation cost achieved by JS is 810.1201. It is also clear that, for this case, JS exceeds SHADE-SF [42] in addition to all other applied algorithms regarding minimization of total cost and solution convergence.
The voltages of load buses (PQ buses) represent very important factor in OPF. Thus, it is necessary to critically addressee the constraints related to voltages of load buses. In IEEE-30 bus system, operating voltages of all load buses must be inside the range from 0.95 to 1.05 p.u.
In this aspect, Fig. 14 presents voltage profiles of PQ buses for Case-3 in addition to Case-4. Fig.14 shows that all PQ buses voltages are within limits after optimization for the two case studies.

E. CASE-5: OPTIMIZED COST AGAINST RESERVE COST
In this case, all parameters are keeping as in Case 3 excluding coefficients of reserve cost. Coefficients of reserve cost for both solar PV and wind powers are changed from K RS = K RW ,1 = K RW ,2 = K R = 4 to K R = 6 with incremental step of one. The coefficient of penalty cost for both solar PV and wind powers is K PS = K PW ,1 = K PW ,2 = K P = 1.5 as used in both Case 1 and Case 3. Fig. 15 indicates the optimized schedule power of all generators in a bar chart form. For clarification, Case 5a represents the case when K R = 4, Case 5b represents K R = 5, and Case 5c represents K R = 6. When the coefficient of reserve cost increases, the optimal scheduled power from wind and solar PV plants reduces as decreasing the scheduled output power needs a lesser amount of spinning reserve. Thermal power generators compensate the shortage of output power from wind and solar PV sources. Hence, the cost of thermal power generator increases as noted from the profile of thermal power generator cost (TPG) in Fig. 16. Costs of solar PV generator (SPG) and wind power generator (WPG) gradually reduces to an amount. WPG cost comprises cost of output power from the two wind generators. Total power cost increases with the increase in the coefficient of reserve cost. In this case, all parameters are keeping as in Case 3 excluding coefficients of penalty cost. Coefficients of the penalty cost for both solar and wind power are increased from K PS = K PW ,1 = K PW ,2 = K P = 1.5 to K P = 3 for Case 6a, K P = 4 for Case 6b and K P = 5 for Case 6c. The coefficient of reserve cost for both solar PV and wind powers is K RS = K RW ,1 = K RW ,2 = K R = 3 as used in both Case 1 and Case 3. Similar to Case 5, Fig. 17 presents the optimized schedule power of all generators in a bar chart form. When the coefficient of penalty cost increases, the optimal scheduled power from wind and solar PV plants increases as increasing the scheduled output power helps in reducing the penalty cost when solar irradiance or wind speed is high. However, the manner, in this case, is different from Case 5 as the increase of scheduled power seems not to uniform for all renewable energy sources. This different can be interpreted by the extremely non-linear relations between the probability distribution functions and reserve or penalty cost of both solar PV and wind power. Fig. 18. shows a progressive increase in wind power generator cost (WPG cost) with a small fluctuation in the cost of solar PV power (SPG cost) due to the fluctuating of its scheduled output power. It is also noted that the cost of thermal power generation (TPG cost) is nearly constant and a steady rise in total cost is observed.      Case 5 and Case 6. The voltages of all generator buses are inside the stated range (0.95 -1.10 p.u.). Fig. 20 illustrates the schedule reactive power from all generators. By observing the limits of reactive power (Q) provided in Table 5, generators TG3 and WG2 operate at their maximum limits of Q capability for several cases. So, it is necessary to consider the constraints on Q during implementation of any optimization algorithm.

G. CASE-7: MINIMIZATION OF TOTAL GENERATION COST WITH RAMP RATE LIMITS OF THERMAL GENERATORS
Case-7 is implemented to optimize the scheduled power of both conventional and renewable power generators to minimize the total generation cost depending on (11) in more practical conditions through considering the effect of ramp rate limits of thermal power generators.
The power generation at the previous hour for the thermal power generators and their down and up ramp-rate limits are presented in Table 2 [61].
For this case, the optimal results of total generation cost, reactive power (Q), control variables, and other important calculated parameters are recorded in Table 7.
A comparison between the convergence characteristics of different OPF optimization methods used in this case study is provided by Fig. 21.
Comparing the results between Table 5 (Case 3) and Table 7 (Case 7), the increase in the total generation cost is clear as expected due to considering the ramp rate limits of thermal generators. The results prove also the effectiveness of the JS algorithm, fast convergence, and high solution quality compared with the other OPF optimization algorithms. The minimum total generation cost reached by JS is 804.10927. Consequently, for this practical case, JS exceeds SHADE-SF in addition to all other applied optimization algorithms regarding minimization of total generation cost and solution convergence. Table 7 demonstrates that all algorithms duly satisfy all system constraints. Fig.22 shows the PQ buses voltage profile for Case-7. From this figure, it can be observed that all PQ buses voltages are within limits for this case study.

H. CASE-8: CASE STUDY CONSIDERING UNCERTAINTIES IN LOAD DEMAND
In this case, a realistic condition of variable load demand is presented. A normal probability density function (PDF)  is used to model the uncertainty in load demand of the system [62]. A scenario-based methodology is implemented to achieve the optimization process at some discrete scenarios (levels) of load demand.   70 and 10, respectively. Fig. 23 includes four-color areas that represent the four considered different scenarios (levels) of system loading (P d ) in this case. The probability of occurrence of a certain loading scenario and its mean are given as [62].  (61) where, P D,i is the mean of i-th level of loading.    Table 8 provides the calculated means (as a percentage of nominal loading i.e. 100% of the network loading, P d ) and probabilities for the four loading scenarios. Under each scenario, the six applied algorithms optimize the objective function of total generation cost depending on (11). The   scheduled power from all the generators of the system is optimized in each scenario.
Comparing the results given Table 5 (Case 3) and Tables 9-12 (Case 8), it is evident that the network operation cost and power loss are justifiably lower under realistic VOLUME 9, 2021 loading scenarios. Fig. 24 shows that JS algorithm outperformed all other applied techniques in terms of minimization of total generation cost and solution convergence.
The bus voltage profiles for all the scenarios using JS algorithm are drawn in Fig. 25. The profiles and the state variables listed in Tables 9-12 clearly demonstrate that the JS algorithm duly satisfies all system constraints.

VI. CONCLUSION
In this paper, a new metaheuristic optimization algorithm based on jellyfish search optimizer has been proposed for providing an optimal solution of the OPF problem incorporating with stochastic wind and solar energy sources in the IEEE-30 bus power system. The uncertainty nature of wind and solar energy sources has been modelled with the help of Weibull and lognormal probability density functions (PDFs), respectively. Total generation cost including all generation sources is optimized with and without the imposition of carbon emission tax and the impact of changing cost coefficients of wind and solar energy sources on the generation cost is also studied.
Some practical cases such as the uncertainty of load demand and considering the ramp rate limits of thermal generators have been studied to prove the validity of the proposed algorithm during the practical cases.
To verify the validity of JS algorithm, four recent optimization algorithms: ABC, CGO, FPA, and GPC are applied for the modified IEEE-30 bus system. The simulation results of JS algorithm are compared with the results obtained from the four optimization techniques and the results of another recent optimization algorithm (SHADE-SF) provided in literature.
Simulation results show the effectiveness of the JS algorithm in solving the OPF problem as it exceeds the other optimization algorithms regarding generation cost minimization and solution convergence during the theoretical and practical conditions.
The simulation results of JS and other optimization algorithms show also that the physical and security constraints are within the limits predefined by the system operator.