Performance Improvement of Multiobjective Optimal Power Flow-Based Renewable Energy Sources Using Intelligent Algorithm

Producing energy from a variety of sources in a power system requires an optimal schedule to operate the power grids economically and efficiently. Nowadays, power grids might include thermal generators and renewable energy sources (RES). The integration of RES adds complexity to the optimal power flow problem due to intermittence and uncertainty. The study suggests a Multi-Objective Search Group Algorithm (MOSGA) to deal with multi-objective optimal power flow integrated with a stochastic wind and solar powers (MOOPF-WS) problem. Weibull and lognormal probability distribution functions (PDFs) are respectively adopted to describe uncertainties in wind speed and solar irradiance. The MOSGA incorporates crowding distance strategies, fast non-dominated sorting, and an archive selection mechanism to define and preserve the best non-dominated solutions. The total cost, real power loss, and emission were defined as the objectives for the MOOPF-WS problem. In the economic aspect, formulated cost modelling includes both overestimation and underestimation situations related to wind and solar power prediction. Further, uncertainty in load demand is represented by a normal PDF and is considered as a special study case due to its novelty. The effectiveness of MOSGA was validated on the IEEE 30-bus and 57-bus systems considering various combinations of objective functions as well as different loading scenarios. Its performance was comprehensively compared with the other three well-regarded multi-objective optimization algorithms including NSGA-II, MOALO, and MOGOA in terms of Spread metric, Hypervolume metric, and the best compromise solutions for all scenarios. The comparisons showed MOSGA was capable of obtaining well-distributed Pareto fronts and producing better quality solutions compared to the others in all tested scenarios. In addition, MOSGA also obtained better solution quality than significant research in the literature for all the comparable cases. These show the superiority of the MOSGA in dealing with the MOOPF-WS problem.


I. INTRODUCTION
The optimal power flow (OPF) was developed by Carpentier [1], and the OPF problem has become a significant function in the production and operation of modern energy The associate editor coordinating the review of this manuscript and approving it for publication was Akshay Kumar Saha . systems [2]. Generation cost minimization is the primary objective of traditional OPF. However, environmental concerns and technical issues have led to the consideration of different objectives, including emissions, real power loss, voltage stability, and voltage profile. Thus, multi-objective OPF (MOOPF) was formulated to optimize a group of objectives simultaneously in power systems. In general, the MOOPF problem considers fossil fuel-powered thermal generators. Renewable energy sources (RES), which are environmentally friendly and sustainable, provide a better solution, reducing the environmental impact and limiting fossil fuel consumption. Studies of OPF integrating RES have become necessary due to their increasing application in power systems.
In the beginning stage of problem exploration, the MOOPF problem was handled by several conventional methods like gradient-based method [3], linear programming [4], quadratic programming [5], Newton-based method [6], and interior point method [7]. Despite having good convergence characteristics and effectively handling inequality constraints, such conventional methods cannot ensure global optimum result [8]. In addition, conventional approaches are designed specifically for each type of MOOPF problem. In other words, each variant of the MOOPF problem is related to a certain solution approach. Further, these approaches face many difficulties when coping with the attributes of the problem such as integer and discrete decision variables; nondifferentiable objective function; and large range of search space. Moreover, the application of these methods to solving the MOOPF problem is very complex [9].
In the later stage of exploration, with the advent of a new field of computation known as soft computing, a large number of optimization algorithms belonging to the metaheuristic class have been proposed. The main advantages of such methods are that they can deal with large-scale search spaces and are less dependent on problem characteristics. Moreover, these algorithms are capable of estimating multiple points in the search domain simultaneously due to their population-based nature. Therefore, metaheuristic methods are more effective in finding optimal solutions for OPF and MOOPF problems in comparison with conventional methods [9]. In the literature, the MOOPF problem has been studied extensively, and many powerful algorithms have been applied with remarkable results. In an attempt to deal with the MOOPF problem, Bouchekara et al. [10] proposed improved colliding bodies optimization (ICBO). The performance of ICBO was improved through iterations with increasing numbers of colliding bodies. Chaib et al. [11] studied a backtracking search optimization algorithm (BSA) considering several objective functions: emissions, voltage profile, voltage stability, and fuel costs. In [12], the MOOPF problem was formulated as multi-objective problems by assigning weighting factors for power loss, voltage stability, deviation of voltage, and emissions. A differential search algorithm (DSA) was then proposed to deal with the MOOPF problem, and a moth swarm algorithm (MSA) was also applied, with fourteen different cases based on the weighting factor approach [13]. Additionally, Biswas et al. [14] applied differential evolution (DE) to handle the MOOPF problem with numerous case studies considering different objectives. The authors dealt with the operational constraints of the system by using a combination of self-adaptive penalty and superiority of feasible solution as constraint handling methods. A multi-objective differential evolution (MDE) technique was implemented by Shaheen et al. [15] for the MOOPF problem on 118-bus and 57-bus systems. Pulluri et al. [16] presented a mixed crossover operator and self-adaptive strategy to improve the DE for the MOOPF problem. In [17], the authors proposed three different strategies to improve the strength Pareto evolutionary algorithm (ISPEA2). The outcomes showed that these strategies enhanced the uniformity and distribution of Pareto fronts. Recently, the multi-objective dimension-based firefly algorithm (MODFA) was suggested by Chen et al. [18] for the MOOPF problem with nine cases in three test systems. Fuzzy affiliation was also used to extract the best compromise solution. Warid et al. [9] introduced a new method by combining a Jaya algorithm with a quasi-oppositional (QO) strategy. The new algorithm, QOMJaya, was implemented for the MOOPF problem.
In the past decade, research studies mainly focused on handling the OPF and MOOPF problems with thermal generators. Recently, the OPF problem integrated with RES has attracted the attention of researchers due to the increasing penetration of RESs into power systems. Henerica et al. [19] proposed an OPF solution for an off-grid hybrid dieselsolar PV-battery system. In [20], the MBFA technique was applied to handle the OPF framework considering hydrothermal-wind (HTW) systems. Jabr et al. [21] introduced into the OPF a stochastic model of wind power, forecasting power with a probability or relative frequency histogram. To address the same problem, Mishra et al. [22] developed the bacterial foraging algorithm (BFA) for the OPF to consider wind generation. Stochastic wind speed and power distribution were defined using the Weibull distribution. In [23], the authors introduced a wind-generated cost model that included the opportunity cost of excesses and shortfalls in wind power. The proposed method was implemented using the IEEE New England system. In [24], a modified bacteria foraging algorithm (MBFA) was applied to simulate uncertainty in wind energy generation in the OPF framework. In [25], the MBFA was implemented for the OPF problem considering combined thermal-wind generation. Moreover, the voltage security aspect was improved using shunt FACTS (STATCOM) devices to compensate for reactive power. In [26], the author formulated the OPF problem for a hybrid thermal-wind-solar-storage system. A two-point estimating method and the GA method were implemented to handle the proposed OPF strategy. In [27], the authors solved a security-constrained OPF with wind plant by using a fuzzy adaptive artificial physics algorithm. Roy and Jadhav [28] studied the OPF with wind generation by way of an Improved artificial bee colony (ABC) method. In [29], a hybrid algorithm (HA) was applied for the OPF problem with thermal-wind generators; a modified BFA technique and mutation techniques of the genetic algorithm (GA) were combined to develop the HA method. A technique based on the DE was suggested in [30] for the OPF problem considering solar and wind powers. Some recent applications of metaheuristic methods for OPF problem in presence of RES such as constrained multi-objective population extremal optimization (CMOPEO) [31], a combination of fitness-distance balance and adaptive guided differential evolution (FDBAGDE) [32], levy interior search algorithm (LISA) [33], grey wolf optimization (GWO) [34], barnacles mating optimizer (BMO) [35], differential evolutionary particle swarm optimization (DEEPSO) [36], and Levy coyote optimization algorithm (LCOA) [37].
A review of the literature has shown that the research on OPF incorporating RESs is indeed encouraging. Although some attempts have been made to address the OPF with the incorporation of RES into power grids, previous research typically only investigated the problem with the singleobjective function of optimal generation cost [21]- [27] or the two-objective function of generation cost and emission [28]- [37]. In addition to economic and environmental factors, technical aspect also needs to be considered when examining the effect of RES on power grids. For this regard, generation cost, real power loss, and emission should be optimized simultaneously in the MOOPF problem. Further, from various combinations of objective functions under different scenario studies, the conflicts between objective functions can be judged. To our knowledge, solving the MOOPF problem with the stochastic RESs using multi-objective algorithms have not been documented in the literature. In addition, many research [38]- [43] considered only solar power or wind power as a renewable source. Very few research [30]- [36] examined both stochastic wind and solar power when solving the MOOPF problem. Further, most of the previous research considered only the invariable load case, which did not reflect the characteristics of the actual load demand. Further study is therefore needed to investigate the MOOPF problem in a system comprising thermal, wind, and solar power generation with the variable loading condition; this constitutes significant motivation for proposing a new technique to determine trade-off solutions for the MOOPF problem considering stochastic wind and solar power along with uncertain load scenario.
The search group algorithm (SGA) is a robust optimization method proposed in [44]. Notably, SGA has the advantage of striking a good balance between exploitation and exploration, providing powerful searchability for finding the optimum solution. Several recent studies have been done to verify the SGA applicability for various optimization problems such as truss optimal voltage regulation in power systems [45], automatic generation control [46], networked control system [47], steel frames optimization [48], and structure optimization [49]. Moreover, SGA has not been applied to deal with the MOOPF problem with the stochastic RESs.
Inspired by the above motivations, this paper aims to suggest a novel multi-objective SGA (MOSGA) for MOOPF incorporating the stochastic wind and solar power (MOOPF-WS) problem. Three approaches are combined into the original SGA to create the proposed MOSGA. The first one named the fast non-dominated sorting technique is applied for determining different non-dominated fronts in a population.
The second one named crowding-distance calculation is to preserve the diversity of non-dominated solutions in a specific front. The last one named Pareto archive selection is used to store non-dominated solutions in an archive. The archive is updated via a selection mechanism after each iteration to avoid erroneous rejection of potential candidate solutions during the optimization process. Further, a decisionmaking process is implemented to obtain the best compromise solution. The proposed MOSGA is implemented to find the optimal operation points of the thermal, wind, and solar generators with different case studies of objective functions.
This study contributes to the literature as follows: • The MOOPF-WS problem was formulated as a constrained multi-objective optimization problem considering the scheduling of thermal power generation and stochastic wind and solar power generation, whereas the references [38]- [43] considered either solar power or wind power as a renewable source. Three vital objective functions of total operational cost, power loss, and emission were concurrently optimized, whereas only a single objective or two objectives were optimized in the references [21]- [37]. The control variables of the problem included most of the important parameters related to the operation of a power transmission network such as the active power output of thermal generators, the active power output of wind and solar generators, voltage magnitudes at generation buses, settings of transformer tap, and reactive power output of shunt VAR compensators. Most of the previous relevant research in [30]- [36] examined the problem with a fewer number of control variables that did not reflect the actual operating condition of the transmission network.
• A new MOSGA was developed to handle the MOOPF-WS problem. Three improvement techniques are integrated into the search mechanism of the original SGA to convert it into a true multi-objective optimization algorithm. The development of MOSGA aims to provision Pareto optimal front and respective trade-offs for conflicting objectives. Meanwhile, previous significant research [30]- [36] used a weighted factors method for transforming from a multi-objective function into a weighted sum function and applied a single-objective optimization algorithm as the optimizer. This workaround only offered a unique optimal solution after a run, which did not reflect the degree of conflict between objectives as well as not provide multiple options in terms of various operating solutions for decision-makers.
• A practical case study relative to uncertainty in load demand was examined to the MOOPF-WS problem where a normal probability distribution function was adopted to describe the variable situation of load demand. Notably, this case study has been investigated for the first time in the MOOPF-WS problem.
• The effectiveness of MOSGA was tested on adapted 30-bus and 57-bus systems incorporating RES with the VOLUME 10, 2022 consideration of various scenarios. Specifically, in the two first scenarios of each network, the MOOPF-WS problem was investigated in the fixed load situation with two simultaneously optimized objectives, whereas the third scenario was studied in the same load condition for the simultaneous optimization of three objectives. For the last scenario of each network, the problem was solved in the variable load condition while minimizing three objectives concurrently.
• The performance of MOSGA was comprehensively compared with three other multi-objective optimization algorithms including NSGA-II, MOALO, and MOGOA for all scenarios. Comparisons were made in terms of Spread metric, Hypervolume metric as well as the best compromise solution. In all scenarios, statistical results showed the MOSGA achieved better convergence and distribution of Pareto optimal solutions than other multiobjective algorithms. Further, when comparing with previous significant research [30]- [34], MOSGA acquired better solution quality in all the comparable cases. Section 2 introduces mathematical models for the generation costs of thermal, wind, and solar generators. Section 3 formulates the MOOPF-WS problem. Section 4 outlines the proposed MOSGA and its implementation for the MOOPF-WS problem. Section 5 details the simulation outcomes, and Section 6 presents the conclusions.

A. COST OF THERMAL GENERATION
The generation cost of thermal generators can be modeled as follows: where a i , b i , and c i are the cost constants of the i th thermal generator that generates output power P TG,i , and N TG denotes the total number of thermal generators.

B. COST OF WIND GENERATION
Being a renewable energy source, wind power has inherent uncertainty related to variable wind power. Due to insufficiency and unavailability, wind turbines cannot effectively provide schedule power as expected in some cases. As a spinning reserve control unit, the independent system operator (ISO) is in charge of handling the deficit amount. Such a situation is referred to as the overestimation of power from renewable sources and should be included in the generation cost because of keeping the spinning reserve. In the contrast, if available power produced by wind turbines is greater than schedule power, the ISO should pay a penalty cost for the surplus amount. This situation is known as underestimation. Therefore, the cost of wind power generation relates to three components of direct cost, reserve cost, and penalty cost [30].
The direct cost of a wind power plant can be expressed as [30]: C w,j (P ws,j ) = g j P ws,j (2) where g j denotes the direct cost coefficient and P ws,j the scheduled power for the j th wind power plant. As above mentioned, there might be another two costs caused by overestimation and underestimation situations. Precisely, the overestimation situation will lead to the reserve cost, whereas the underestimation one will lead to the penal cost. The wind power plant's reserve cost can be determined as [30]: where P wav,j , K Rw,j , and f w (p w,j ) represent the actual available power, reserve cost constant, and probability density function for the j th wind power plant, respectively.
The wind power plant's penalty cost can be determined as [30]: where P wr,j and K Pw,j denote rated output power and the penalty cost constant for the j th wind power plant, respectively.
The wind speed distribution is defined using the Weibull probability density function (PDF) with two parameters [25], [50]. Hence, the probability of wind speed (v) is expressed as: where k is the shape coefficient and c is the scale coefficient for the PDF. The mean of the Weibull distribution can be stated as: The wind turbine power output can be defined as: where p wr represents the rated output power of the wind turbine; v in , v r , and v out respectively denote the cut-in speed, rated speed, and cut-out speed of wind turbine.
48382 VOLUME 10, 2022 According to Eq. (7), in a couple of zones of wind speeds, the variable wind power is discrete. Specifically, the wind power is zero for wind speed v smaller than v in or for v greater than v out . The wind power is equal to the rated value when v is between v r and v out . Probabilities for these discrete zones can be calculated based on Weibull cumulative distribution function (cdf) as follows [51]: For wind speed being between v in and v r , the wind power is a continuous variable. The continuous zone is related to the probability as follows: Given the probabilities of available wind power at different operating zones, Eqs. (3) and (4) can be expanded to include the PDF of wind power as in Eqs. (11) and (12), as shown at the bottom of the page, respectively.
Remark 1: The process of handling the stochastic wind power can be summarized as follows: Firstly, the wind speed distribution according to the Weibull probability density function (PDF) [41], [52]- [54] is established using Eqs. (5) and (6). With the known parameters of Weibull shape (k) and scale (c), the Monte-Carlo method is adopted for generating the wind speed frequency distribution and Weibull fitting. Next, based on the power output function of a turbine in Eq. (7) and the mathematical equations of Weibull distribution, the available wind power probabilities corresponding to different wind speed zones can be defined by Eqs. (8) - (10). Lastly, thanks to the probability functions for available wind power estimation at different operating zones, we can compute the cost of stochastic wind power using Eqs. (11) and (12).
Overall, the generation cost of wind generators can be defined as: where N WG represents the number of wind generators.

C. COST OF SOLAR GENERATION
Besides wind source, solar PV is another well-known renewable source. Generation of a solar PV plant also includes a summation of direct cost, reserve cost, and penalty cost components. The direct cost of the solar PV plant can be modeled as in Eq. (14) [30]: where h k and P ss,k denote the direct cost constant and scheduled power for the k th solar PV plant. The solar PV plant' reserve cost can be obtained using the below relation [30]: where P sav,k and K Rs,k are the actual available power and the reserve cost constant of the k th solar PV plant, respectively; f s (P sav,k < P ss,k ) denotes the occurrence probability of solar power shortage from the schedule power (P ss,k ); and E(P sav,k < P ss,k ) signifies the expectation of solar power being less than P ss,k .
C Pw,j (P wav,j − P ws,j ) VOLUME 10, 2022 The solar PV plant' penalty cost is as follows [30]: where K Ps,k are the penalty cost constant relative to the k th solar PV plant; f s (P sav,k > P ss,k ) denotes the probability that solar power is more than the schedule power (P ss,k ); and E(P sav,k > P ss,k ) signifies the expectation of solar power being more than P ss,k . Lognormal PDF is used to represent the solar irradiance (G s ) distribution. Hence, the probability of solar irradiance (G s ) is defined by the following equation [55]: where µ is the mean and σ is the standard deviation for the PDF.
The mean of the lognormal PDF can be determined as follows: The solar irradiance (G s ) to energy conversion for solar PV can be determined as follows [56]: where P sr denotes the rated output power of solar PV plant, R c is the certain irradiance point, and G std denotes solar irradiance in the standard environment. The reserve cost in Eq. (15) is rewritten as follows [30]: where P sn− signifies that available power falls short of scheduled power, f sn− denotes the relative frequency of occurrence of P sn− , and N − signifies the number of pairs (P sn− , f sn− ) created for the PDF. Similarly, the penalty cost in Eq. (16) is defined as [30]: (21) where P sn+ signifies that available power is greater than scheduled power, f sn+ denotes the relative frequency of occurrence of P sn+ , and N + signifies the number of pairs (P sn+ , f sn+ ) created for the PDF. Remark 2: Similarly, the process of handling the stochastic solar power can be summarized as follows: The lognormal PDF having characteristics as in Eqs. (17) and (18) is used to express the distribution of solar irradiance. Via Monte Carlo simulation, we can obtain the lognormal fitting and frequency distribution of solar irradiance. Next, we calculate the solar irradiance to energy conversion for solar PV plant using Eq. (19). Given the distribution of actual available solar power, the reserve and penalty costs related to stochastic solar power can be computed using Eqs. (20) and (21), respectively.
Overall, the generation cost of solar generators is as follows: where N SG represents the number of solar PV generators.

III. PROBLEM FORMULATION
The MOOPF-WS problem is set out to optimize predefined objectives while meeting several constraints by defining the optimum values of control variables. Therefore, the MOOPF-WS problem may be expressed accordingly: Minimize: Subject to: where F m (u,x) is the m th objective function; u and x are the vectors of control and state variables, respectively; h(u,x) denotes the set of equality constraints and g(u,x) represents the inequality constraints.

A. PREDEFINED OBJECTIVES 1) TOTAL COST MINIMIZATION
The total cost includes generation costs for the thermal, wind, and solar PV generators, which is defined as a summation of Eqs. (1), (13), and (22):

2) POWER LOSS MINIMIZATION
Real power loss can be expressed by the following equation: where N l denotes the number of transmission lines, G q(ij) denotes the transfer conductance between buses i and j; V i and V j represent the voltage magnitudes at buses i and j, respectively; δ i and δ j represent voltage angles at buses i and j, respectively.

3) EMISSION MINIMIZATION
Producing electricity using conventional fossil fuels would release harmful emissions into the environment. The total emission caused by thermal generators is defined as follows: where  (29) in which N g denotes the total number of generation buses, N t denotes the number of transformers, and N c denotes the number of switchable capacitors.

2) STATE VARIABLES
The set of state variables used is as follows: -P G,1 : real power output at the slack generator.
-Q G : reactive power output of generation buses.
-V L : voltage magnitudes at the load buses. -S L : power flow in the transmission lines. Thus, the vector of state variables (x) is expressed as: in which N d signifies the number of load buses. (32) in which N b denotes the number of buses; P D,i and Q D,i represent real and reactive power demand at bus i, respectively; G ij represents the transfer conductance and B ij the susceptance between bus i and bus j.

IV. MULTI-OBJECTIVE SEARCH GROUP ALGORITHM
The MOSGA is created by integrating fast non-dominated sorting, crowding distance strategies, and Pareto archive selection into the SGA. The original SGA, the new integration strategies, and the proposed MOSGA are detailed as below.
A. SGA SGA is a metaheuristic method developed in [44]. An introduction of the basic SGA is described as follows. To start the optimization algorithm, a population P can be generated randomly as the following equation: where n pop is the population size and n denotes the number of design variables; P ij denotes the j th control variable of the i th individual of P, U [0, 1] denotes a random number ranging from 0 to 1; x j,min and x j,max respectively denote the lower and upper limits of the j th control variable. The objective function value for each individual is computed after population initialization. A number of good individuals are then extracted from population P to create a search group R using a standard tournament selection. Further information on the tournament selection is given in [57]. The search group R is mutated at every iteration to improve global exploration. Depending on the rank in search group R, the worst member is chosen for mutation. The aim is to explore new areas in the search domain. Therefore, new individuals are mutated as follows: in which x j,mut denotes the j th control variable of a specified mutated individual, R :,j signifies the j th column of the search group matrix, E denotes the mean value, t is the coefficient that determines the position of a newly created individual, ε denotes a random variable, and σ denotes the standard deviation.
In SGA, a set comprising each search group member and the respectively generated individuals is a family. Therefore, after the search group is created, each member generates a family by perturbation as follows: where α denotes the perturbation constant, which is reduced by the iteration: in which b means the coefficient, which is defined by a combination of the linear function. SGA comprises global and local phases. The principal aim in the global phase is to explore the entire design domain by creating a new search group from the best member of each  family. To exploit the region of the best current individual, new search group is formed from the best n g individuals from all the families in the local phase. Algorithm 1 in Table 1 describes the SGA's pseudocode.

B. FAST NON-DOMINATED SORTING TECHNIQUE
For each solution p of a set S, two entities are computed [58]: -Domination count n p : the number of solutions that dominate the solution p. -S p : the set of solutions dominated by solution p. The solutions with a zero domination count (n p ) are contained in the first non-dominated rank (P 1 ). Afterwards, for each solution p with zero domination count, the process visits each solution q of set S p and decreases the domination count of each solution q by one. If a member q has a zero domination count, it is located in a different list Q. These members are also placed in the second front (P 2 ). This process continues with each member of Q to determine the third rank (P 3 ) and all non-dominated ranks.

C. CROWDING DISTANCE COMPUTATION
Crowding distance is computed to estimate the density of solutions around a specific solution i in a specific non-dominated front. The population needs to be sorted by the values of the objective function in ascending order. Then, an infinite distance value for each objective function is assigned to the solutions with minimum and maximum objective values, i.e., the boundary solutions. A distance value is assigned to all other intermediate solutions; this value is the absolute normalized difference between the objective values of two neighboring solutions (i+1 and i -1). This computation is employed for other objectives. The value for crowding distance is computed as the total of the individual distance values for each objective. Before calculating the crowding distance, each objective function is normalized. Fig. 1 shows a schematic of the crowding distance computation.
The crowded-comparison operator (≺ n ) is used based on non-dominated level (r) and crowding distance (d) to define the better solution: According to Eq. (48), the solution that belongs to the better non-dominated level is prioritized; if both solutions are in the same level, the one with the better value for crowding distance is prioritized.

D. PARETO ARCHIVE SELECTION
The Pareto archive is utilized with the best non-dominated solutions to maintain a non-dominated set. In the MOSGA, all generated solutions are placed in the advanced archive. Two archives are obtained (current and advanced), and these are combined. The combined archive exceeds the fixed size (n pop ), and so this archive is truncated using the selection strategy suggested by Deb et al. [58] for the next step of the algorithm. Fig. 2 outlines this procedure.
Fast non-dominated sorting is used first to separate the combined archive into non-dominated ranks (P 1 , P 2 , . . . , P n ). The solutions in the first non-dominated front (P 1 ) are chosen first for the new archive. If the first front (P 1 ) is smaller than the archive, the remainder of the archive is selected from other non-dominated fronts in rank order (P 2 , . . . , P n ). The process is implemented until the Pareto archive reaches its limit, assuming that front F k is the final non-dominated front, after which no further sets can be added. Solutions from the last front (P k ) are chosen in descending order of the crowding distance to select the number of solutions for the new Pareto archive accurately.

E. THE PROPOSED MOSGA
The MOSGA begins the procedure by randomly generating an initial population P. Objective values are calculated for all individuals of P. Based on the non-dominated ranks of P, a tournament selection is later employed to choose n g best individuals from P to form an initial search group R for the subsequent two processes (mutation and generation of the families). Afterwards, the advanced archive is yielded, and the current and advanced archives are combined. The selection strategy is implemented to choose n pop best individuals for the new Pareto archive. Finally, a crucial stage of the MOSGA is to select a new search group. This process in both stages is carried out using a tournament selection. The new search group is generated by the best member of each family in the global phase. Meanwhile, the selection mechanism is adjusted in the local phase so that the best n g individuals from the new Pareto archive is extracted to generate the new search group. The MOSGA process executes until the stopping condition is met. The MOSGA's pseudocode (Algorithm 2) is presented in Table 2.

2) OBJECTIVE FUNCTION VALUE
In the MOOPF-WS problem, each objective function value is calculated as follows: (50) in which F m denotes the m th objective for the MOOPF-WS; K p , K q , K v , and K s denote the penalty constants. In this study, penalty constants are set to 10 6 .
The limitation values of the state variables are defined according to Eq. (51): where x lim and x denote the limited values and calculated values of Q G , V L , and S L , respectively. The power flow problem was addressed using the Matpower 6.0 toolbox [59].

3) BEST COMPROMISE SOLUTION
A decision-making method based on a fuzzy membership function is used to define the best compromise solution. The membership function µ ij represents the degree of satisfaction of the i th solution for the j th objective [60]:  where min(F j ) and max(F j ) are the minimum and maximum values of the j th objective, respectively. The normalized membership function of the i th solution can be defined as: where n pf denotes the number of objective functions and n obj denotes the number of non-dominated solutions, respectively. The best compromise solution is the one with a maximum value of the normalized membership function.

4) OVERALL PROCEDURE
The MOSGA implementation for the MOOPF-WS is as follows: Step 1: Define detailed data of the test system and acceptable ranges of control variables; Step 2: Set the controlling parameters for MOSGA, including population size (n pop ), number of search group members (n g ), number of mutations (n mut ), perturbation factor (α), Pareto archive size (N ), maximum number of iterations (Iter max ), and global iteration ratio (GIR); Step 3: Create the initial population P as in Eq. (44); Step 4: Execute a power flow calculation using Matpower 6.0 to evaluate objective function values (F 1 , F 2 , F 3 ) for each individual of P; Step 5: Sort all the individuals in the initial population P using fast non-dominated sorting and crowding distance computations and place them into a Pareto archive; Step 6: Choose the best n g solutions from P to generate the initial search group R k . Set k = 0; Step 7: Begin the main loop, k = k+ 1; Step 8: Perform the mutation phase for n mut individuals using Eq. (45); Step 9: Create families (F i ) for each search group member using Eq. (46); Step 10: Store all newly generated solutions into an advanced archive. Combine the advanced and current archives; Step 11: Choose the best N solutions from the combined archive as the new Pareto archive; Step 12: Select the new search group in two stages: -Global stage: select the best member of each family to create search group R k+1 ; -Local stage: select n g best solutions from the archive to create search group R k+1 .
Step 13: Determine α k+1 according to Eq. (47); Step 14: If k < Iter max , return to Step 7; or else, go to Step 15; Step 15: Solutions achieved: non-dominated solutions in the archive set; Step 16: Define the best compromise solution based on the defined approach; The diversity of non-dominated solutions of a specified algorithm is provided by the spread ( ) indicator, which shows the spread of extreme solutions and the distribution of solutions in the generated non-dominated set. This parameter is defined in the following equation [61]: where is the generated Pareto optimal front, E i is the i th extreme solutions in true Pareto front, and m represents the number of objective functions. The smaller value of the indicator provides better diversity (i.e., the better extent of spread and distribution) in a non-dominated set.

2) HYPERVOLUME
This metric computes the volume related to members of a non-dominated set in the objective space. The hypervolume (HV) can be computed [62]: where v i is a hypercube for each solution i ∈ formed with a reference point W . The HV values after normalization lie within the range [0,1]. It is desirable that the algorithm has a large HV value.

V. SIMULATION RESULTS
The MOSGA was implemented to deal with the MOOPF-WS problem on the adapted IEEE 30-bus and 57-bus systems incorporating RES. Table 3 summarizes the scenarios considered in this study. The initial parameters of MOSGA, including n pop , n g , n mut , α k , and GIR, were chosen as 100, 20, 5, 2, and 0.3, respectively, for all test systems. The optimization process of the MOSGA was independently run thirty times for each scenario. The number of function evaluations (NFEs) was set as 30,000 and 50,000 for the 30-bus and 57-bus systems, respectively. The Pareto archive size was kept at 100 for all trials. Moreover, MOSGA was compared with multi-objective ant lion optimizer (MOALO) [64], multi-objective grasshopper optimization algorithm (MOGOA) [65], and nondominated sorting genetic algorithm II (NSGA-II) [58]. The population size for all algorithms was fixed at 100. In MOGOA, the maximum and minimum values of decreasing coefficient (c max and c min ) were set to 1 and 0.00004,  respectively. NSGA-II used a mutation probability of 0.5, a crossover probability of 0.9, and distribution indexes for mutation and crossover operators of η m = 20 and η c = 20, respectively. Table 4 summarizes the modified IEEE 30-bus system incorporating the RES used in this study; its topology diagram is displayed in Fig. 3. Of note, locations of RESs were chosen similar to the references [30] to make a fair results comparison. In this context, some conventional generators were replaced by respective renewable generators. Table 5 details the fuel cost and emission factors for thermal generators. Tables 6-7 give PDF parameters and cost constants for wind and solar PV generators. After performing 8,000 Monte-Carlo scenarios, Figs. 4-5 depict wind frequency distributions for two wind farms. Similarly, Fig. 6 displays the frequency distribution of solar irradiance for a solar PV plant. Stochastic power output from the solar PV plant is presented in Fig.  7. For the initial scenario, the system has a total cost of 834.6627 $/h, a real power loss of 5.7866 MW, and an emission of 0.2694 ton/h. VOLUME 10, 2022

1) SCENARIO 1: SIMULTANEOUS OPTIMIZATION OF TOTAL COST AND ACTIVE POWER LOSS
The MOSGA was executed to optimize the total cost and real power loss concurrently in scenario 1. Fig. 8 portrays the best Pareto optimal front associated with the best HV metric for scenario 1, which shows a trade-off between the two objectives considered. Table 8 reports the simulation results found by MOSGA, NSGA-II, MOALO, and MOGOA, including optimal values of design variables and objective functions for the best compromise solution. Table 8 also indicates that the MOSGA offered a total cost of 799.1777 $/h and a real power loss of 3.2970 MW for the best compromise solution. Thus, total cost and real power loss reductions for this scenario  Therefore, MOSGA is capable of finding a better optimal solution for this case study.

2) SCENARIO 2: SIMULTANEOUS OPTIMIZATION OF TOTAL COST AND EMISSION
In this scenario, the MOSGA minimized the total cost and emission concurrently. Fig. 9 demonstrates the Pareto front    obtained by MOSGA for this scenario. As stated above, the best Pareto optimal front is the one that corresponds to the best value for the HV metric. The curve in Fig. 9 shows the conflicting nature of total cost and emission due to the equation characteristics of objective functions. Table 8 displays simulation output obtained from the MOSGA, NSGA-II, MOALO, and MOGOA for the best compromise solution. The total cost and emission for the best compromise solution yielded by MOSGA were 784.7058 $/h and 0.2788 ton/h, respectively. These correspond to a reduction in total cost of 5.9853% and an increase in emissions of 3.4754%, as compared to the initial scenario. Compared with NSGA-II,   MOSGA acquired the best compromise solution with lower total cost and higher emission. For a comprehensive comparison, MOSGA not only found a better result than MOALO and MOGOA in both objectives, but it also obtained the best value for the total cost objective among other algorithms for best compromise solution. Hence, MOSGA can provide a good-quality best compromise solution for scenario 2.

3) SCENARIO 3: SIMULTANEOUS OPTIMIZATION OF TOTAL COST, ACTIVE POWER LOSS, AND EMISSION
The total cost, real power loss, and emission were minimized concurrently by way of the MOSGA in this scenario. Fig. 10 portrays the Pareto optimal front found by the MOSGA, indicating the relations between the three objectives. The best compromise solution was also extracted through the decisionmaking method. Table 10 represents simulation outcomes for the best compromise solutions yielded by MOSGA, NSGA-II, MOALO, and MOGOA. Table 10 also shows that the best compromise solution provided a total cost of 807.6874 $/h, a power loss of 3.0442 MW, and emission of 0.1342 ton/h, corresponding to a total cost reduction of  3.2319%, a real power loss reduction of 47.3919%, and an emission reduction of 50.1767%, respectively, in comparison with the initial scenario. It can be seen that MOSGA dominated MOGOA in all three objectives in the best compromise solution. Moreover, MOSGA outperformed NSGA-II and MOALO in two out of three objectives. Therefore, MOSGA had great potential in finding a better compromise solution.
According to the results obtained from this scenario, all three objective functions indicate a remarkable improvement using the MOSGA. Precisely, all the objectives were met simultaneously to the maximum extent, an optimal result for the best compromise solution. Moreover, an essential constraint in the MOOPF-WS problem is that load bus voltage must be maintained within a specific range; Fig. 11 shows voltage profiles of load buses for scenarios 1-3, and it can be seen that the voltages are all within specified boundaries. 48392 VOLUME 10, 2022

4) SCENARIO 4: SIMULTANEOUS OPTIMIZATION OF TOTAL COST, ACTIVE POWER LOSS, AND EMISSION (UNCERTAIN LOAD DEMAND CASE STUDY)
This case study considers a practical variable situation of load demand. Normal PDF is utilized to reflect the uncertain nature of load [66]. Next, we apply the scenario-based method to implement the process of optimization at some discrete load levels. Fig. 12 portrays the load demand uncertainty in relation to the graphic form of normal distribution, in which the horizontal axis deputizes the percentage of network loading. The PDF parameters include a mean (µ l ) and a standard deviation (σ l ) of 70 and 10, respectively.   The color-filled areas manifest four different scenarios of network load demand (P l ) considered. For a certain scenario, the mean loading and its occurrence probability can be estimated by the following equations [66].
Probability of the i th loading scenario is calculated by: in which P lo l,i and P up l,i are the lower and upper bounds of the i th loading scenario. Next, mean of the i th loading VOLUME 10, 2022     Table 11 provides information on the estimated means and probabilities for all loading scenarios. For a specific scenario, the demands in all nodes are adjusted by multiplying them by the percent loading, given the nominal loading of the network of 100%. At each scenario, the proposed MOSGA algorithm optimizes the multi-objective function as formulated in Eq. (50). Hence, the control variables of scheduled power from all generators, voltage magnitudes at generation nodes, transformer tap settings, and VAR compensators' output are optimized in each scenario of loading. Of note, the case study with the implementation of OPF at certain time intervals for various objectives optimization is similar to the practical operation condition of an electrical network.
Regarding the boundaries of control variables, they were set the same as in the fixed loading case. Simulation results related to various scenarios of loading of the 30-bus network are tabulated in Table 12 using MOSGA optimizer. Figs. 13-16 depict the Pareto optimal fronts found by the   MOSGA for scenarios 4.1-4.4, respectively. As for each loading scenario, the optimized values of the total cost (TC sce ), active power loss (PL sce ), and emission (E sce ) are also reported in Table 12. Afterwards, the expected values of the total cost (ETC), active power loss (EPL), and emission (EE) over all scenarios can be computed according to Eqs. (60)-(62) [67]: (62) in which, n sce and δ sce are the total numbers of load scenarios and the likelihood of a scenario, respectively (Table 11).
From the simulated outcomes in Tables 10 and 12, the total operational cost, power loss, and emission achieved by the uncertain loading case were justifiably lower than those achieved by the fixed loading case. The change in network loading resulted in a significant change in schedule power of generators, especially in solar PV and wind  power generators. Specifically, the optimal schedule power of generators increased when the load demand of network increased. To satisfy the increasing load demand, solar PV and wind generators were preferred to be scheduled. This may contribute to the notable improvement level of network performance profiles in examined discrete loading scenarios as well as of expected performance profiles over all scenarios. Compared with the fixed loading case, the uncertain loading case proposed a more flexible and effective scheduling of generators that led to lower total cost and power loss, and emission. Besides MOSGA, three multi-objective optimization algorithms including NSGA-II, MOALO, and MOGOA were also applied for results comparison in this case study. A comparison of results optimized by the four methods is summarized in Table 13. Clearly, MOSGA dominated MOALO and MOGOA in two out of three objectives in the best compromise solutions for all considered discrete loading scenarios. These helped MOSGA achieve lower three expected values (i.e., total cost, power loss and emission) over all scenarios than those of MOALO, whereas MOSGA achieved lower expected values of total cost and emission over all scenarios than those of MOGOA. Compared with NSGA-II, MOSGA obtained the best compromise solution with the better two out of three objectives for the second and third discrete loading scenarios while it excelled at only one objective for the first and fourth loading scenarios. Notably, VOLUME 10, 2022   when considering the overall performance across all loading scenarios, MOSGA was superior in two out of three expected performance profiles over NSGA-II. Further, the voltage profiles at buses of the 30-bus network for all scenarios after the optimization by MOSGA are plotted in Fig. 17. The solutions obtained by MOSGA led to a significant improvement in bus voltage profiles for all loading scenarios while ensuring the permissible limits of bus voltage. Table 14 describes the 57-bus system incorporating RES; its topology diagram is shown in Fig. 18. Notably, locations of RESs were selected the same as the references [36] by replacing some conventional generators with respective renewable generators. Table 15 details the fuel cost and emission factors  of thermal generators. Tables 16-17 provide the PDF parameters and cost constants of wind and solar PV generators. Figs. 19-20 depict wind frequency distributions for two wind farms after Monte-Carlo simulations with 8000 scenarios obtained for each. Fig. 21 displays the frequency distribution of solar irradiance for a solar PV plant. Stochastic power output from the solar PV plant is depicted in Fig. 22. For the initial scenario, the system has a total cost of 51795.8561 $/h, a real power loss of 27.5627 MW, and an emission of 2.6474 ton/h.

1) SCENARIO 5: SIMULTANEOUS OPTIMIZATION OF TOTAL COST AND ACTIVE POWER LOSS
In this scenario, the MOSGA method simultaneously optimized total cost and real power loss objectives. Fig. 23 presents the Pareto optimal front obtained by the MOSGA for scenario 5. Similar to scenario 1, scenario 5 also shows the conflicting relationship between the total

2) SCENARIO 6: SIMULTANEOUS OPTIMIZATION OF TOTAL COST AND EMISSION
The total cost and emission were simultaneously optimized using the MOSGA in this scenario. The Pareto optimal front yielded by MOSGA is described in Fig. 24. Table 18 shows the simulation outputs from four algorithms for scenario 6. The total cost and emission were reduced to 30575.27 $/h and 1.0043 ton/h for the best compromise solution, corresponding  to 40.9697% and 62.0652% reductions. Therefore, both objectives were significantly improved upon compared to the initial scenario. In this case, MOSGA obtained better results than MOALO and MOGOA for the total cost and emission objectives. Meanwhile, the optimal solution of MOSGA was better than NSGA-II in the objective of total cost only. Hence, MOSGA had great potential in finding better compromise solution.

3) SCENARIO 7: SIMULTANEOUS OPTIMIZATION OF TOTAL COST, ACTIVE POWER LOSS, AND EMISSION
In scenario 7, the MOSGA simultaneously minimized the total cost, real power loss, and emission. Fig. 25 portrays the Pareto optimal front found by MOSGA. Table 19 gives the simulation results generated by the MOSGA, NSGA-II, MOALO, and MOGOA for the best compromise solution in this scenario. For MOSGA, the total cost, real power loss, and emission were reduced to 30663.58 $/h, 11.8439 MW, and 1.0846 ton/h, respectively. Hence, MOSGA found a total cost reduction of 40.7992% and a power loss reduction of 57.0292%, along with a reduction of 59.0327% in emission for the best compromise solution. Although MOSGA  was not able to directly dominate NSGA-II, MOALO, and MOGOA, MOSGA obtained better total cost and real power loss values than the others in the best compromise solution. Hence, MOSGA was more advantageous in handling the MOOPF-WS problem with the large-scale system.
It can be seen that the implementation of the MOSGA for the MOOPF-WS problem offers a notable improvement for all three objective function values. Fig. 26 depicts the voltage profiles of load buses for scenarios 5-7, which shows that the voltage limits were all within acceptable limitations. Similarly, experiment outcomes for the uncertain load case of the 57-bus network obtained by MOSGA are presented in Table 20. Figs. 27-30 show the Pareto optimal fronts found by the MOSGA for scenarios 8.1-8.4, respectively. Comparing the results between Table 19 and Table 20 shows the uncertain case of load demand offered justifiably lower network performance profiles of total cost, power loss, and emission than those of the fixed loading case. Again, the scheduling of generators from the uncertain load case was more suitable with the variable load demand. Solar PV and wind power generators were preferable to schedule due to their positive impacts in reducing network operation cost,   power loss, and emission. To confirm the effectiveness of MOSGA, the obtained best compromise solutions related to four scenarios of loading were compared with those of NSGA-II, MOALO, and MOGOA (Table 21). Obviously, MOSGA showed its domination in two of three objectives when compared with NSGA-II, MOALO, and MOGOA for all discrete loading scenarios. As a result, MOSGA acquired lower three expected values of the total cost, power loss, and emission over all scenarios than MOALO. In comparison with NSGA-II and MOGOA, MOSGA obtained lower two expected values of total cost and emission over all scenarios. Moreover, for all scenarios, bus voltage profiles of the 57-bus network were significantly enhanced and completely satisfied with the limits after the optimization by MOSGA (Fig. 31).

C. STATISTICAL COMPARISONS AND ANALYSIS
This section demonstrates the comparisons among MOSGA, NSGA-II, MOALO, and MOGOA based on and HV metrics. The initial parameters of all algorithms were retained at the same values as in the preceding section. Table 22 shows statistical results attained by the MOSGA, NSGA-II, MOALO, and MOGOA for the indicator for scenarios 1-8. Boxplots of the statistical analysis for metric are depicted in Fig. 32. As can be seen, the MOSGA provided the best solutions for the metric, with the best distribution for Pareto optimal fronts among the four methods for all scenarios. Hence, the proposed MOSGA showed superiority over the others for scenarios 1-8.

1) SPREAD ( ) METRIC
2) HYPERVOLUME (HV) METRIC Table 23 gives the optimization results of four methods for the HV metric for scenarios [1][2][3][4][5][6][7][8] shows the boxplots of the HV metric for scenarios 1-8. From Table 23 and Fig. 33, MOSGA was better than the other methods for the average value of the HV metric, which indicated that solutions yielded by the MOSGA had the best convergence and diversity. Therefore, MOSGA found a better quality Pareto optimal front in comparison to NSGA-II, MOALO, and MOGOA. As seen in Figs. 33, MOSGA was also superior to the other methods in terms of stability and robustness.

3) COMPUTATIONAL TIME
In this study, MOSGA, NSGA-II, MOALO and MOGOA were executed on a 2.6 GHz 4-core computer with 16 GB of RAM. Table 24 presents the average computational times of MOSGA, NSGA-II, MOALO and MOGOA for all scenarios. The computational times of MOSGA were significantly faster than those of NSGA-II, MOALO, and MOGOA for all case studies. Hence, MOSGA outperformed the other methods for the solution quality, robustness, and computational times in all cases.

VI. CONCLUSION
This study suggested a new MOSGA for the MOOPF-WS problem, in which thermal, wind and solar PV generators were integrated into the grid. Different PDFs were used to model uncertainties relative to RES and load demand as well. The MOSGA integrated fast non-dominated sorting and crowding distance methods in defining the non-dominated ranks and densities of the solutions obtained. Furthermore, the Pareto archive selection strategy was applied for the distribution maintenance of the non-dominated solutions. The MOSGA was implemented successfully to the adapted 30-bus and 57-bus systems incorporating RES where various combinations of objective functions and different loading scenarios were examined. The findings showed that system performance is significantly improved for the best compromise solution in both fixed and variable loading scenarios. Notably, in scenario 7, the total cost, real power loss, and emission were decreased by 40.7992%, 57.029%, and 59.0327%, respectively, in comparison with the initial scenario. For the variable loading case, the expected performance profiles over all discrete loading scenarios were significantly lower than the respective ones found by the fixed loading case. This proved that solving OPF in the variable load situation offered more flexible and effective scheduling of generators than in the fixed load situation, resulting in lower expected performance profiles. Moreover, the performance of the MOSGA was compared with NSGA-II, MOALO, and MOGOA using and HV indicators. In all cases, the MOSGA showed superiority in convergence and diversity of Pareto optimum solutions compared with the other three methods. As a result, the MOSGA could detect a broader range of non-dominated solutions for all objective functions. Furthermore, when comparing with significant research in the literature, MOSGA obtained better solution quality in all the comparable cases. These revealed the strengths of the MOSGA and validated its potential in dealing with the MOOPF-WS problem. In future studies, the MOOPF-WS problem can be expanded in some specific aspects such as connection consideration of smallhydro source, electric vehicles, and flexible AC transmission system (FACTS) devices.