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

Abstract

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.

Introduction

The optimization of power flow was initially developed by Carpentier in 1962¹. Afterward, several methods have been developed for addressing the optimal power flow (OPF) problem. The OPF is used to minimize power losses, maintain voltage stability, optimize generating costs, and eliminate gas emissions. The physical limits of the power network, which include the need to comply with power generator capability, buses’ voltage, capacities of transmission lines, power cable flows, and any other technical requirements, typically constrain this optimization. This may seem like a complex issue, particularly in high-power systems. Therefore, particular measures should be taken to prevent exceeding these physical boundaries. The classical OPF only consists of fossil fuel-fired conventional generating sources, which creates an exceedingly mixed integer, non-linear, and non-convex optimization issue^2,3,4. The increasing inclusion of renewable energy into electrical networks necessitates the inclusion of its uncertain character in OPF studies because of the accompanying issues throughout the operational and planning stages. Several traditional optimization methods have been developed to deal with the OPF challenge. These techniques include quadratic programming, non-linear programming, mixed-integer linear programming, and interior-point techniques^5,6. Certain strategies have been successfully employed in the industry because of their quick convergence and ability to provide the optimum solution. However, these optimization techniques necessitate first linearizing the optimization function. On the other hand, some heuristic optimization strategies have been proposed as a potential solution to address this issue⁷. For this reason, the OPF is solved using a variety of heuristic approaches.

The OPF problem was addressed using a sequential GA solution approach in combination with a simple genetic algorithm (SGA) to acquire a suitable control variable resolution without violating system constraints⁸. In Ref.⁹, a dependable and effective Tabu search best approach has been proposed and assessed on the IEEE 30-bus power network. Numerous earlier studies have relied on differential evolution to rectify the OPF problem. These studies have quick convergence characteristics and are appropriate for OPF problems with complex variables. Nevertheless, there is a significant chance that they will converge to a local instead of a global optimal solution^10,11,12. In several challenging OPF issues, particle swarm optimization has been applied. The premature convergence is a major disadvantage of classical PSO, as it is with many heuristic techniques^13,14,15. Grey wolf optimization¹⁶, artificial bee colony¹⁷, flower pollination algorithm¹⁸, cuckoo search optimization¹⁹, crow search algorithms²⁰, success history-based adaptive differential evolution algorithm²¹, group search optimization²², JAYA algorithm²³, moth swarm algorithm²⁴, golden ratio optimization method²⁵, and Aquila optimizer²⁶, barnacle mating optimizer²⁷, mayfly algorithm²⁸, coronavirus herd immunity algorithm²⁹, and weighted mean of vectors (INFO) algorithm³⁰ are just a few of the meta-heuristic population-based algorithms that have been employed in the past few years to solve OPF problems.

Various adjustments have been made to metaheuristic optimization approaches in the literature to address the issue of early convergence and provide an improved solution for the OPF problem, such as modified JAYA³¹, enhanced bacteria foraging algorithm (MBFA)³², SHADE-SF³³, modified grasshopper optimization³⁴, improved rao-2 algorithm³⁵, boosted quasi-reflection jellyfish optimization algorithm³⁶, hybrid cross entropy-cuckoo search algorithm³⁷, hybrid TLTFWO³⁸ through the integration between the teaching and learning algorithm and turbulent flow of water algorithm, hybrid Mayfly algorithm and Aquila optimizer³⁹. Accordingly, this study aims to develop a recent optimization technique named white shark optimization (WSO) to tackle the OPF, considering several real-world scenarios and the uncertainties associated with the generation and demand.

Freshly, the WSO algorithm was developed by Malik et al.⁴⁰ in 2022 and applied for handling most complex optimization challenges, such as solving uncertain optimal power flow⁴¹ and distributed generation optimal allocation⁴². However, WSO has some drawbacks, such as a slow convergence rate and an imbalance between the exploration and exploitation phases. In the literature, some studies have been conducted to overcome such limitations. In Ref.⁴³, the authors proposed a method for adjusting the force control parameters of the WSO by including a chaotic generator to enhance the exploitation capabilities of the algorithm. Further, the authors in Ref.⁴⁴ provided a suggested methodology involves adjusting the probability parameters of WSO to align with the optimization process and effectively synchronize all phases of the algorithm’s search process. Furthermore, they incorporated wave theory to elucidate the equation governing the trajectory motion of fluid particles inside the micro amplitude wave theory. The exploration process is also enhanced by incorporating the spiral search technique from the whale optimization algorithm. In Ref.⁴⁵, the authors proposed a new hybrid WSO and whale optimization algorithm to improve the stochastic behavior of the WSO algorithm for specifying the appropriate parameters of Li-ion battery Shepherd model equivalent circuits. Also, WSO is hybridized with the equilibrium optimizer for utilizing IOT for power scheduling problems⁴⁶. In this work, the Gaussian barebones (GB) and quasi-oppositional-based learning (QOBL) strategies are incorporated into the original WSO algorithm to enhance its convergence speed and accuracy while addressing the complicated optimal power flow problem.

The developed MWSO is evaluated via 23 benchmark functions, which include unimodal, high-dimensional multi-modal, and fixed high-dimensional multi-modal functions, and a comparison with other six rivals is conducted using different statistical analysis. These algorithms comprise particle swarm optimization (PSO)⁴⁷, whale optimization algorithm (WOA)⁴⁸, salp swarm algorithm (SSA)⁴⁹, Kepler optimization algorithm (KOA)⁵⁰, nutcracker optimizer algorithm (NOA)⁵¹, and the traditional WSO. Then, the MWSO algorithm is employed to solve the optimal power flow problem on the modified IEEE 30-bus and 57-bus power networks, considering different real-world scenarios. Eventually, the key effort of this research can be listed as follows:

• 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 and discussions of the 23 benchmark testing functions. Finally, “Conclusion” section concludes the findings and future recommendations of the paper.

Different cost models

In this study, some modifications are applied on the IEEE 30-bus test network to include wind and solar plants. At buses 5 and 11, the two thermal plants have been replaced by two wind power plants, and the thermal plant at bus 13 has also been replaced by a solar PV plant³³. The IEEE 57-bus system is also reformed by changing the thermal plants at buses 2 and 6 with two wind plants and changing the thermal plant at bus 9 with solar PV plant⁵². The data of the wind and solar plants of the IEEE 30-bus system and the IEEE 57-bus system are provided in Supplementary Material Tables 1A and 4A, respectively. This section will provide a detailed explanation of the production costs of each power source in the IEEE 30-bus power system. Since the production costs of the IEEE 57-bus power system follow the same procedure, they will not be explained here.

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.

$${C}_{Th}({P}_{Th})=\sum \limits_{i=1}^{{N}_{Th}}{a}_{Thi}+{b}_{Thi}{P}_{Thi}+{c}_{Thi}{{P}^{2}}_{Thi} +\left|{d}_{Thi}\times {\text{sin}}({e}_{Thi}\times ({P}_{Thi}^{min}-{P}_{Thi}))\right|,$$

(1)

where ${P}_{Thi}$ is the output power of the $i\text{-th}$ thermal plant, while ${a}_{Thi}$, ${b}_{Thi}$, and ${c}_{Thi}$ indicate the cost coefficients of the $i\text{-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}_{Thi}^{min}$ denotes the minimum amount of power produced from the $i\text{-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:

$${C}_{directwj}={dw}_{j}{P}_{schwj},$$

(2)

where ${P}_{schwj}$ represents the intended wind power of the $j\text{-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:

$${C}_{reserve\, wj}={K}_{reswj}\left({P}_{schwj}-{P}_{available\, wj}\right)={K}_{reswj}\int \limits_{0}^{{P}_{sch wj}}\left( {P}_{schwj}-{ P}_{wind\, j}\right){f}_{wind\, j}\left({P}_{wind\, j}\right)d{P}_{wind\, j},$$

(3)

where ${K}_{reswj}$ corresponds to the reserve cost coefficient for the $j\text{-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\text{-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):

$${C}_{penalty\, wj}={K}_{pen\, wj}({P}_{available\, wj}-{P}_{schwj})= {K}_{pen\, wj}\int \limits_{{P}_{schwj}}^{{P}_{rated \,wj}}({P}_{wind\, j}-{P}_{schwj}){f}_{wind\, j}({P}_{wind\, j})d{P}_{wind\, j},$$

(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})$:

$${f}_{{Wd}_{v}}\left({Wd}_{v}\right)=\left(\frac{k}{c}\right)+{\left(\frac{{Wd}_{v}}{c}\right)}^{\left(k-1\right)}{e}^{{-\left({{Wd}_{v}}/{c}\right)}^{k}}\,for\, 0<{Wd}_{v}<\infty ,$$

(5)

where the letters $k$ and $c$, respectively, stand for scale and form factors. Weibull distribution’s mean is calculated as follows:

$${M}_{weibull}=c\times\Gamma \left(1+{k}^{-1}\right).$$

(6)

The gamma function, which is represented by the sign $\Gamma$ in (6), is provided by:

$$\Gamma \left(x\right)=\int \limits_{0}^{\infty }{e}^{-t}{t}^{x-1}dt.$$

(7)

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.

Figure 1

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

Figure 2

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

The wind speed influences the wind plant’s output power of. According to Ref.³³, the following is the formula for wind turbine power output:

$${P}_{wind}\left({Wd}_{v}\right)=\left\{\begin{array}{*{20}l}0\, for\, {Wd}_{v}<{Wd}_{vin} \,and \,{Wd}_{v}>{Wd}_{vout} & \quad (discrete \,region )\\ {P}_{rated w}\times \left(\frac{({Wd}_{v}-{Wd}_{vin}}{({Wd}_{vr}-{Wd}_{vin}}\right)\,for\, {Wd}_{vin}\le {Wd}_{v}\le {Wd}_{vr} & \quad \left(continous\, region \right)\\ {P}_{rated\, w} \,for\, {Wd}_{vr}\le {Wd}_{v}\le {Wd}_{vout} &\quad \left(discrete \,region \right)\end{array}\right..$$

(8)

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}_{rated w}$.

It is possible to establish the probability of output power from wind plant in the discrete region as follows⁵⁴:

$${f}_{wind}\left({P}_{wind}\right)\left\{{P}_{wind}=0\right\}=1-{\text{exp}}\left[-{\left(\frac{{Wd}_{vin}}{c}\right)}^{k}\right]+{\text{exp}}\left[-{\left(\frac{{Wd}_{vout}}{c}\right)}^{k}\right],$$

(9)

$${f}_{wind}\left({P}_{wind}\right)\left\{{P}_{wind}={P}_{rated w}\right\}={\text{exp}}\left[-{\left(\frac{{Wd}_{vr}}{c}\right)}^{k}\right]-{\text{exp}}\left[-{\left(\frac{{Wd}_{vout}}{c}\right)}^{k}\right].$$

(10)

Regarding the continuous zone, the following formula can be used to determine the probabilities for the power that the wind plant will produce:

$${f}_{wind}\left({P}_{wind}\right)=\frac{k\left({Wd}_{vr}-{Wd}_{vin}\right)}{{c}^{k}\times {P}_{rated w}}{\left[{Wd}_{vin}+\frac{{P}_{wind}}{{P}_{rated w}}\left({Wd}_{vr}-{Wd}_{vin}\right)\right]}^{k-1}\times {\text{exp}}\left[{-\left(\frac{{Wd}_{vin}+\frac{{P}_{wind}}{{P}_{rated w}}\left({Wd}_{vr}-{Wd}_{vin}\right)}{c}\right)}^{k}\right].$$

(11)

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).

$${C}_{directsk}={ds}_{k}{P}_{schsk},$$

(12)

where ${P}_{schsk}$ represents the intended solar power of the $k\text{-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:

$${C}_{reserve\, sk}={K}_{ressk}\left({P}_{schsk}-{P}_{available \,sk}\right)= {K}_{ressk}\times {f}_{solar k}\left({P}_{available \,sk}<{P}_{sch sk}\right)\times \left[({P}_{schsk}-E({P}_{available \,sk}<{P}_{schsk})\right].$$

(13)

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}_{solar k}({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:

$${C}_{penalty sk}={K}_{pen sk}\left({P}_{available \,sk}-{P}_{schsk}\right)= {K}_{pensk}\times {f}_{solar k}\left({P}_{available \,sk}>{P}_{schsk}\right)\times \left[(E({P}_{available \,sk}<{P}_{sch sk})-{P}_{schsk}\right],$$

(14)

where, ${K}_{pensk}$ signifies the penalty cost coefficient, ${f}_{solar k}({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 remaining 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⁵⁵.

$${f}_{I}\left(I\right)=\frac{1}{I\sigma \sqrt{2\pi }}{\text{exp}}\left\{\frac{-{\left({\text{ln}}x-\mu \right)}^{2}}{{2\sigma }^{2}}\right\}for I>0.$$

(15)

As, the irradiance probability is denoted by ${f}_{I}(I)$, the mean of the lognormal distribution is denoted by $\mu$, while the standard deviation is signified by $\sigma$, respectively. While the mean of lognormal, ${M}_{lgn}$ is calculated by:

$${M}_{lgn}={\text{exp}}\left(\mu +\frac{{\sigma }^{2}}{2}\right).$$

(16)

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.

Figure 3

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

Consequently, the sun irradiation vs. the energy conversion of the solar plant can be presented as follows⁵⁶:

$${P}_{solar}\left(I\right)=\left\{\begin{array}{l}{P}_{solarr}\left(\frac{{I}^{2}}{{I}_{std}{R}_{c}}\right)\, for\, 0<I<{R}_{c}\\ {P}_{solarr}\left(\frac{I}{{I}_{std}}\right)\, for\, I\ge {R}_{c}\end{array}\right.,$$

(17)

where in this formula, ${I}_{std}$ signifies the solar irradiance when the environment is a standard i.e. (800 W/m²), the symbol ${R}_{c}$ specifies a specific value of irradiance (120 W/m²), 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:

$${C}_{reserve\, sk}={K}_{ressk}\left({P}_{schsk}-{P}_{available \,sk}\right)= {K}_{ressk}\times \sum_{n=1}^{{N}^{-}}\left[{P}_{schsk}-{P}_{sn-}\right]\times {f}_{sn-},$$

(18)

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.

Figure 4

Available real solar PV power at bus 13.

While, the penalty cost that are previously stated in (14) can be rewritten as follows:

$${C}_{penalty sk}={K}_{pensk}\left({P}_{available \,sk}-{P}_{schsk}\right)={K}_{pensk}\sum_{n=1}^{{N}^{+}}\left[{P}_{sn+}-{P}_{schsk}\right]\times {f}_{sn+},$$

(19)

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}_{sch sk}$) 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 coefficient 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 (${\text{CAE}})$, thus Eq. (20) provides an estimation for these emissions.

$$CAE=\left(\sum \limits_{i=1}^{{N}_{Th}}\left({\varphi }_{Th,i}+{\Psi }_{Th,i}\;{P}_{Thi}+{\omega }_{Th,i}\;{{P}^{2}}_{Thi}\right)+{\tau }_{Th,i}\;{e}^{({\zeta }_{Th,i}\;{P}_{Thi})}\right).$$

(20)

Here, ${\varphi }_{Th,i}$, ${\Psi }_{Thi}$, ${\omega }_{Th,i}$, ${\tau }_{Th,i}$, and ${\zeta }_{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.

$${C}_{CE}=CAE\times {C}_{Taxc}.$$

(21)

Objective functions and system constraints

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:

$${F}_{1} =\text{ min}\left({C}_{Th}\left({P}_{Th}\right)+\sum \limits_{j=1}^{{N}_{WP}}\left({dw}_{j}{P}_{schwj}+{K}_{reswj}\left({P}_{schwj}-{P}_{available\, wj}\right)+{K}_{penwj}({P}_{available\, wj}-{P}_{schwj})\right)+\sum \limits_{j=1}^{{N}_{SP}}\left({ds}_{k}{P}_{schsk}+{K}_{ressk}\left({P}_{schsk}-{P}_{available \,sk}\right)+{K}_{pensk}({P}_{available \,sk}-{P}_{schsk})\right)\right).$$

(22)

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:

$${F}_{2}=min\left({C}_{Th}\left({P}_{Th}\right)+\sum \limits_{j=1}^{{N}_{WP}}\left[{dw}_{j}{P}_{schwj}+{K}_{reswj}\left({P}_{schwj}-{P}_{available\, wj}\right)+{K}_{penwj}\left({P}_{available\, wj}-{P}_{schwj}\right)\right]+\sum \limits_{j=1}^{{N}_{SP}}\left[{ds}_{k}{P}_{schsk}+{K}_{ressk}\left({P}_{schsk}-{P}_{available \,sk}\right)+{K}_{pensk}\left({P}_{available \,sk}-{P}_{schsk}\right)\right]+\left[CAE\times {C}_{Tax}\right]\right).$$

(23)

Minimizing the carbon emissions $({F}_{3})$

${F}_{3}$ Is to minimize the total carbon emissions of the thermal plants in (20). Therefore, it can be written as:

$${F}_{3}=min\left(\sum \limits_{i=1}^{{N}_{Th}}\left[{\varphi }_{Th,i}+{\Psi }_{Th,i}\;{P}_{Thi}+{\omega }_{Th,i}\;{{P}^{2}}_{Thi}+{\tau }_{Th,i}\;{e}^{({\zeta }_{Th,i}\;{P}_{Thi})}\right]\right).$$

(24)

Minimizing the total power losses $({F}_{4})$

The power losses of the power system can be formulated as follows:

$${P}_{loss}=\left(\sum_{i=1}^{NTL}\sum_{j\ne 1}^{NTL}\left[{G}_{ij}{V}_{i}^{2}+{V}_{j}^{2}-2{V}_{i} {V}_{j}{\text{cos}}({\delta }_{ij})\right]\right),$$

(25)

${\text{where }\delta }_{ij}$ expresses the difference in voltage angles between buses $i$ and $j$; $NTL$ signifies the number of transmission lines; and ${G}_{ij}$ expresses the transfer conductance. Consequently, ${F}_{4}$ can be formulated as follows:

$${F}_{4}=min\left(\sum_{i=1}^{NTL}\sum_{j\ne 1}^{NTL}\left[{G}_{ij}{V}_{i}^{2}+{V}_{j}^{2}-2{V}_{i} {V}_{j}{\text{cos}}({\delta }_{ij})\right]\right).$$

(26)

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:

$$P_{Gi} = P_{Di} + V_{i} \mathop \sum \limits_{j = 1}^{NB} V_{j} \left[ {G_{ij} {\text{cos}}\left( {\delta_{ij} } \right) + B_{ij} {\text{sin}}\left( {\delta_{ij} } \right)} \right]{ }\quad \forall i \in NB{ },$$

(27)

$$Q_{Gi} = Q_{Di} + V_{i} \mathop \sum \limits_{j = 1}^{NB} V_{j} \left[ {G_{ij} {\text{sin}}\left( {\delta_{ij} } \right) - B_{ij} {\text{cos}}\left( {\delta_{ij} } \right)} \right]\quad \forall i \in NB,$$

(28)

where $NB$ stands for the network buses number. ${{B}_{ij}\text{ and }G}_{ij}$ stand for the susceptance and conductance among bus $i$ and bus $j$, respectively. ${\delta }_{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}\, {\text{and}}\, {P}_{Di}$, respectively, while the reactive components of the consumed and generated power are expressed by ${Q}_{Di}{ \,{\text{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), and (31), respectively, where ${N}_{G}$ represents the number of network’s generators.

$${P}_{Gi}^{min}\le {P}_{Gi}\le {P}_{Gi}^{max}, i=1,\ldots ,{N}_{G},$$

(29)

$${Q}_{Gi}^{min}\le {Q}_{Gi}\le {Q}_{Gi}^{max}, i=1,\ldots, {N}_{G},$$

(30)

$${V}_{Gi}^{min}\le {V}_{Gi}\le {V}_{Gi}^{max}, i=1,\ldots ,{N}_{G}.$$

(31)

Transformer limits

$${T}_{t}^{min}\le {T}_{t}\le {T}_{t}^{max}, t=1,\ldots ,{N}_{T}.$$

(32)

Shunt compensator limits

$${Q}_{Cc}^{min}\le {Q}_{Cc}\le {Q}_{Cc}^{max}, c=1,\ldots ,{N}_{C},$$

(33)

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:

$${{{\text{P}}}_{{\text{Thi}}}-{\text{P}}}_{{\text{Thi}}}^{0}\le {{\text{Ur}}}_{{\text{i}}},\text{ if power generation rises },$$

(34)

$${{\text{P}}}_{{\text{Thi}}}^{0}-{{\text{P}}}_{{\text{Thi}}}\le {{\text{Dr}}}_{{\text{i}}},\text{ if power generation reduces},$$

(35)

where,${P}_{Thi}^{0}$ signifies the previous hour’s output power from the thermal plant $( i )$. ${Ur}_{i} {\text{and}} {Dr}_{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:

$${{\text{V}}}_{{{\text{LB}}}_{{\text{p}}}}^{{\text{min}}}\le {{\text{V}}}_{{{\text{LB}}}_{{\text{p}}}}\le {{\text{V}}}_{{{\text{LB}}}_{{\text{p}}}}^{{\text{max}}},\text{ p}=1,\dots .\dots \dots ,{{\text{N}}}_{{\text{LB}}},$$

(36)

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:

$${V}_{d}=\sum_{p=1}^{{{\text{N}}}_{{\text{LB}}}}\left|{V}_{LBp}-1\right|.$$

(37)

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.

$${S}_{{L}_{q}}\le {S}_{{L}_{q}}^{max} , q=1,\ldots ,{N}_{L}.$$

(38)

Proposed MWSO methodology

Original WSO algorithm: an overview

This section clarifies a short description of the mathematical formulation of the original WSO⁴⁰, 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).

$$w=\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{l}{w}_{1}^{1}\\ {w}_{1}^{2}\end{array}& \begin{array}{l}{w}_{2}^{1}\\ {w}_{2}^{2}\end{array}\end{array}& \begin{array}{l}\cdots \\ \dots \end{array}& \begin{array}{l}{w}_{d}^{1}\\ {w}_{d}^{2}\end{array}\\ \begin{array}{cc}\vdots & \vdots \end{array}& \ddots & \vdots \\ \begin{array}{cc}{w}_{1}^{3}& {w}_{2}^{3}\end{array}& \cdots & {w}_{d}^{n}\end{array}\right],$$

(39)

where $w$ symbolizes the position of all sharks in the searching space, $d$ indicates the number of a chosen variables 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.

$${v}_{k}^{j}=\mu \left[{v}_{k-1}^{j}+{p}_{1}\left({w}_{gbest,k-1}-{w}_{k}^{j}\right)\times {c}_{1}+{p}_{2}\left({w}_{best}^{{v}_{k-1}^{j}}-{w}_{k-1}^{j}\right)\times {c}_{2}\right],$$

(40)

where $j=\text{1,2},\dots ,n$, represents the white shark’s index for $n$ search agents; the new velocity of the jth white shark in ${(k-1)}{\text{th}}$ step is denoted by ${v}_{k-1}^{j}$; ${w}_{gbest,k-1}$ is the global best solution obtained so far in the iteration ${(k-1)}{\text{th}}$; ${w}_{best}^{{v}_{k-1}^{j}}$ symbolizes the jth best known position for the swarm and ${v}_{k-1}^{j}$ 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 governing the influence of ${w}_{gbest,k-1}$ and ${w}_{best}^{{v}_{k-1}^{j}}$ on ${w}_{k-1}^{j}$ are represented by ${p}_{1}$ and ${p}_{2}$, respectively, which are computed by Eqs. (42) and (43); eventually, $\mu$ is the constriction factor to represent the convergence characteristics of the white sharks, formulated as per Eq. (44).

$$v=\left[n\times rand(1,n)\right]+1,$$

(41)

$${p}_{1}={p}_{ub}+\left({p}_{ub}-{p}_{lb}\right)\times {e}^{-{\left(\frac{4a}{A}\right)}^{2}},$$

(42)

$${p}_{2}={p}_{lb}+\left({p}_{ub}-{p}_{lb}\right)\times {e}^{-{\left(\frac{4a}{A}\right)}^{2}}.$$

(43)

Here, $a$ and $A$ are the current and maximum iterations, respectively; ${p}_{lb}$ and ${p}_{ub}$ are respectively the minimum and maximum velocity to accomplish good movement for the white sharks, and their values are 0.5 and 1.5, respectively.

$$\mu =\frac{2}{\left|2-\tau -\sqrt{{\tau }^{2}-4\tau }\right|},$$

(44)

where, $\tau$ 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.

$${w}_{k}^{j}=\left\{\begin{array}{*{20}l}{w}_{k-1}^{j}\to \oplus {w}_{o}+u.x+l.y & \quad rand<mv\\ {w}_{k-1}^{j}+\frac{{v}_{k-1}^{i}}{f} & \quad rand\ge mv\end{array},\right.$$

(45)

where ${w}_{k}^{j}$ 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); $\to$ 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).

$$x=sgn\left({w}_{k-1}^{j}-u\right)>0,$$

(46)

$$y=sgn\left({w}_{k-1}^{j}-l\right)<0,$$

(47)

$${w}_{o}=\oplus \left(x,y\right),$$

(48)

where $\oplus$ is a bit-wise XOR operation.

$$f={f}_{min}+\frac{{f}_{max}-{f}_{min}}{{f}_{max}+{f}_{min}} .$$

(49)

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.

$$mv=\frac{1}{{d}_{o}+{e}^{(A/2-a)/{d}_{1}}},$$

(50)

where ${d}_{o}$ and ${d}_{1}$ are respectively two constant positive numbers that utilized to manage the behaviour of exploration 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.

$${\stackrel{`}{w}}_{k}^{j}={w}_{gbest, k-1}+{r}_{1}{\overrightarrow{D}}_{w} sgn\left({r}_{2}-0.5\right) {r}_{3}<s,$$

(51)

where ${\stackrel{`}{w}}_{k}^{j}$ is the updated position of the jth white sharks with respect to the location of the prey; $sgn\left({r}_{2}-0.5\right)$ 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 ${\overrightarrow{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).

$${\overrightarrow{D}}_{w}=\left|rand\times \left({w}_{gbest, k-1}-{w}_{k-1}^{j}\right)\right|,$$

(52)

where $rand$ is a random number within the range [0,1]; ${w}_{k-1}^{j}$ represents the current position of the white shark in respect to the best position, ${w}_{gbest, k-1}$.

$$s=\left|1-{e}^{-\left(ga/A\right)}\right| ,$$

(53)

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:

$${w}_{k}^{j}=\frac{{w}_{k-1}^{j}+{\stackrel{`}{w}}_{k}^{j}}{2rand}.$$

(54)

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.

Figure 5

The complete flowchart of the original WSO algorithm.

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.

Table 1 The pseudo code of the MWSO algorithm.

Figure 6

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⁵⁷. 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:

$${w}_{k,GB}^{j}=\left\{\begin{array}{*{20}l}G\left(\frac{\left({w}_{gbest, k-1}+{w}_{k}^{j}\right)}{2},\left|{w}_{gbest, k-1}-{w}_{k}^{j}\right|\right) &\quad if \,rand<R\\ {w}_{k}^{j1}+{r}_{4}\left({w}_{k}^{j2}-{w}_{k}^{j3}\right)& \quad otherwise\end{array}\right.,$$

(55)

where ${w}_{k,GB}^{j}$ 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 $\left(\left({w}_{gbest, k-1}+{w}_{k}^{j}\right)/2\right)$ and standard deviation $\left({w}_{gbest, k-1}-{w}_{k}^{j}\right)$; $j1$, $j2$, and $j3$ are three randomly chosen individuals that are diverse from the current individual, $j$; ${w}_{k}^{j}$ 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⁵⁸, 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⁵⁹ to help clarify the notion of opposition-based learning.

Opposite number

If $z$ is a random number in the range $[{a}_{1}, {a}_{2}]$, then its opposite one may be written as:

$${z}^{o}={a}_{1}+{a}_{2}-z.$$

(56)

Opposite point

If $p({z}_{1}, {z}_{2},\dots ,{z}_{j},\dots ,{z}_{d})$ is a point in a searching space with a $d$-dimensional in which the ${z}_{j}\in [{a}_{1,j}, {a}_{2,j}]$, then its opposite one $op({z}_{1}^{o}, {z}_{2}^{o},\dots ,{z}_{j}^{o},\dots ,{z}_{d}^{o})$ is expressed by the following relation:

$${z}_{j}^{o}={a}_{1,j}+ {a}_{2,j}-{z}_{j}.$$

(57)

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⁶⁰. 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.

$${z}^{qo}=rand\left[\left(\frac{{a}_{1}+{a}_{2}}{2}\right),\left({a}_{1}+{a}_{2}-z\right)\right].$$

(58)

Quasi-opposite point

Similarly, the quasi-opposite point, $qop({z}_{1}^{qo}, {z}_{2}^{qo},\dots ,{z}_{j}^{qo},\dots ,{z}_{d}^{qo})$ in a searching region with a $d$-dimensional problem is calculated as per Eq. (59).

$${z}_{j}^{qo}=rand\left[\left(\frac{{a}_{1,j}+{a}_{2,j}}{2}\right),\left({a}_{1,j}+{a}_{2,j}-{z}_{j}\right)\right].$$

(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.

Experimental results of CEC 2017 testing functions

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.

Table 2 Results of the MWSO and comparison against other competitors on the CEC2017 benchmark.

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, 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.

Figure 7

Convergence characteristics of the developed HRSOAPOA and other competitors for CEC 2017 benchmark functions.

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.

Figure 8

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

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.

Table 3 Wilcoxon Rank test of MWSO vs compared methods for CEC2017.

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

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:

$${\text{Minimize}}:\text{ F}\left({\text{x}},{\text{u}}\right),$$

$$\text{Subject to}:\text{ h}\left({\text{x}},{\text{u}}\right)\le 0,$$

$${\text{g}}\left({\text{x}},{\text{u}}\right)=0,$$

where ${\text{x}}$ and ${\text{u}}$ denote the dependent (state) variables and the independent (control) variables, respectively. While, ${\text{F}}\left({\text{x}},{\text{u}}\right)$ represent the objective functions of the OPF. The objective functions are constrained by set of equality constraints which are represented by ${\text{g}}\left({\text{x}},{\text{u}}\right)$ and set of inequality constraints which are represented by ${\text{h}}\left({\text{x}},{\text{u}}\right)$, 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 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.

Case#1: minimizing the total cost

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 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.

Table 5 Findings of cases#1, 4, 5, and 6.

Figure 9

Case 1—solution convergence.

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.

Figure 10

Optimal scheduled active powers of all generators with different reserve cost coefficient.

Figure 11

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

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.

Figure 12

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

Figure 13

Impact of various penalty cost coefficients (PCC) on various costs.

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.

Figure 14

Case 4-solution convergence.

Case#5: minimizing carbon emissions

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.

Figure 15

Case#5-solution convergence.

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 the WSO. The voltage profile of load buses voltage for the Cases 1, 4, 5, and 6 is indicated by Fig. 17. The voltage profile shows that all voltages of load buses are within the allowed values.

Figure 16

Case#6-Solution Convergence.

Figure 17

Voltage profile of load buses—Cases#1, 4, 5, and 6.

Case#7: ramp rate of thermal generators

The limits of ramp rate for thermal generators can change the optimal solution of optimal power flow problem, thus this case study is dedicated to study their impact on the OPF problem. The input factors of this case are the same as in Supplementary Material Table 1A, while the output power at the preceding hour and the ramp rate limits for each thermal generator are indicated in Supplementary Material Table 2A. The simulation findings for this situation are indicated in Table 6. MWSO has achieved the lowest cost of production when compared to WSO, and it is faster to converge, as indicated in Fig. 18. Most notably, the total production cost increased from 781.6393 in Case 1 to 803.6681 in this case, as was to be expected, as the limits of operation of the thermal generators were changed in this case. Additionally, the voltage profile of load buses voltage for Case 7 is indicated by Fig. 19. It demonstrates that every voltage of the load buses is within the allowed range.

Table 6 Case#7’s findings.

Figure 18

Case#7—solution convergence.

Figure 19

Voltage profile of load buses—Case#7.

Case#8: uncertainty of load demand

Another important factor that may influence the solution of the optimal power flow problem is the uncertainty in the load demand, so this case study was dedicated to figuring out the OPF problem using the proposed 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⁶¹ as shown in Fig. 20. The selected values of the standard deviation $({\upsigma }_{{\text{ld}}})$ and the mean $({\upmu }_{{\text{ld}}})$ for the normal PDF are 10 and 70, respectively. Each loading level (scenario) has a probability of occurrence, this probability (${\Delta }_{{\text{ld}},{\text{i}}}$) can be calculated as follows:

Figure 20

Normal PDF of network loading.

$${\Delta }_{ld,i}=\int \limits_{{P}_{ld,i}^{low}}^{{P}_{ld,i}^{high}}\frac{1}{{\sigma }_{ld}\surd 2\pi }exp\left[\frac{{\left({P}_{ld}-{\mu }_{ld}\right)}^{2}}{{2{\sigma }_{ld}}^{2}}\right]d{P}_{ld} .$$

(60)

Here, ${P}_{ld}$ denotes the system loading, while ${P}_{ld,i}^{high}$ and ${P}_{ld,i}^{low}$ 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:

$${P}_{ld,i}^{-}=\frac{1}{{\Delta }_{ld,i}}\int \limits_{{P}_{ld,i}^{low}}^{{P}_{ld,i}^{high}}\left(\frac{1}{{\sigma }_{ld}\surd 2\pi }exp\left[\frac{{\left({P}_{ld}-{\mu }_{ld}\right)}^{2}}{{2{\sigma }_{ld}}^{2}}\right]\right)d{P}_{ld }.$$

(61)

The estimated means (in percentages of nominal system loading, ${{\text{P}}}_{{\text{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.

Table 7 Means and probabilities of different loading scenarios.

Table 8 Findings of Case#8.

Figure 21

Voltage profile of load buses—Case#8.

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.8 MW for active power and 336.4 MVA for reactive power. The simulation findings of this case are listed in Table 9.

Table 9 Findings of Case#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, compared to the original WSO. The voltage 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.

Figure 22

Case#9—solution convergence.

Figure 23

Voltage profile of load buses—Case#9.

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 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.

Table 10 Findings of Case#10.

Figure 24

Case#10—solution convergence.

Figure 25

Voltage profile of load buses—Case#10.

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).

Table 11 Statistical analysis.

Table 12 Results of Wilcoxon signed rank test.

$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, the stronger the argument against the null hypothesis (${H}_{O}$). The results show that MWSO outperformed WSO in all cases, as the p-value is lower than 0.05 and there are no null hypotheses.

Conclusion

This paper has introduced a modified white shark optimization (WSO) algorithm for optimizing power flow problems. The modified algorithm incorporates Gaussian barebones and quasi-oppositional learning mechanisms to improve its performance. The MWSO algorithm is tested using the CEC2017 benchmark functions and compared against six other efficient algorithms. The results show superior performance, making it well-suited for addressing power flow optimization problems, especially in renewable energy sources with intermittent output and fluctuating load demands. The paper introduces probabilistic models for solar and wind power using Weibull and lognormal PDFs, and presents a normal PDF-based probabilistic model for load demand. The MWSO algorithm is applied to solve the power flow optimization problem in two modified IEEE standard test systems. In the IEEE 30-bus system, it is used to minimize the total generation cost, both with and without considering carbon tax on emissions, while simultaneously minimizing active power losses. The study investigates the impact of varying reserve and penalty costs for overestimating and underestimating wind and solar power output, four different load scenarios, and the influence of imposing ramp rate limits of thermal generators on the optimal power flow problem. To validate the robustness of the proposed algorithm in more complex systems, the IEEE 57-bus network is also modified and subjected to the MWSO algorithm for solving the power flow optimization problem. The simulation results, statistical analysis, and the Wilcoxon signed rank test confirm the superiority and effectiveness of the MWSO algorithm in addressing power flow optimization problems.