A modified white shark optimizer for optimal power flow considering uncertainty of renewable energy sources

This paper presents a novel approach to solve the optimal power flow (OPF) problem by utilizing a modified white shark optimization (MWSO) algorithm. The MWSO algorithm incorporates the Gaussian barebones (GB) and quasi-oppositional-based learning (QOBL) strategies to improve the convergence rate and accuracy of the original WSO algorithm. To address the uncertainty associated with renewable energy sources, the IEEE 30 bus system, which consists of 30 buses, 6 thermal generators, and 41 branches, is modified by replacing three thermal generators with two wind generators and one solar PV generator. And the IEEE 57-bus system, which consists of 57 buses, 7 thermal generators, and 80 branches, is also modified by the same concept. The variability of wind and solar generation is described using the Weibull and lognormal distributions, and its impact on the OPF problem is considered by incorporating reserve and penalty costs for overestimation and underestimation of power output. The paper also takes into account the unpredictability of power consumption (load demand) by analyzing its influence using standard probability density functions (PDF). Furthermore, practical conditions related to the thermal generators, such as ramp rate limits are examined. The MWSO algorithm is evaluated and analyzed using 23 standard benchmark functions, and a comparative study is conducted against six well-known techniques using various statistical parameters. The results and statistical analysis demonstrate the superiority and effectiveness of the MWSO algorithm compared to the original WSO algorithm for addressing the OPF problem in the presence of generation and demand uncertainties.

• Introducing a modified white shark optimization algorithm (MWSO) by incorporating Gaussian barebones and quasi-oppositional-based learning to enhance exploration capabilities and improve convergence rates compared to the original WSO.• Validating the effectiveness of the MWSO algorithm by applying it to 23 benchmark functions and comparing its performance against efficient competitors using various statistical metrics.• Modifying the IEEE 30-bus to include wind and solar power plants, and utilizing both the MWSO and original WSO algorithms to address the optimal power flow (OPF) problem through four different objective functions.• Conducting practical scenarios that consider the uncertainty of generation and demand, as well as ramp rate limits of thermal power plants, and analyzing the results obtained from the proposed MWSO algorithm and the original WSO algorithm in these simulation scenarios.• Using a modified IEEE 57-bus system to demonstrate the scalability of the proposed MWSO.
The obtained results clearly demonstrate the superiority and dominance of the developed MWSO algorithm over the traditional WSO algorithm in effectively addressing the OPF problem.
The outstanding portions of the present study are: "Different cost models" section outlines the different cost models that include thermal, wind, and solar power costs."Objective functions and system constraints" section presents the various OPF objective functions and corresponding constraints.Then, "Proposed MWSO methodology" section presents the modified algorithm (MWSO).The simulation results, comprising real-world case studies using the MWSO and WSO methods, are given in "Simulation results and discussion" section, in addition to the statistical analysis using the Wilcoxon signed rank test.Also, this section includes the experimental results

Cost of thermal power
The produced thermal power charges a cost that can be calculated using (1), where the valve point impact of thermal plants has been taken into consideration while calculating the cost of thermal power to provide more accurate values.
where P Thi is the output power of the i-th thermal plant, while a Thi , b Thi , and c Thi indicate the cost coefficients of the i-th thermal plant.N Th indicates the number of thermal plants, while d Thi and e Thi indicate the coefficients of valve point loading, and P min Thi denotes the minimum amount of power produced from the i-th thermal plant.The values of all mentioned coefficients in this equation are listed in Supplementary Material Table 2A.

Components of wind power cost
In contrast to thermal power, wind power is subject to considerable uncertainty.Accordingly, the cost of production using wind is computed differently, as stated below.

Direct component
The power that is intended to be generated by wind turbines has a direct cost that can be estimated as follows: where P schwj represents the intended wind power of the j-th wind plant and dw j indicates the coefficient of its direct charge.

Uncertain components
Given the variable character of wind power, two scenarios are possible.The first of these scenarios comes about if the actual production of wind turbines is less than what was anticipated to be produced.This is known as overestimation, and a commitment to the spinning reserve must be made to compensate for it.According to that, a reserve cost is required, which is computed as follows: where K reswj corresponds to the reserve cost coefficient for the j-th wind power plant and P available wj signifies the actual available power from the same plant.The PDF of the wind power from the j-th wind plant is signed as f wind j .
In the second scenario, the amount of electrical power actually provided by the wind turbines may be greater than what was anticipated.If it is not possible to use the extra electrical power, traditional generators' output must be reduced.A penalty fee equal to the excessive power is due from ISO.The definition of the penalty cost corresponding to a wind plant can be clarified by (4): where, K penwj signifies the penalty cost coefficient, and P rated wj states to the rated power of a wind plant (j).

Probabilistic power of wind plants
In this part, the probabilistic power of wind plants, the term " f wind j (P wind j ) " in (3) and (4), will be determined.The Weibull probability density function (PDF) works well with the wind speed distribution 32,53 .Following the Weibull PDF, the following formula is utilized for calculating the probability of wind speed (Wind v ): (1) C Th (P Th ) = N Th i=1 a Thi + b Thi P Thi + c Thi P 2  Thi + d Thi × sin(e Thi × (P min Thi − P Thi )) , C reserve wj = K reswj P schwj − P available wj = K reswj P schwj 0 P schwj − P wind j f wind j P wind j dP wind j , (4) C penalty wj = K pen wj (P available wj − P schwj ) = K pen wj P rated wj P schwj (P wind j − P schwj )f wind j (P wind j )dP wind j , where the letters k and c , respectively, stand for scale and form factors. Weibull distribution's mean is calculated as follows: The gamma function, which is represented by the sign Ŵ in (6), is provided by: After conducting 8000 Monte-Carlo simulation scenarios, Figs. 1 and 2 reveal the frequency distribution of the wind based on Weibull fitting for the wind plant at bus 5 and the wind plant at bus 11, respectively.The applied values of the Weibull distribution have been listed in Supplementary Materials Table 1A.
The wind speed influences the wind plant's output power of.According to Ref. 33 , the following is the formula for wind turbine power output: In this formula, the cut-in speed is shown by Wd vin , the cut-out speed is indicated by Wd vout , and the rated wind speed is indicated by Wd vr .The wind turbine's rated power is shown by the variable P ratedw .
It is possible to establish the probability of output power from wind plant in the discrete region as follows 54 : Regarding the continuous zone, the following formula can be used to determine the probabilities for the power that the wind plant will produce: (5)

Components of solar power cost
The total cost of producing electricity from solar system may be broken down into a direct cost and an uncertainty cost, much like the cost of wind power.

Direct component
The direct component of solar power cost is estimated using (12).
where P schsk represents the intended solar power of the k-th wind plant and ds k indicates the coefficient of its direct cost.

Uncertain components
Similar to how wind energy is estimated, the cost of producing power from solar plants is determined in both overestimation and underestimation scenarios.Consequently, the reserve cost of solar power in case of overestimating is determined by the following formula: As, K ressk indicates the reserve charge coefficient for the solar plant (k) , and P available sk signifies the available output of the same solar plant, while f solark (P available sk < P schsk ) signifies the deficiency existence probability in the production of solar plant, and the expectation of being the output of solar plant below the P schsk is denoted by E(P available sk < P schsk ) .And the penalty cost of solar power in case of underestimating is determined by: where, K pensk signifies the penalty cost coefficient, f solark (P available sk > P schsk ) expresses the probability of the unused solar power produced from the solar plant (k) , while E(P available sk < P schsk ) signifies the expected remain- ing output power from the solar plant (k).

Probabilistic power of solar plant
The variable solar irradiance (I) impacts the output power of the solar plant.Equation ( 15) provides a clarification on how the probability of sun irradiance is calculated based on the lognormal PDF 55 .(11) C penaltysk = K pensk (P available sk − P schsk ) = K pensk × f solark (P available sk > P schsk ) × [(E(P available sk < P schsk ) − P schsk ], As, the irradiance probability is denoted by f I (I) , the mean of the lognormal distribution is denoted by µ , while the standard deviation is signified by σ , respectively.While the mean of lognormal, M lgn is calculated by: In this regard and after performing 8000 Monte-Carlo scenarios, the lognormal distribution for the solar plant is shown in Fig. 3.The applied values of the lognormal distribution are listed in Supplementary Material Table 1A.
Consequently, the sun irradiation vs. the energy conversion of the solar plant can be presented as follows 56 : where in this formula, I std signifies the solar irradiance when the environment is a standard i.e. (800 W/m 2 ), the symbol R c specifies a specific value of irradiance (120 W/m 2 ), and P solarr refers to the rated power of the solar PV system.The reserve charge of the solar power that are stated in ( 13) can be rewritten after determining the probabilities of solar power as follows: where P sn− signifies the unavailability of solar power (lesser than the schedule power) as indicated by the left half plane of the schedule power of the solar plant ( P schsk ) inside Fig. 4. f sln− signifies the relative frequencies of the P sn− occurrence.N − signifies the number of discrete bins on the left plane of the schedule power of the solar plant.
While, the penalty cost that are previously stated in ( 14) can be rewritten as follows: where, P sn+ signifies the surplus of solar power (higher than the schedule power) as indicated by the right half plane of the schedule power of the solar plant ( P schsk ) provided in Fig. 4. f sln+ gives the relative frequencies of the P sn+ occurrence.N + signifies the number of discrete bins on the right plane of the schedule power of the solar plant.
( The reserve charge coefficient ( K res ) for the two wind plants and the solar plant are constant value equals to 3 for all of them, and the penalty charge coefficient ( K pen ) equals to 1.5 for all of them.The direct charge coef- ficient of the two wind plants (dw j ) and the solar plants (ds k ) equals to 1.75, while the direct charge coefficient for the solar plant equals to 1.6.Unless otherwise specified, these values will be utilised in the case studies in "Conclusion" section.

Cost of carbon emissions
The thermal plants are sources for carbon emissions ( CAE) , thus Eq. ( 20) provides an estimation for these emissions.
Here, ϕ Th,i , Thi , ω Th,i , τ Th,i , and ζ Th,i signify emissions coefficients of the thermal plants.The values of these coefficients are given in Supplementary Materials Tables 2A and 3A for the two systems.These emissions in tonnes are translated into cost, C CE in $/h as per Eq. ( 21), where, C Taxc is the tax of emissions in $/tonne.

Objective functions
Minimizing the total cost of production with and without enforcing a tax on carbon emissions, minimizing carbon emissions, and minimizing power losses are the objective functions of this optimal power flow model.The different objective functions of this model can be expressed as follows.
Minimizing the total cost of production without enforcing tax on carbon emissions (F 1 ) F 1 Can be formulated by combining all the above-mentioned different costs in "Different cost models" section.Therefore, F 1 can be written as: Minimizing the total cost of production with enforcing tax on carbon emissions (F 2 ) F 2 Can be formulated by adding all the costs in F 1 into the emissions cost that was expressed in (21).Therefore, F 2 can be written as: ϕ Th,i + � Th,i P Thi + ω Th,i P 2 Thi + τ Th,i e (ζ Th,i P Thi ) .
( F 3 Is to minimize the total carbon emissions of the thermal plants in (20).Therefore, it can be written as: Minimizing the total power losses (F 4 ) The power losses of the power system can be formulated as follows: where δ ij expresses the difference in voltage angles between buses i and j ; NTL signifies the number of transmis- sion lines; and G ij expresses the transfer conductance.Consequently, F 4 can be formulated as follows:

Problem constraints
The problem of OPF is constrained by some conditions that must not be violated.Some of these constraints are equality and the others are inequality constraints.

Equality constraints
The equality constraints are dedicated to ensuring that the generated (active and reactive) power in the system equals to the consumed (actives and reactive) power in addition to the power loss: where NB stands for the network buses number.B ij and G ij stand for the susceptance and conductance among bus i and bus j , respectively.δ ij is the voltage angle of bus i minus the voltage angle of bus j .The real components of the produced and consumed power at bus i are expressed by P Gi and P Di , respectively, while the reactive components of the consumed and generated power are expressed by Q Di and Q Gi , respectively.

Inequality constraints
These constraints define upper and lower limits for the operation of system components such as the generators, transmission lines, and load buses.
Generator limits.Lower and higher limits govern the functioning of all power generators in the network.There are limits for the active power, reactive power, and voltage of the generator as showed by ( 29), (30) ϕ Th,i + � Th,i P Thi + ω Th,i P 2 Thi + τ Th,i e (ζ Th,i P Thi ) . ( Vol where N T and N C refer to the number of transformers and shunt compensators in the network, respectively.
Limits of ramp rate.The ramp-rate restrictions of thermal generators can be identified in the following manner: where,P 0 Thi signifies the previous hour's output power from the thermal plant (i) .Ur i andDr i signify the up and down limits of ramp-rate for the i-th thermal plant, respectively.
Limits of load buses.The voltages of load buses are constrained by lower and upper limits that can be clarified as follows: where N LB denotes the load buses number in the grid.There is another important parameter related to the load buses, which is the voltage deviation of load buses which can be calculated as follows: It is a measure of the quality system's voltage.This indicator is defined as the total deviation of all load buses buses in the system from the nominal value of 1 p.u.
Capacity of transmission lines.The capacity of the transmission lines in the network should not exceed a certain limit.This constraint can be clarified as per Eq.(38), where N L denotes the of transmission lines number in the grid.

Original WSO algorithm: an overview
This section clarifies a short description of the mathematical formulation of the original WSO 40 , which was designed to depict how white sharks behave when foraging.This involves pursuing and chasing prey.The great white shark is capable of locating prey (i.e., a food source) at the depths of the ocean.The location of the food supply in a specific search area is unknown, though.In this situation, white sharks must conduct extended searches to find food sources in the ocean's depths.The three activities of the great white sharks to identify prey (i.e., the best food source) are as follows: (1) moving towards prey is based on the waves' hesitancy as a result of their movement.In this situation, the white shark utilizes a wavy movement to locate prey using its related sense of hearing and smell.As well as in (2) its haphazard quest for prey in the ocean's depths.Great white sharks do this by swimming in the direction of their prey and remaining close to the best one, and (3) white shark behavior in seeking the nearest prey.In doing so, the great white shark mimics fish school behavior by swimming toward the best white shark that is close to the optimum prey.Based on such actions, all white shark sites are adjusted around the global possible solutions if the prey is not discovered properly.The mathematical expressions for these activities are as follows.

Initialization
The WSO algorithm is a population-based algorithm like several metaheuristic optimization techniques.The candidate solutions to an optimization problem with an n population size (i.e., white shark) and d dimensional space are expressed as per Eq. ( 39). (32) Thi − P Thi ≤ Dr i , if power generation reduces, where w symbolizes the position of all sharks in the searching space, d indicates the number of a chosen vari- ables for a specific problem.

Movement towards prey
White sharks spend the majority of their time seeking and chasing prey because they are creatures with a strong need to survive.They often employ all available tactics to follow and track prey while utilizing their exceptional hearing, sight, and smell skills.A white shark moves to its prey in an undulating motion that may be expressed by Eq. ( 40) when it locates its prey based on the hesitation of the waves it hears during the movement of the prey.
where j = 1,2, . . ., n , represents the white shark's index for n search agents; the new velocity of the jth white shark in (k − 1)th step is denoted by v j k−1 ; w gbest,k−1 is the global best solution obtained so far in the iteration (k − 1)th ; best symbolizes the jth best known position for the swarm and v j k−1 shows the optimal position of white sharks' index vector, expressed using Eq. ( 41); c 1 and c 2 are random number between [0,1]; White shark forces govern- ing the influence of w gbest,k−1 and w v j k−1 best on w j k−1 are represented by p 1 and p 2 , respectively, which are computed by Eqs. ( 42) and ( 43); eventually, µ is the constriction factor to represent the convergence characteristics of the white sharks, formulated as per Eq. ( 44).
Here, a and A are the current and maximum iterations, respectively; p lb and p ub are respectively the mini- mum and maximum velocity to accomplish good movement for the white sharks, and their values are 0.5 and 1.5, respectively.
where, τ represents the acceleration factor whose value is 4.125.

Movement towards optimal prey
Excellent white sharks spend the majority of their time looking for prospective prey, whether the location of the prey is ideal or not.Accordingly, the locations of white sharks are continually varying.When they either hear the waves produced by the motion of prey or detect its scent, they usually proceed toward the prey.In this situation, Eq. ( 45) illustrates the motion of white sharks as they proceed toward the prey.
where w j k represents the new position of the white sharks in the in (k − 1) th iteration; x and y symbolize one dimensional binary vectors represented by Eqs. ( 46) and ( 47); w o is a logical vector given by Eq. ( 48); → is a negation operator; the lower and upper limits of the search area are represented by l and u , respectively; f is the frequency of the wavy movement of the white sharks, which is expressed as per Eq. ( 49); rand is a random number within the interval [0,1]; finally, the white shark's movement force, denoted by mv, grows with iteration, as the shark gets closer to its prey, as expressed in Eq. ( 50).
where ⊕ is a bit-wise XOR operation.(40) Vol.:(0123456789) Here, the maximum and minimum frequency of the undulating movement of the white sharks are represented by f max and f min , respectively.Whose values are 0.75 and 0.007, respectively.
where d o and d 1 are respectively two constant positive numbers that utilized to manage the behaviour of explora- tion and exploitation phases.mv demonstrates how the white shark's acute sense of hearing and smell improves with repetition.

Movement towards the best white shark
Great white sharks have the ability to sustain their proximity to their best ones that is closer to prey.Equation ( 51) provides a formalization of this phenomenon.where ' w j k is the updated position of the jth white sharks with respect to the location of the prey; sgn(r 2 − 0.5) provides either 1 or − 1 to reverse the direction of the search; r 1 , r 2 , and r 3 are normally distributed numbers between [0,1]; the distance between the prey and white shark is illustrated by − → D w , as represented by Eq. ( 52); s demonstrates the efficacy of olfaction and vision in white sharks as they track other white sharks that are in proximity to ideal food, which is formulated as given in Eq. ( 53).
where rand is a random number within the range [0,1]; w j k−1 represents the current position of the white shark in respect to the best position, w gbest,k−1 .
where g is a positive constant number, whose value is 0.0005 to manage the behaviour of exploration and exploitation stages.

Fish school behavior
The first two ideal candidates were kept, and the location of other white sharks was modified in accordance with these optimal positions to mathematically imitate the behavior of the school of white sharks.The following formula is presented to characterize white shark schooling behavior: Equation (54) shows that white sharks may adjust their location such that it matches that of the best white shark, which has now moved into an ideal location, extremely near to its meal.Great white sharks (i.e., search agents) would end up in a location in the search space that is almost ideal for their prey.Collective actions, such as schooling fish behavior or white sharks migrating to find the best white shark, are indicative of WSO and expand the potential for more fruitful exploration and exploitation.The complete flowchart for the original WSO algorithm is represented in Fig. 5.

Modified WSO algorithm
In the original form of the WSO, the great white sharks move toward their prey spot using a single approach, which may cause the algorithm to miss additional favorable positions in the surrounding area.In other words, the white shark optimization (WSO) algorithm tends to lose diversity as its evaluation progresses, which can cause challenges with convergence speed and accuracy.Additionally, the WSO has been applied for solving some of the complex optimization challenges as reported in "Introduction" section; however, it has some limitations, such as unbalanced exploration and exploitation, a propensity to become stuck in suboptimal local areas, and a sluggish convergence speed.Therefore, in this study, a new version of the WSO is introduced to overcome these issues, which incorporates the Gaussian barebones (GB) and opposition-based learning (QOBL) strategies.The developed MWSO algorithm is exhibited as follows: Initially, the population is randomly generated, and the optimal individual in the current individual is identified based on the objective function.Then, update the position of the white sharks using Eqs.(51) and (54).The GB is adopted to find superior positions for the updated white sharks.Eventually, the QOBL mechanism is employed to update the individuals again.The increasing diversity of the population improves the exploitation ability and enhances the convergence speed and accuracy of the MWSO algorithm.The strategies of the GB and QOBL are illustrated in the following subsections.Table 1 exhibits the complete pseudo code of the proposed MWSO.Eventually, Fig. 6 shows the complete flowchart of the proposed MWSO algorithm. '

Gaussian barebones strategy
As previously stated, during the later phase of evaluating the WSO algorithm, the variety of the white shark optimization (WSO) algorithm diminishes, leading to potential issues with convergence speed and accuracy.The GB strategy facilitates the selection of the most suitable direction for white sharks and enables them to steadily progress towards the optimal solution, hence preventing premature convergence to local optima.Consequently, once the positions of all search agents have been updated, the inclusion of GB's randomization features into WSO is employed to augment the variety of the population.This maintains a balance between the algorithm's local exploitation and its capability for global search, resulting in enhanced convergence speed.The GB strategy is derived from the bare-bones PSO (BBPSO) algorithm 57 .In this strategy, the parameter R is utilized to guide each individual.If the chance of random generation is less than R, the individual's location is updated using the Gaussian distribution in the next assessment; alternatively, the concept of differential evolution is included.Eventually, the GB strategy may be expressed as follows: where w j k,GB represents the new position for the jth white shark using the GB mechanism; w gbest,k−1 symbolizes the global best solution obtained so far in the iteration (k − 1)th ; r 4 is random number within the interval [0,1]; G denotes the Gaussian distribution with mean w gbest,k−1 + w j k /2 and standard deviation w gbest,k−1 − w j k ; j1 , j2 , and j3 are three randomly chosen individuals that are diverse from the current individual, j ; w j k is the updated position of jth white shark using fish school behaviour in the current iteration kth.

Quasi-oppositional movement strategy
The opposition-based learning (OBL) method, initially proposed by Tizhoosh 58 , can speed up convergence and boost solution quality by simultaneously considering both the current solution and the exact opposite one.According to probability theory, each answer has a 50% chance of being better than the other.The best of the two inverse solutions is picked as the candidate solution to improve the evolutionary algorithms' search efficiency.The OBL approach has been successfully used across a wide range of challenges.Definitions of opposing numbers and opposite points are provided in this work 59 to help clarify the notion of opposition-based learning.(55)   Nevertheless, it's worth noting that OBL has certain development processes and that several researches have demonstrated that quasi-opposition-based learning (QOBL) is more successful than OBL 60 .So, the definitions of the quasi-opposite number and point are as below: Quasi-opposite number.A random number z in the search region [a 1 , a 2 ] has a quasi-opposite number z qo , which may be written as.( 56) Table 1.The pseudo code of the MWSO algorithm.determine the distance between the prey and white shark using Eq. ( 52).

22:
if == 1; then 23: updated position of the white sharks with respect to the location of the prey using Eq. ( 51).

24:
Else 25: updated position of the white sharks with respect to the location of the prey using Eq. ( 54). ) in a searching region with a d-dimensional problem is calculated as per Eq. ( 59).
The QOBL methodology may be utilized not only during the initialization phase but also throughout the evolutionary phase of a WSO algorithm for updating the individuals.In the current study, the solution obtained by the Gaussian barebones process utilizing Eq. ( 55) has the potential to be substituted with a quasi-opposite solution.

Simulation results and discussion
In this section, the performance of the MWSO algorithm is evaluated and quantified using the 23rd standard benchmark functions.These benchmark functions comprise seven unimodal functions, six high-dimensional multi-modal functions, and eight fixed high-dimensional multi-modal functions.Furthermore, the application of the MWSO algorithm is applied to solving the optimal power flow of the modified IEEE 30-bus and 57-bus power systems, considering several real-world scenarios.In this study, known metaheuristic algorithms like SSA, PSO, NOA, KOA, and WOA, as well as the original WSO are utilized to evaluate MWSO's performance.In this situation, the default parameters of these rivals are utilized as recommended by the designer of the algorithms.For a fair comparison of outcomes, each algorithm is performed thirty times, with a population size of 30 and a maximum number of iterations of 1000.Several statistical metrics, such as best, mean, standard deviation, rank, and Wilcoxon rank-sum p-values are employed in this study.Specifically, "Windows 10 (64bit)" is used for the OS, "CPU Core i7 with 16 GB of RAM" is the hardware setup, and "MatLab 2020b" is the analytical tool of choice. (58)

Mathematical validation
The improved WSO method demonstrates superior performance in finding optimal solutions for unimodal functions, namely F1, F3, F4, and F7, beating the competing optimization algorithms, as seen in Table 2.In other words, the MWSO algorithm surpasses the original algorithm in these test functions.Conversely, the KOA method exhibits the most unfavourable results.Concerning the type of high-dimensional multimodal test functions, namely F8 through F13, the MWSO algorithm tends to get trapped in local optimum solutions for F8, F12, and F13.The modified MWSO algorithm demonstrates superior performance compared to the standard WSO and other rivals in achieving global solutions for various fixed multimodal testing functions i.e., F14-F23.Furthermore, it can be noticed that the MWSO algorithm offers the best overall rank, exceeding all the efficient and recent competing algorithms.Therefore, the findings concluded that the MWSO algorithm yields improved outcomes in tackling such optimization issues.

Convergence curve
With a dimension of 30 for the unimodal and multimodal functions, Fig. 7 visually analyzes the convergence rate of the developed MWSO method across the CEC 2017 benchmark testing functions.As can be observed in this figure, the proposed algorithm converges more quickly than other methods, particularly for the functions F1, www.nature.com/scientificreports/F3, F7, F4, F9, F10, and F11.The functions F2, F6, and F13 are examples of situations in which MWSO becomes trapped in a nearly optimal state.Additionally, the developed algorithm exceeds the original one for evaluating the best optimal solutions for all testing functions.In any case, the modified optimizer often yields better results with fewer iterations.Due to its quick convergence and improved accuracy, the MWSO approach is an efficient optimization tool for handling increasingly complex optimization scenarios.

Boxplot behaviour
Figure 8 depicts the boxplot curves derived from the MWSO optimizer and its competing counterparts.The figure visually represents the distribution of data across all functions from the CEC 2017 dataset.The whiskers on the boxplots represent the minimum and maximum values achieved by the algorithms.A tight box plot is indicative of a significant level of data consensus.Specifically, the MWSO method exhibits a lack of outliers throughout a set of more than ten functions, namely F7, F9, F10, F11, F16, F17, F19, F21, F22, and F23.Upon evaluating the boxplots of the majority of the testing functions, it becomes evident that the MWSO optimizer has a superior distribution characterized by lower values.The MWSO approach continuously exhibits better performance compared to other existing optimization methods, hence confirming its enhanced usefulness.

p-value-based statistical analysis
The statistical significance of the findings acquired by the MWSO algorithm and other competing algorithms is assessed via the use of the Wilcoxon rank-sum test.This test is employed to demonstrate that the observed performance was not only attributable to random chance.The analysis is carried out using a substantial threshold of 5% for all testing functions.Table 3 presents a summary of the outcomes obtained from the MWSO algorithm compared to the competing algorithms, as evaluated by the Wilcoxon test.It can be noticed from this table that the developed algorithm differs significantly from the other optimizers, in which the p-value is less than the significance level of 0.05, indicating that MWSO outperforms the original WSO and the others in terms of attaining optimum solutions and a higher convergence rate.Furthermore, the results of the 23 functions are shown in Table 4, using ANOVA analysis, the Friedman test, and the Kruskal test.According to the table, the MWSO demonstrates significantly greater efficacy in comparison to the six competing alternatives, as shown by p-values below 0.05.

Application of MWSO for OPF problem
In this subsection, the performance of the MWSO algorithm is compared to the performance of the original WSO in several real-world scenarios to determine whether the proposed algorithm is more successful at solving the OPF problem.
The following mathematical model can be used to express the OPF problem to be solved by the MWSO: where x and u denote the dependent (state) variables and the independent (control) variables, respectively.While, F(x, u) represent the objective functions of the OPF.The objective functions are constrained by set of equality constraints which are represented by g(x, u) and set of inequality constraints which are represented by h(x, u) , as previously presented in "Objective functions and system constraints" section.The control variables of the IEEE 30 bus system are considered the scheduled power of the thermal generators except the swing generator (at bus Vol:.( 1234567890) Table 3. Wilcoxon Rank test of MWSO vs compared methods for CEC2017.www.nature.com/scientificreports/1), the scheduled output power of the two wind plants, the scheduled output power of the solar PV plant, and the voltages of all generator buses, while the control variables of the IEEE 57 bus system are similar to the IEEE 30-bus system in addition to the reactive powers of the shunt compensators and the tap changer steps of branch transformers.The cases from 1 to 8 are conducted on the modified IEEE 30 bus power network, while Cases 9 and 10 are dedicated for solving the OPF problem in the IEEE 57 bus power network.The simulation process for these real -world cases is achieved through using the MATLAB software.The objective function of this case is to minimize the total production cost from all power sources in the system.The formulation of this objective is based on (22).The values of input parameters required for this case are summarized in Supplementary Material Tables 1A and 2A, while all simulation findings are recorded in Table 5.In comparison to the outcomes of the other techniques utilized, it was discovered through analysis of the findings in Table 5 that the MWSO produced power at the lowest cost, which came to $/h 781.6393.In addition, it was found that the suggested technique has the best convergence for the solution weighed against the WSO, as shown  www.nature.com/scientificreports/ in Fig. 9. Furthermore, looking at the control variables' limits and the network constraints, all values are within the acceptable limits, as indicated in Table 5.

Case#2: changing the value of reserve cost coefficient
In Case #1, the reserve cost coefficient for the two wind plants and the solar plant was constant value equals to 3 for all of them.In this case, the value of this coefficient will be changed from 3 to 5, 6, and 7 to study the effect of this change on the optimal cost of production.The value of penalty charge coefficient for both wind and solar plants is constant in this case at 1.5.The optimal schedule of output power of all generators is determined at each value of the reserve cost coefficient as a subcase.This optimal schedule is highlighted by Fig. 10.As anticipated, an escalation in the reserve charge coefficient led to a drop in the planned output of wind and solar power facilities.This drop can be explained as decreasing the schedule of renewable power will decrease the reserve charge in the event of overestimation.In contrast, the schedule of thermal power will increase due to reducing the schedule of renewable power.Consequently, the cost of production from renewable energy will decrease, while the cost of production from thermal generators will increase and the total cost of production from all generators will increase as indicated in Fig. 11.

Case#3: changing the value of penalty cost coefficient
This case is similar to Case 2, the only difference is varying the penalty cost coefficient with maintain the reserve cost coefficient constant to study the influence of changing the penalty charge coefficient on the schedule of power in the system and the associated costs of production.The values of the penalty cost coefficient are changed from 1.5 in Case 1 to be 2.5, 3.5, and 4.5, respectively.Each value from these values is considered as a subcase, and the results of schedule power and production cost is obtained.To analyse these results, the schedule powers from all generators are illustrated in Fig. 12.As anticipated, with escalating the penalty charge coefficient, the schedule power from renewable energy resources increases to minimize the penalty fees in the event of underestimation of renewable power.This increase in the schedule of renewable energy resulted in reducing the schedule power from thermal plants.This change in the scheduled powers will consequently be translated into the cost of production as shown in Fig. 13, where the cost of wind and solar power increases, while the cost of thermal power decreases, but the total cost of production will increase.

Case#4: impact of forcing carbon tax
The MWSO is applied to examine the impact of placing a tax on emissions from thermal energy generation in this scenario.The objective function is minimizing the total production cost with the existing carbon tax based on (23).All input parameters are set to the same values as in Supplementary Materials Tables 1A and 2A, except for the carbon tax, which is set at $20 per tonne.The purpose of imposing a tax on emissions is to reduce energy production from thermal sources and increase reliance on renewable energy sources.To ensure that the imposition of this tax achieved its goal, the results of this case, which are listed in Table 5, were examined, and it was observed that production from thermal energy sources was actually reduced while production from renewable energy sources increased compared to the first case in which no tax was imposed.As in the previous case study, the MWSO achieved the lowest production cost ($/h 810.3348) with the fastest solution convergence, as shown in Fig. 14 as well as all values for constraints inside the acceptable range.This case study was assigned to employ the suggested strategy (MWSO) according to Eq. ( 24) to lessen emissions since the system under investigation uses three thermal energy sources that emit a significant amount of greenhouse gases.In this situation, lowering emissions is the main objective, regardless of the cost of production.Therefore, it is evident from Table 5 that the emissions are minimized, while the total cost increases compared to Case#1.It is also noted that MWSO has outperformed the original WSO in minimizing the carbon emissions and convergence characteristics as illustrated in Fig. 15.
Case#6: minimizing power losses Another important objective of OPF is minimizing the active power losses.This objective is performed in this case study according to (26).The obtained result of this case is indicated in Table 5 and Fig. 16.It is observed from these outcomes that the minimum power loss is achieved by the MWSO with fast convergence compared to www.nature.com/scientificreports/method (MWSO) in this situation.This helps assess how well the proposed method works for figuring out the OPF problem in some complex scenarios that include changes in both the source and the load.For modelling the uncertainty of load demand, a normal PDF is used 61 as shown in Fig. 20.The selected values of the standard deviation (σ ld ) and the mean (µ ld ) for the normal PDF are 10 and 70, respectively.Each loading level (scenario) has a probability of occurrence, this probability ( ld,i ) can be calculated as follows: Here, P ld denotes the system loading, while P high ld,i and P low ld,i denotes the upper and lower limits of the loading level.While the mean of occurring a certain loading level ( P − ld,i ) can be determined as follows: The estimated means (in percentages of nominal system loading, P ld ) and the likelihoods for the four loading scenarios are indicated in Table 7.The outcomes of solving the OPF in this case using the WSO and the MWSO are listed in Table 8.The outcomes demonstrate once more how much more successful the MWSO is in this more complicated case when compared to the conventional WSO.For this case, the voltage profile of load buses through the four different loading scenarios is indicated by Fig. 21.It demonstrates that every voltage of the load buses is within the allowed range.

Case#9: minimization of total production cost in IEEE 57-bus system
This case study was created to determine the validity of the MWSO in tackling the problem of OPF in the most complex systems by minimizing the cost of power production in the standard IEEE-57 system.Based on (22), the objective function is the same as in the IEEE-30 bus system.The system restrictions are identical to those of the IEEE-30 bus system.The IEEE-57 bus system has been upgraded to include four thermal generators linked at buses 1 (swing), 3, 8, and 12, two wind plants linked at buses 2 and 6, and a solar PV plant linked at bus 9.The cost and emission coefficients of thermal generators in this system are detailed in Supplementary Material Table 3A, while the parameters of Weibull and lognormal PDF are provided in Supplementary Material Table 4A.The load of this system is 1250.8MW for active power and 336.4 MVA for reactive power.The simulation findings of this case are listed in Table 9.
The findings of this complicated case clearly prove the success of the MWSO in minimizing the total cost of production with high convergence characteristics, as shown in Fig. 22

Case#10: minimization of total production cost with carbon tax in IEEE 57-bus system
In this case, more complicated objective function is defined for minimizing the total cost with enforcing a tax on carbon emissions of thermal generator in the IEEE 57-bus system.The formulation of the objective function is the same as in (23).The parameters of this case are the same as in Case #9, and the carbon tax is set at $20 per tonne.This case is performed for 10 runs with 600 iterations for each run.The findings of this case are presented in Table 10.The convergence curves of the MWSO and WSO for this case are illustrated by Fig. 24.It also proofs www.nature.com/scientificreports/ that the MWSO has minimized the total cost with fast convergence compared to the original WSO.The voltages of the load buses are also within the allowed limits as shown by Fig. 25.

Statistical analysis
In this section, a statistical summary is presented in Table 11 for the case studies from 1 to 10 among the implemented simulation runs for each case.In addition to that, Wilcoxon signed rank test is carried out to compare between the MWSO and WSO as presented in Table 12.The column H O in this table specifies whether or not the null hypothesis is correct.The effectiveness of the two algorithms is statistically similar for the study instance if the null hypothesis is true (i.e., H O = Yes, with a threshold of significance = 0.05).
R+ is the sum of the rankings for runs in which MWSO exceeds WSO, while R− denotes the rankings for runs in which WSO exceeds MWSO.The p-value establishes the importance of results.The lower the p-value,

0Figure 1 .
Figure 1.Wind speed distribution and Weibull fitting for wind plant at bus 5.

Figure 2 .
Figure 2. Wind speed distribution and Weibull fitting for wind plant at bus 11.

Figure 3 .
Figure 3. Distribution of irradiance and lognormal fitting for solar PV at bus 13.

Figure 5 .
Figure 5.The complete flowchart of the original WSO algorithm.

Figure 6 .
Figure 6.The complete flowchart of the proposed MWSO algorithm.

F15Figure 8 .
Figure 8. Box plot of the modified MWSO and other rivals 23rd testing functions.

Figure 11 .
Figure 11.Impact of various reserve cost coefficients (RCC) on various costs.

Figure 12 .Figure 13 .
Figure 12.Optimal scheduled active powers of all generators with different penalty cost coefficient.

Figure 20 .
Figure 20.Normal PDF of network loading.

Available real power (MW) from solar PV at bus 13 0 0.05 0.1 0.15 0.2 0.25 Relative frequency Schedule solar power Figure 4. Available
Th (P Th ) + real solar PV power at bus 13.

Table 1 . The pseudo code of the MWSO algorithm
Optimal solutions.3:Initializethe parameters of WSO.4: Initialize the population using Eq.(39).5: Initialize the velocity of the initial populations.6: Evaluate the position of the initial populations.

Table 2 .
Results of the MWSO and comparison against other competitors on the CEC2017 benchmark.*Bold face highlights the best obtained solutions.

Table 4 .
Outcomes from the ANOVA, Friedman, and Kruskal tests.

Table 7 .
Means and probabilities of different loading scenarios.profile of load buses of the IEEE 57 bus network is indicated by Fig.23.It shows that all load buses voltages are within the allowed values.

Table 12 .
Results of Wilcoxon signed rank test.