Multi-objective Snow Ablation Optimization Algorithm: An Elementary Vision for Security-Constrained Optimal Power Flow Problem Incorporating Wind Energy Source with FACTS Devices

Abstract

This study delves into the exploration of a novel Multi-objective Snow Ablation Optimizer (MOSAO) algorithm, tailored for addressing expansive Optimal Power Flow (OPF) challenges inherent in intricate power systems. These systems are often complemented with the integration of renewable energy modalities and the state-of-the-art Flexible AC Transmission Systems (FACTS). Building upon the foundational framework of a previously documented single-objective Snow Ablation Optimizer, we have evolved it into the MOSAO paradigm. This transformation is achieved by harnessing the potency of non-dominated sorting coupled with the crowding distance strategy. The task of OPF magnifies in complexity when integrating renewable energy resources due to their inherent unpredictability and intermittent nature. As the modern power landscape evolves, FACTS devices are witnessing an increasing deployment to mitigate network demand and alleviate congestion issues. Within the ambit of this research, we've incorporated a stochastic wind energy source, working synergistically with an array of FACTS instruments. These encompass the static VAR compensator, thyristor-controlled series compensator and thyristor-driven phase shifter, all operating within the confines of an IEEE-30 bus framework. Strategic placement and calibration of these FACTS devices aim to optimize the system by minimizing the cumulative fuel expenditure. The capricious essence of wind as an energy source is elegantly depicted through the lens of Weibull probability density graphs. To distil the optimal middle-ground solutions, we've employed a fuzzy decision-making matrix. When benchmarking our findings against those derived from other esteemed optimization algorithms, we observe a notable distinction. The results from the modified IEEE-30 bus system accentuate the superior convergence, diversity and distribution attributes of MOSAO, especially when scrutinizing power flows. The MOSAO source code is available at: https://github.com/kanak02/MOSAO.

1 Introduction

Multi-objective optimization challenges, which are characterized by constraints, stand at the forefront of today’s problem-solving domain. Unlike their single-objective counterparts, which typically converge to a single optimal solution, multi-objective optimizations offer a diverse landscape of potential solutions. These solutions are often termed Pareto Fronts (PFs). A robust multi-objective optimization strategy should consistently yield solutions that span the range of these PFs, ensuring they represent feasible and optimal responses to the given problems. The complexity of these strategies lies in balancing numerous objectives. Traditionally, Metaheuristic (MH) algorithms are benchmarked against well-defined, straightforward optimization challenges. However, the unique and varied requirements of engineering design distinguish it from more conventional search problems. Tailoring and refining these algorithms to meet the intricate demands of engineering tasks is crucial for effective optimization.

In the realm of power systems, the optimal power flow (OPF) challenge has been a persistent issue for researchers. The integration of unconventional energy sources into the electrical grid in the modern era adds to this complexity [1,2,3]. Recent advancements in the field of unconventional energy sources, such as the development of magnetoelectric current sensors for wireless current measurement [4], further complicate the OPF challenge. A significant contribution was made by Ullah et al. [5], who introduced a hybrid strategy combining phasor particle swarm optimization with gravitational search, specifically targeting OPF issues in integrated wind and solar energy systems. This approach aligns with the trends in OPF challenges, where new technologies like the nine-phase open-end winding PMSMs are being integrated [6].

Elattar [7] developed a mathematical model to address the OPF problem within a combined heat and power framework, incorporating stochastic wind energy sources. This method's effectiveness was proven on an IEEE 30-bus test system, especially when compared to traditional techniques. Similarly, Man-Im et al. [8] explored an enhanced particle swarm optimization (PSO) approach, integrating chaotic mutation with stochastic weights, to tackle a multi-objective OPF problem centered on wind energy systems. Their research, which included IEEE 30- and 118-bus test systems, contributed significantly to the field. In parallel, studies [9,10,11] like have explored risk-constrained stochastic scheduling in energy hubs, highlighting the increasing complexity and multi-dimensionality of OPF problems. Salkuti [12] explored the glowworm swarm optimization technique to address multi-objective OPF challenges in advanced power systems utilizing wind energy. Extensive testing on IEEE 30- and 300-bus systems affirmed the value of this approach.

Kathiravan et al. [13], inspired by nature, applied the flower pollination algorithm to OPF problems in coal, wind, and solar energy matrices. Their extensive simulations on IEEE 30-bus and Indian utility 30-bus systems highlighted the versatility of their method. Duman et al. [14] employed the differential evolutionary particle swarm optimization (DEEPSO) for OPF challenges related to manageable wind and photovoltaic (PV) energy sources. Testing on IEEE 30-, 57-, and 118-bus systems revealed the robustness of the DEEPSO approach, especially when compared to other contemporary strategies. This team also emphasized integrating FACTS devices, such as the thyristor-controlled phase shifter (TCPS) and the thyristor-controlled series capacitor (TCSC), to solve the OPF problem. They introduced chaotic maps along with a refined version of the PSOGSA technique, yielding promising results [15]. Biswas et al. [16] presented a comprehensive solution to the OPF challenge that incorporates coal-based, wind, solar, and small-hydro energy sources. Their use of the multi-objective evolutionary algorithm on successive iterations of the IEEE 30-bus systems demonstrated the extensive potential and adaptability of their methodology.

Chen et al. [17] employed the Constrained Multi-objective Population Extremal Optimization methodology to tackle the OPF challenges associated with integrating wind and solar energy systems. The field also saw advancements through the Hybrid Modified Imperialist Competitive Algorithm paired with the SQP Algorithm [18], Success-History-Based DE emphasizing feasible solutions [19], Evolutionary PSO [20], and the Symbiotic Organisms Search [21]. Pandya and Jariwala's [22,23,24,25] multi-dimensional investigations into both single and composite objective OPF challenges, utilizing an array of sustainable energy sources, leveraged the latest metaheuristic tools. Biswas et al. [26] provided a comprehensive analysis of FACTS devices, wind turbines, and coal-based power systems using the SHADE technique.

As the “No Free Lunch” (NFL) theorem [27] reminds us that no single metaheuristic can solve all real-world problems. This realization has spurred the development and refinement of various metaheuristic methods [28, 29]. For instance, the multi-objective dragonfly algorithm [30] adapts its focus during different stages based on population evolution characteristics. Despite its effectiveness, it requires time-consuming calculations of two indicators for each selection, leading to variable running times. The multi-objective ions motion optimizer [31], known for its fast convergence and single indicator use, often results in poor population distribution and may not be effective for certain problems. Similarly, the multi-objective ant lion optimizer [32], effective in handling problems with uniform solution distributions on the Pareto front, struggles with problems having uneven distributions and faces limitations due to its simplistic approach to weight vector production.

The proposed Multi-objective Snow Ablation Optimization Algorithm (SAO) [33] presents two key advantages: efficient processing with fewer objectives and the use of secondary criteria, like crowding distance, for performance enhancement. However, these methods face challenges with increased objective numbers and generally slow convergence speeds. A critical challenge in the SAO is balancing convergence and divergence. SAO algorithms filter solutions based on quality criteria, with elite non-dominated solutions stored in distinct archives. Strategies for crafting these archives, often with a predefined maximum capacity, are a subject of ongoing scientific discourse. A notable feature of the SAO is its ability to autonomously determine the optimal population size and end the process. Despite its distinctive characteristics, the SAO faces both advantages and challenges in the evolving optimization landscape. Some of the pros of the SAO Algorithm are

• The SAO algorithm can achieve convergence swiftly, effortlessly marrying the concepts of exploration and exploitation.

• With its relaxed convergence mechanism, the SAO algorithm is less inclined to get entrapped in local minima.

• With a user-friendly design, the algorithm is easily adaptable to varied problem sets.

However the traditional SAO falters demonstrating oscillatory convergence trends or getting trapped in deceptive optima, especially high dimensional and/or multimodal problems.

As the energy landscape evolves, there is a growing need to intensify research across both traditional and alternative energy sources. In this regard, the present study makes the following contributions:

• 1. A concentrated effort on single and MO OPF problems that seamlessly amalgamate traditional and non-traditional energy streams along with FACTS devices.

• 2. Spotlight on defining appropriate PDFs that encapsulate the inherent randomness of wind energy units.

• 3. The study introduces the Non-dominated Sorting Snow Ablation Optimizer (NSSAO) technique as a promising solution mechanism for OPF challenges involving stochastic non-traditional energy sources.

• 4. The empirical prowess of the MOSAO algorithm is dissected through rigorous comparative studies.

For readers delving deeper into this research, the subsequent sections are organized as follows: Sect. 2 elucidates the mathematical underpinnings of coal-powered energy, wind energy and FACTS devices. Section 3 underscores the targeted objectives warranting optimization. Sections 4 and 5 provides a detailed walkthrough of the multi-objective SAO approach, complete with relevant examples. The subsequent Sect. 6 offers a comprehensive analysis, juxtaposing numerical outcomes, while Sect. 7 rounds off the study with final observations and insights.

2 Mathematical Representations

The research detailed in this article offers an evolved take on the standard IEEE 30-bus testing framework. By seamlessly integrating wind turbines and FACTS devices into the system, we've presented an updated model, thoroughly detailed in Table 1. Figure 1 offers a visual representation of the apparatus used throughout the analysis. In this illustration, the FACTS devices, complete with their ratings, are demarcated with dotted lines, signifying their optimal positions after the refinement process. Following this, the subsequent section delves deeper into the economic dimensions, elaborating on the expenditure patterns linked with traditional coal-driven production units in comparison to those leveraging non-conventional energy sources [26, 34, 35].

Table 1 Test system key features for analysis [26]

Fig. 1

Adapted IEEE 30-bus scheme with wind power units and FACTS policies [26]

2.1 Cost of Coal-Based Power Units

The generation cost is computed using a generalized quadratic equation, as represented in Eq. (1).

$${C}_{T0}({P}_{TG})=\sum_{i=1}^{{N}_{TG}}{a}_{i}+{b}_{i}{P}_{TGi}+{c}_{i}{P}_{TGi}^{2} \quad \mathrm{in \; \$}/{\text{h}},$$

(1)

valve point effect included as;

$${C}_{T}\left({P}_{TG}\right)=\sum_{i=1}^{{N}_{TG}}{a}_{i}+{b}_{i}{P}_{TGi}+{c}_{i}{P}_{TGi}^{2}+\left|{d}_{i}\times {\text{sin}}\left({e}_{i}\times \left({P}_{TGi}^{{\text{min}}}-{P}_{TGi}\right)\right)\right| \quad \mathrm{in\; \$}/{\text{h}}.$$

(2)

The pricing constants for coal-based energy, along with emission rate constants under different scenarios, are presented in reference [26].

2.2 Toxic Gas Emanation

Coal-powered plants emit harmful gases into the environment. Thus, the release of these toxic emissions can be quantified in terms of tons per hour.

$$\text{F2} = \text{t,}E=\sum_{i=1}^{{N}_{TG}}\left[({\alpha }_{i}+{\beta }_{i}{P}_{TGi}+{\gamma }_{i}{P}_{TGi}^{2})\times 0.01+{\omega }_{i}{e}^{({\mu }_{i}{P}_{TGi})}\right] \quad \text{in ton}/{\text{h}}.$$

(3)

Emission constants for harmful gases from coal-powered plants are sourced from reference [25].

2.3 Direct Cost of Stochastic Non-conventional Sources Plants

Integrating non-traditional energy sources into the power grid poses unique challenges due to their unpredictable nature. The task of managing these volatile energy sources falls on the shoulders of the Independent System Operator (ISO). Consequently, private operators need to enter agreements with the ISO or the overarching grid to deliver a predetermined amount of power. The ISO is mandated to ensure the promised power is consistently delivered. If these alternative energy farms falter in maintaining the agreed-upon output, the ISO is held accountable for any power deficiencies. In such cases, they may have to fall back on spinning reserves as a contingency measure. Resorting to spinning reserves introduces additional expenses for the ISO, a situation commonly termed as overestimation of non-traditional energy sources. On the other hand, if these sources generate more energy than anticipated, this surplus might go unutilized, potentially leading to wastage. In such scenarios, the ISO might face penalty charges for this surplus. In essence, the financial implications tied to these energy sources can be categorized into three main costs: the cost of the scheduled power, the additional expenditure resulting from overestimation and reliance on spinning reserves and the penalties incurred due to underutilization of generated energy. The direct cost linked to the wind farms is demonstrated with the ${P}_{ws}$ scheduled power from the same sources as;

$${C}_{w}\left({P}_{ws}\right)={g}_{w}{P}_{ws}.$$

(4)

2.4 Indeterminate Non-conventional Sources of Wind Power Cost

The unpredictable behavior of wind can result in wind farms occasionally delivering less energy than anticipated. When there's an upswing in demand, reliance on spinning reserves becomes essential to meet the pre-determined power quotas. There are instances when the actual energy produced by wind farms falls short of the expected levels. This shortfall in power, when stemming from an uncertain source like wind, is often termed as “overstated power.” In an effort to manage these uncertainties and ensure consistent power delivery to consumers, the grid's Independent System Operator (ISO) resorts to spinning reserves. The financial implication associated with deploying backup generators to cover this overstated power is termed the “reserve cost.” The reserve cost for wind-based units can be defined as follows:

$${C}_{Rw}\left({P}_{ws}-{P}_{wav}\right)={K}_{Rw}\left({P}_{ws}-{P}_{wav}\right)={K}_{Rw}{\int }_{0}^{{P}_{ws}}\left({P}_{ws}-{p}_{w}\right){f}_{w}\left({p}_{w}\right){{\text{d}}p}_{w}.$$

(5)

There are instances when wind farms might produce more energy than what is anticipated, presenting a scenario contrary to overestimation. This situation is commonly referred to as “underestimated power.” Without proper mechanisms in place to regulate output from coal-based units, this surplus energy could go unused. In such cases, the Independent System Operator (ISO) faces penalties for the excess generation. The cost associated with this penalty for wind-based energy is defined as follows:

$${C}_{Pw}\left({P}_{wav}-{P}_{ws}\right)={K}_{Pw}\left({P}_{wav}-{P}_{ws}\right)={K}_{Pw}{\int }_{{P}_{ws}}^{{P}_{wr}}\left({p}_{w}-{P}_{ws}\right){f}_{w}\left({p}_{w}\right){{\text{d}}p}_{w}.$$

(6)

2.5 Uncertainty Models of Stochastic Wind Units

In the revamped IEEE-30 system, thermal generators originally located at buses 5 and 11 have been replaced with wind power units. The Weibull model's constants, specifically the scale ($c$) and shape ($k$), are outlined in Table 2.

Table 2 PDF constants of wind power plants [26]

Using 8000 Monte–Carlo simulations, the Weibull curve and wind frequency distributions were plotted, as depicted in Fig. 2 for the wind facility at bus 5 and Fig. 3 for the one at bus 11. The guidelines highlight the intricacies of wind turbine design and establish the highest turbulence class IA, designed to function optimally at an average annual wind speed of 10 m/s at the turbine's hub height. The scale ($c$) and shape ($k$), parameters of the wind farms are especially significant, given that their peak Weibull average is set near 10. It’s a well-accepted fact that wind speed distribution adheres to the Weibull Probability Density Function (PDF). The probability of a specific wind speed 'v' meters/sec following the Weibull PDF, with scale ($c$) and shape ($k$), parameters, can be determined using the subsequent equation:

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

(7)

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

(8)

gamma function $\Gamma (x)$ is expressed in (9):

Fig. 2

Basic Structure and model of TCSC [26]

Fig. 3

Model of TCPS [26]

$$\Gamma (x)={\int }_{0}^{\infty }{e}^{-t}{t}^{x-1}{\text{d}}t.$$

(9)

2.6 Average power calculation for Wind plants

The wind unit connected at bus number 5 aggregates the outputs from the 25 turbines in the farm. Each of these turbines boasts a power rating of 3 MW. However, the actual output of a wind turbine is influenced by the prevailing wind speed, causing it to fluctuate. To represent the turbine's power output as a function of wind velocity (v), one can employ the subsequent formulas:

$${p}_{w}\left(v\right)=\left\{\begin{array}{ll}0, &{\text{for}} \; v\langle {v}_{{\text{in}}} \; \mathrm{ and}\; v\rangle {v}_{{\text{out}}}\\ {p}_{wr}\left(\frac{v-{v}_{{\text{in}}}}{{v}_{r}-{v}_{{\text{in}}}}\right) &{\text{for}} \;{v}_{{\text{in}}}\le v\le {v}_{r}\\ {p}_{wr} &{\text{for}} \; <v\le {v}_{{\text{out}}}\end{array}\right..$$

(10)

The Enercon E82-E4 design specification is referred to for the 3-MW wind turbine. v_in = 3 m/s, v_r = 16 m/sec and v_out = 25 m/s are the various speeds.

2.7 Wind Power Probabilities Calculation

Wind generation can be unpredictable within specific wind speed intervals. If the wind speed exceeds the cut-out speed or falls below the cut-in speed, no power is generated. However, between these speeds, specifically from the rated to the cut-out speed, the turbine generates a consistent amount of power. The probabilities for these intervals can be characterized as follows:

$${f}_{w}\left({p}_{w}\right)\left\{{p}_{w}=0\right\}=1-{\text{exp}}\left[-{\left(\frac{{v}_{in}}{c}\right)}^{k}\right]+{\text{exp}}\left[-{\left(\frac{{v}_{out}}{c}\right)}^{k}\right],$$

(11)

$${f}_{w}\left({p}_{w}\right)\left\{{p}_{w}={p}_{wr}\right\}={\text{exp}}\left[-{\left(\frac{{v}_{r}}{c}\right)}^{k}\right]-{\text{exp}}\left[-{\left(\frac{{v}_{out}}{c}\right)}^{k}\right].$$

(12)

Within the range of the wind's cut-in speed and its rated speed, wind power generation remains stable. The probability of this consistent zone can be described as follows:

$${f}_{w}\left({p}_{w}\right)=\frac{\beta \left({v}_{r}-{v}_{{\text{in}}}\right)}{{\alpha }^{\beta }*{p}_{wr}}{\left[{v}_{{\text{in}}}+\frac{{p}_{w}}{{p}_{wr}}\left({v}_{r}-{v}_{{\text{in}}}\right)\right]}^{\beta -1}{\text{exp}}\left[-{\left(\frac{{v}_{{\text{in}}}+\frac{{p}_{w}}{{p}_{wr}}\left({v}_{r}-{v}_{in}\right)}{\alpha }\right)}^{\beta }\right].$$

(13)

2.8 Thyristor-Controlled Series Compensator (TCSC) Modelling

Illustrated in Fig. 2 is the fundamental configuration of the TCSC. This design integrates a fixed series capacitor (${X}_{{\text{C}}}$) with a thyristor-controlled reactor (${X}_{{\text{L}}}$). To enable the TCSC to act as an adjustable capacitive reactance, it’s ensured that reactance ${X}_{{\text{C}}}$ is less than ${X}_{{\text{L}}}$. Adjustments to the thyristor's firing angle (α) modify the inductive reactance. When the inductive reactance reaches its peak, the minimum corresponding capacitive reactance emerges (signifying an open circuit for the inductive branch). Thus, the cumulative reactance of the TCSC, encompassing a fixed capacitive reactance ${X}_{{\text{C}}}$ and a modifiable inductive reactance ${X}_{{\text{L}}}$ based on the angle (α), can be expressed as follows: ${X}_{{\text{TCSC}}}(\alpha )$ can be written as

$${X}_{{\text{TCSC}}}\left(\alpha \right)=\frac{{X}_{{\text{C}}}{X}_{{\text{L}}}\left(\alpha \right)}{{X}_{{\text{L}}}\left(\alpha \right)-{X}_{{\text{C}}}}=-j{X}_{{\text{C}}}.$$

(14)

Figure 2 illustrates the static model of the TCSC, positioned between buses m and n. Once the TCSC has been incorporated, functioning in its variable capacitive reactance mode, the modified reactance (${X}_{{\text{eq}}}$) of the transmission line can be defined as:

$${X}_{{\text{eq}}}={X}_{mn}-{X}_{{\text{TCSC}}}=\left(1-\tau \right){X}_{mn},$$

(15)

where,

$$\tau =\frac{{X}_{{\text{TCSC}}}}{{X}_{mn}}.$$

(16)

The power flow equations of the line incorporating TCSC are written as;

$$\begin{aligned} {P}_{mn} & ={V}_{m}^{2}{g}_{mn}-{V}_{m}{V}_{n}{g}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n})\\ & \quad -{V}_{m}{V}_{n}{b}_{mn}{\text{sin}}\left({\delta }_{m}-{\delta }_{n}\right),\end{aligned}$$

(17)

$$\begin{aligned} {Q}_{mn} & =-{V}_{m}^{2}{b}_{mn}-{V}_{m}{V}_{n}{g}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n})\\ & \quad +{V}_{m}{V}_{n}{b}_{mn}{\text{cos}}\left({\delta }_{m}-{\delta }_{n}\right),\end{aligned}$$

(18)

$$\begin{aligned} {P}_{nm} & ={V}_{n}^{2}{g}_{mn}-{V}_{m}{V}_{n}{g}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n})\\ & \quad +{V}_{m}{V}_{n}{b}_{mn}{\text{sin}}\left({\delta }_{m}-{\delta }_{n}\right),\end{aligned}$$

(19)

$$\begin{aligned} {Q}_{nm} & =-{V}_{n}^{2}{b}_{mn}+{V}_{m}{V}_{n}{g}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n})\\ & \quad +{V}_{m}{V}_{n}{b}_{mn}{\text{cos}}\left({\delta }_{m}-{\delta }_{n}\right),\end{aligned}$$

(20)

$${\text{Where}}; \; {g}_{mn}=\frac{{r}_{mn}}{{r}_{mn}^{2}+({x}_{mn}-{x}_{c}{)}^{2}},$$

(21)

$${b}_{mn}=-\frac{{x}_{mn}-{x}_{c}}{{r}_{mn}^{2}+({x}_{mn}-{x}_{c}{)}^{2}}.$$

(22)

2.9 Model of Thyristor-Controlled Phase Shifter (TCPS)

Figure 3 showcases the schematic of the TCPS situated along the transmission line linking buses m and n. With the introduction of the phase shift angle ϕ by the TCPS, the equations governing the power flow through this line can be formulated as follows:

$$\begin{aligned} {P}_{mn} & =\frac{{V}_{m}^{2}{g}_{mn}}{{{\text{cos}}}^{2}\phi }-\frac{{V}_{m}{V}_{n}}{{\text{cos}}\phi }\left[{g}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n}+\phi )\right.\\ &\quad \left.+{b}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n}+\phi )\right],\end{aligned}$$

(23)

$$\begin{aligned} {Q}_{mn} & =-\frac{{V}_{m}^{2}{b}_{mn}}{{{\text{cos}}}^{2}\phi }-\frac{{V}_{m}{V}_{n}}{{\text{cos}}\phi }[{g}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n}+\phi )\\ &\quad -{b}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n}+\phi )],\end{aligned}$$

(24)

$$\begin{aligned} {P}_{nm} & ={V}_{n}^{2}{g}_{mn}-\frac{{V}_{m}{V}_{n}}{{\text{cos}}\phi }[{g}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n}+\phi )\\ & \quad -{b}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n}+\varphi )],\end{aligned}$$

(25)

$$\begin{aligned} {Q}_{nm} & =-{V}_{n}^{2}{b}_{mn}+\frac{{V}_{m}{V}_{n}}{{\text{cos}}\phi }[{g}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n}+\phi )\\ & \quad +{b}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n}+\varphi )].\end{aligned}$$

(26)

The inserted actual and reactive power of TCPS at bus $m$ and $n$ is

$$\begin{aligned} {P}_{ms} & =-{g}_{mn}{V}_{m}^{2}{{\text{tan}}}^{2}\phi -{V}_{m}{V}_{n}{\text{tan}}\phi [{g}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n})\\ &\quad -{b}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n})],\end{aligned}$$

(27)

$$\begin{aligned} {Q}_{ms} & ={b}_{mn}{V}_{m}^{2}{{\text{tan}}}^{2}\phi +{V}_{m}{V}_{n}{\text{tan}}\phi [{g}_{mn}{\text{cos}}({\delta }_{m}-{\delta }_{n})\\ & \quad +{b}_{mn}{\text{sin}}({\delta }_{m}-{\delta }_{n})],\end{aligned}$$

(28)

$${P}_{ns}=-{V}_{m}{V}_{n}{\text{tan}}\phi \left[{g}_{mn}{\text{sin}}\left({\delta }_{m}-{\delta }_{n}\right)+{b}_{mn}{\text{cos}}\left({\delta }_{m}-{\delta }_{n}\right)\right],$$

(29)

$${Q}_{ns}=-{V}_{m}{V}_{n}{\text{tan}}\phi \left[{g}_{mn}{\text{cos}}\left({\delta }_{m}-{\delta }_{n}\right)-{b}_{mn}{\text{sin}}\left({\delta }_{m}-{\delta }_{n}\right)\right].$$

(30)

2.9.1 Model of Static VAR Compensator (SVC)

The basic circuit architecture and the SVC model are depicted in Fig. 4. It is made up of a thyristor-controlled reactor (${X}_{{\text{L}}}=\omega L)$ and a fixed capacitor ${X}_{{\text{C}}}=1/\omega C)$. By changing the thyristor firing angle $(\alpha )$, the reactance can be changed. The equivalent susceptibility is computed as;

Fig. 4

Basic structure and model of SVC [26]

$${B}_{{\text{eq}}}={B}_{{\text{L}}}\left(\alpha \right)+{B}_{{\text{C}}},$$

(31)

$${\text{where}}; \; {B}_{{\text{L}}}\left(\alpha \right)=-\frac{1}{\omega L}\left(1-\frac{2\alpha }{\pi }\right),{B}_{{\text{c}}}=\omega \times C.$$

(32)

The reactive power offered by SVC:

$${Q}_{{\text{SVC}}}=-{V}_{m}^{2}\times {B}_{{\text{SVC}}}.$$

(33)

3 Objectives of Optimization

Optimal power flow (OPF) aims to achieve the ideal distribution of active and reactive power. In this context, several objectives related to wind power flow optimization are integrated. Here are a few instances detailing these objectives:

3.1 Reducing Overall Costs while Using Non-conventional Energy Sources

Objective 1

The primary objective is to reduce the total associated costs. The holistic generation cost is derived by combining the coal-based unit's expense with direct, reserve and penalty fees associated with non-traditional energy sources. Thus, the cumulative cost encompassing both coal-based and wind energy production can be expressed as:

Minimize

$${F1=C}_{{\text{Tot}}}={C}_{T}\left({P}_{TG}\right)+\left[{C}_{w}\left({P}_{ws}\right)+{C}_{Rw}\left({P}_{ws}-{P}_{wav}\right)+{C}_{Pw}\left({P}_{wav}-{P}_{ws}\right)\right].$$

(34)

3.2 Reduction of Voltage Variation with the Use of Non-conventional Energy Sources

Maintaining bus voltage is paramount for both safety and operational excellence. To enhance the voltage profile, voltage deviations at the PQ bus should be confined to a deviation from 1.0 per unit. The goal function can be derived as follows:

Objective 2: Minimize

$$F2=\sum_{i=1}^{{N}_{pq}}\left|{v}_{i}-1.0\right|.$$

(35)

3.3 Minimization of APL Including Non-conventional Energy Sources

The optimization of actual power losses ${P}_{{\text{LOSS}}}$ (MW) maybe calculated by:

Objective 3: Minimize

$${F3=P}_{{\text{LOSS}}}=\sum_{i=1}^{NB}{P}_{Gi}-\sum_{i=1}^{NB}{P}_{Di}.$$

(36)

3.4 Enhancement of VSI Including Non-conventional Energy Sources

The L_max index serves as a critical metric to assess the voltage stability margin of each bus, ensuring that the voltage remains within acceptable bounds during regular operations. This index can be determined using the subsequent equation:

$${L}_{j}=\left|1-\sum_{i=1}^{{N}_{g}}{F}_{ji}\frac{{V}_{i}}{{V}_{j}}\right| \quad \forall \; j=\mathrm{1,2},\dots \dots ,NL,$$

(37)

$${F}_{ji}={-[{Y}_{1}]}^{-1}\left[{Y}_{2}\right],$$

(38)

$$F4=L={\text{max}}\left({L}_{j}\right) \quad \forall \; j=\mathrm{1,2},\dots \dots ,NL.$$

(39)

3.5 Minimization of Entire Gross Cost Including Non-conventional Energy Resources

When comparing Objectives 1 and 3, it’s evident that while one scenario results in a considerably higher generation cost, the other leads to increased loss, can be articulated in this manner:

$${F5=C}_{{\text{Tot}}}+{P}_{{\text{LOSS}}}\times \mathrm{10,000}\times 0.10.$$

(40)

3.6 Equality Constraints

Power flow equations establish equality constraints, indicating that both active and reactive power produced within a system must meet the load requirements and compensate for system losses.

$${P}_{Gi}-{P}_{Di}-{V}_{i}\sum_{j=1}^{NB}{V}_{j}\left[{G}_{ij}{\text{cos}}({\delta }_{ij})+{B}_{ij}{\text{sin}}({\delta }_{ij})\right]=0 \quad \forall \; i\in NB,$$

(41)

$${Q}_{Gi}-{Q}_{Di}-{V}_{i}\sum_{j=1}^{NB}{V}_{j}\left[{G}_{ij}{\text{sin}}({\delta }_{ij})-{B}_{ij}{\text{cos}}({\delta }_{ij})\right]=0 \quad \forall \; i\in NB.$$

(42)

3.7 Inequality Constraints

Inequality constraints represent the operational limits of devices, as well as the safety margins for transmission lines and PQ buses.

Generator bounds:

$${P}_{TGi}^{{\text{min}}}\le {P}_{TGi}\le {P}_{TGi}^{{\text{max}}} \quad \forall \; i\in {N}_{TG},$$

(43)

$${P}_{\text{ws}}^{{\text{min}}}\le {P}_{ws}\le {P}_{\text{ws}}^{{\text{max}}},$$

(44)

$${Q}_{TGi}^{{\text{min}}}\le {Q}_{TGi}\le {Q}_{TGi}^{{\text{max}}} \quad \forall \; i\in {N}_{TG},$$

(45)

$${Q}_{\text{ws}}^{{\text{min}}}\le {Q}_{ws}\le {Q}_{\text{ws}}^{{\text{max}}},$$

(46)

$${V}_{Gi}^{{\text{min}}}\le {V}_{Gi}\le {V}_{Gi}^{{\text{max}}}, \quad i=1\text{,}\dots , NG.$$

(47)

Security bounds:

$${V}_{{L}_{p}}^{{\text{min}}}\le {V}_{{L}_{p}}\le {V}_{{L}_{p}}^{{\text{max}}} , \quad p=1\text{,}\dots , NL,$$

(48)

$${S}_{{l}_{q}}\le {S}_{{l}_{q}}^{{\text{max}}}, \quad q=1, \dots, nl.$$

(49)

FACTS devices bounds:

$${\tau }_{{\text{TCSC}}{\text{m}}}^{{\text{min}}}\le {\tau }_{{\text{TCSC}}m}\le {\tau }_{{\text{TCSC}}{\text{m}}}^{{\text{max}}} \quad \forall \; m\in {N}_{{\text{TCSC}}},$$

(50)

$${\phi }_{{\text{TCPS}}{\text{n}}}^{{\text{min}}}\le {\tau }_{{\text{TCPS}}n}\le {\tau }_{{\text{TCPS}}{\text{n}}}^{{\text{max}}}\forall n\in {N}_{{\text{TCPS}}},$$

(51)

$${Q}_{{\text{SVC}}j}^{{\text{min}}}\le {Q}_{{\text{SVC}}j}\le {Q}_{{\text{SVC}}j}^{{\text{max}}} \quad \forall \; j\in {N}_{{\text{SVC}}}.$$

(52)

The power generation constraints for coal-based and wind energy units are delineated in Eqs. (43) and (44), respectively. Following this, the constraints related to the VAR power generation of these units are represented in Eqs. (45) and (46). $NG$ represents the total number of voltage regulator buses. The voltage boundaries for PV buses are articulated in Eq. (47) and for PQ buses in Eq. (48). $NL$ denotes the total count of PQ buses in the system. Meanwhile, Eq. (49) describes the load-bearing limits of the transmission lines when $NL$ represents the aggregate number of such lines. Lastly, the operational boundaries for devices like TCSC, TCPS and SVC are captured in Eqs. (50)–(52), respectively.

4 Snow Ablation Optimizer Algorithm

In this section, we delve into the foundation of the Snow Ablation Optimizer (SAO), which is rooted in the transformational properties and behaviors of snow, including sublimation and melting processes. Subsequently, we will delineate the mathematical framework of this algorithm before presenting the SAO's pseudo-code and discussing its time complexity.

4.1 Source of Inspiration

The conceptual foundation for the SAO is derived from the intricate processes of snow's melting and sublimation. In subsequent sections, we will elaborate on the initial phase, exploration phase, exploitation phase and the dual-population strategy embedded within the SAO algorithm shown in Fig. 5.

Fig. 5

Schematic diagram of inspiration source

4.2 Initial Phase

Within the SAO, the iterative process commences with the formation of a randomly generated swarm. As delineated in Eq. (53), the entire swarm can generally be depicted as a matrix featuring N rows and Dim columns. Here, “N” signifies the size of the swarm, whereas “Dim” alludes to the dimensionality encompassing the solution space.

$$Z=L+\theta \times \left(U-L\right)={\left[\begin{array}{ccccc}{z}_{\mathrm{1,1}}& {z}_{\mathrm{1,2}}& \cdots & {z}_{1,Dim-1}& {z}_{1,Dim}\\ {z}_{\mathrm{2,1}}& {z}_{\mathrm{2,2}}& \cdots & {z}_{2,Dim-1}& {z}_{2,Dim}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {z}_{N-\mathrm{1,1}}& {z}_{N-\mathrm{1,2}}& \cdots & {z}_{N-1,Dim-1}& {z}_{N-1,Dim}\\ {z}_{N,1}& {z}_{N,2}& \cdots & {z}_{N,Dim-1}& {z}_{N,Dim}\end{array}\right]}_{N\times Dim}.$$

(53)

It is important to note that within this context, “$L$” and “$U$” denote the lower and upper boundaries of the solution space, correspondingly. Additionally, $\theta$ symbolizes a value that is randomly generated within the [0, 1] range.

4.3 Exploration Stage

This segment elaborately discusses the investigation methodology employed in SAO. As snow or its liquid derivative transitions into vapor, the search agents exhibit a heightened decentralization characteristic, brought about by erratic movements. In this research, we employ Brownian motion as a representation of this circumstance. Recognized as a probabilistic phenomenon, Brownian motion finds extensive applications in mimicking the foraging strategies of various fauna, the unpredictable and ceaseless motion of particles and the volatile tendencies observed in stock prices, among other scenarios:

$${f}_{{\text{BM}}}\left(x;\mathrm{0,1}\right)=\frac{1}{\sqrt{2\pi }}\times {\text{exp}}\left(-\frac{{x}^{2}}{2}\right).$$

(54)

Figure 6 showcases the one-dimensional distribution of Brownian motion alongside the trajectories within the two and three-dimensional search realms. As explicated in Fig. 6, the utilization of dynamic and homogeneous step lengths allows the Brownian motion to probe potential zones within the search space effectively, thereby adequately illustrating the dispersal of vapor within the search region. The equation to determine positions throughout the exploration phase is illustrated subsequently:

$$Z_{i} \left( {t + 1} \right) = Elite\left( t \right) + BM_{i} \left( t \right) \otimes \left( {\theta_{1} \times \left( {G\left( t \right) - Z_{i} \left( t \right)} \right) + \left( {1 - \theta_{1} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)} \right),$$

(55)

Fig. 6

Brownian motion

In this context, the term ${Z}_{i}(t)$ signifies the $i$th member during the $t$th iteration, while $B{M}_{i}(t)$ stands for a vector encompassing random values grounded in Gaussian distribution, symbolizing the Brownian motion. The symbol ‘$\otimes$’ is utilized to denote element-wise multiplications and ${\theta }_{1}$ exemplifies a numerical value randomly selected from the range [0, 1]. Moreover, $G(t)$ is a notation that depicts the optimum solution at the current stage and $Elite\left(t\right)$ identifies a randomly chosen individual from a collection of the most proficient members within the swarm. Meanwhile, $\overline{Z}\left( t \right)$ illustrates the geometric center of the entire swarm's positioning. The related mathematical delineations are elucidated as follows:

$$\overline{Z}\left( t \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} Z_{i} \left( t \right),$$

(56)

$${\text{Elite}}\left(t\right)\in \left[G\left(t\right),{Z}_{{\text{second}}}\left(t\right),{Z}_{{\text{third}}}\left(t\right),{Z}_{c}\left(t\right)\right].$$

(57)

In this framework, ${Z}_{{\text{second}}}(t)$ and ${Z}_{{\text{third}}}(t)$ are used to symbolize the individuals who are ranked as second and third best within the present population, correspondingly. The representation ${Z}_{c}(t)$ identifies the central position of members who have fitness evaluations positioned within the upper 50 percentile. For clarity in this research, individuals who are positioned within the top half in terms of fitness values are termed as ‘leaders’. Moreover, the computation of ${Z}_{c}(t)$ is facilitated through the mathematical formula demonstrated in Eq. (58),

$${Z}_{c}(t)=\frac{1}{{N}_{1}}\sum_{i=1}^{{N}_{1}}{Z}_{i}\left(t\right),$$

(58)

wherein ${N}_{1}$ denotes the leader count, equating to half the magnitude of the entire swarm and ${Z}_{i}(t)$ signifies the $i$th highest ranking leader. Consequently, throughout each iteration, the individual designated as $Elite(t)$ is randomly singled out from a group embodying the present.

Figure 7 visually depicts the cross terms -${\theta }_{1}\times (G(t)-{Z}_{i}(t))$ and $\left( {1 - \theta_{1} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)$ in 2-dimensional parametric space. The variable ${\theta }_{1}$ plays a pivotal role in moderating the trajectory towards the present optimum individual as well as the movement directed towards the leaders' centroid position. The amalgamation of the aforementioned two interacting elements primarily serves to mirror the collaborative dynamics existing between individuals.

Fig. 7

Schematic diagram of the cross term in SAO

4.4 Exploitation Stage

In this section, we delineate the focused approach of the SAO algorithm, particularly emphasizing its exploitative attributes. Contrary to employing a highly decentralized method in navigating the solution space, the search agents are prompted to honing in on premium solutions in proximity to the present optimal outcome, akin to the pattern observed when snow transitions into water via melting processes. To illustrate this snow melting transition, we employ one of the most revered snowmelt models, the degree-day approach as highlighted in recent studies. A broad representation of this methodology is showcased below:

$$M=DDF\times \left(T-{T}_{1}\right).$$

(59)

Herein, the term ‘$M$’ stands for the rate of snowmelt, a critical element in mimicking the melting phase during the exploitation phase. The symbol ‘$T$’ indicates the mean daily thermal reading, while ${T}_{1}$ is synonymous with the foundational temperature, typically pegged at 0, leading to the equation:

$$M=DDF\times T.$$

(60)

In this context, $DDF$ signifies the degree-day coefficient, fluctuating between 0.35 and 0.6. Within each iterative process, the mathematical formula to adjust the $DDF$ value is exhibited as follows

$$DDF=0.35+0.25\times \frac{{e}^{\frac{t}{{t}_{{\text{max}}}}}-1}{e-1},$$

(61)

where ${t}_{{\text{max}}}$ demarcating the endpoint criteria. Figure 8 delineates the trajectory of the $DDF$ variable.

Fig. 8

Trends of the $DDF$

$$M=\left(0.35+0.25\times \frac{{e}^{\frac{t}{{t}_{max}}}-1}{e-1}\right)\times T(t),T(t)={e}^{\frac{-t}{{t}_{{\text{max}}}}}.$$

(62)

Proceeding with the SAO strategy, the rate of snowmelt is derived from the succeeding equation:

$$Z_{i} \left( {t + 1} \right) = M \times G\left( t \right) + BM_{i} \left( t \right) \otimes \left( {\theta_{2} \times \left( {G\left( t \right) - Z_{i} \left( t \right)} \right) + \left( {1 - \theta_{2} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)} \right)$$

(63)

As we advance to the exploitation phase within the SAO framework, we introduce the equation responsible for updating positions, outlined below. Here, ‘$M$’ connotes the rate of snowmelt and ${\theta }_{2}$ represents a random variable selected from the interval [− 1, 1], serving as a conduit for interaction amongst the entities. In this stage, with the help of the cross terms—${\theta }_{2}\times (G(t)-{Z}_{i}(t))$ and, $\left( {1 - \theta_{2} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)$. During this juncture, the participants are predisposed to scrutinize potent sectors, drawing from insights acquired from the currently superior search agent and the aggregate central positioning within the swarm.

4.5 Dual-Population Mechanism

Achieving a balanced approach between exploitation and exploration is crucial in the realm of metaheuristic algorithms. As highlighted in Sect. 4.1, a portion of water that originates from snow melting can transition into steam, initiating an exploration phase. This implies that, as time progresses, the propensity of entities to undertake sporadic movements characterized by a high degree of decentralization escalates, urging the algorithm to increasingly venture into various sectors of the solution space. In our research, a two-fold population framework has been crafted to mirror this scenario. The aggregate population and the two derived subgroups are labeled as $P$, ${P}_{a}$ and ${P}_{b}$ respectively. Furthermore, the respective magnitudes of $P$, ${P}_{a}$ and ${P}_{b}$ are represented as $N$, ${N}_{a}$ and${N}_{b}$. ${P}_{a}$ is entrusted with spearheading the exploration endeavours, while ${P}_{b}$ is designated to orchestrate the exploitation initiatives. As the iterations progress, there is a noted diminution in the size of ${P}_{b}$, concurrently witnessing an augmentation in the constituency of ${P}_{a}$. In conclusion, the comprehensive equation for updating positions within the SAO algorithm is illustrated below:

$$Z_{i} \left( {t + 1} \right) = \left\{ {\begin{array}{*{20}l} {Elite\left( t \right) + BM_{i} \left( t \right) \otimes \left( {\theta_{1} \times \left( {G\left( t \right) - Z_{i} \left( t \right)} \right) + \left( {1 - \theta_{1} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)} \right),} \hfill & {i \in index_{a} } \hfill \\ {M \times G\left( t \right) + BM_{i} \left( t \right) \otimes \left( {\theta_{2} \times \left( {G\left( t \right) - Z_{i} \left( t \right)} \right) + \left( {1 - \theta_{2} } \right) \times \left( {\overline{Z}\left( t \right) - Z_{i} \left( t \right)} \right)} \right),} \hfill & {i \in index_{b} } \hfill \\ \end{array} } \right.$$

(64)

As indicated in Eq. (53), the cumulative population is embodied as a position matrix. Therefore, in Eq. (64), ‘$inde{x}_{a}$’ and ‘$inde{x}_{b}$’ symbolize sets of indices which incorporate the row numbers of individuals residing in ${P}_{a}$ and ${P}_{b}$, within the overarching position matrix, respectively. The time complexity of SOA each part is presented in the following: initialization in position updating of individuals $O(N*Dim)$, fitness evaluation $O(N)$, and fitness sorting $O(N*logN)$. Among them, $Dim$ refers to the dimensionality of the solution space. Then the time complexity can be estimated in the following: $O(N*Dim+N*{t}_{max}*(logN+Dim+1))$.The entire process of the SAO optimizer is summarized in Algorithm 1.

Algorithm 1:

Snow ablation optimizer algorithm

5 Multi-objective Snow Ablation Optimizer Algorithm

5.1 Basic Definitions of Multi-objective Optimization

MO optimization is a strategy that permits the simultaneous consideration of multiple objective functions for any given problem. The core components of any MO are detailed below.

$$\left.\begin{array}{ll}Minimize:F\left(\overrightarrow{x}\right)=\left[{f}_{1}\left(\overrightarrow{x}\right),{f}_{2}\left(\overrightarrow{x}\right),\dots ,{f}_{o}\left(\overrightarrow{x}\right)\right] & \\ {\text{Subjected to}}: & \\ {y}_{i}\left(\overrightarrow{x}\right)=0, & i=\mathrm{1,2},\dots ,n\\ {z}_{i}\left(\overrightarrow{x}\right)\ge 0,& i=\mathrm{1,2},\dots ,p\\ {lb}_{i}\le {x}_{i}\le {ub}_{i}, & i=\mathrm{1,2},\dots m\end{array}\right\}$$

(65)

In this context, ${y}_{i}\left(\overrightarrow{x}\right)$ represents the ith equality constraint, while ${z}_{i}\left(\overrightarrow{x}\right)$ symbolizes the $i$th inequality constraint. The variables ‘$n$’ and ‘$p$’ respectively indicate the number of equality and inequality constraints. The vector $F\left(\overrightarrow{x}\right)=\left[{f}_{1}\left(\overrightarrow{x}\right),{f}_{2}\left(\overrightarrow{x}\right),\dots ,{f}_{o}\left(\overrightarrow{x}\right)\right]$ signifies the decision or optimization variables. The lower and upper boundaries for these decision variables are denoted by ‘lb’ and ‘ub’. Here, ‘$m$’ designates the total count of decision variables and ‘o’ represents the number of objective functions. Traditional relational operators become insufficient for comparing solutions within the search space for MO problems. Instead, the concept of Pareto optimality emerges as a pivotal tool for such comparisons. The foundational definitions of MO problems can be outlined as follows, emphasizing the role of the Pareto optimality operator in solution assessments.

Definition 1:

Pareto Optimality (PO).

The solution $\overrightarrow{x}\in X$ is called PO if and only if:

$$\nexists\, \overrightarrow{y}\in X | F\left(\overrightarrow{y}\right)\prec F\left(\overrightarrow{x}\right)$$

(66)

Definition 2:

Pareto Dominance (PD).

Let two different vectors be represented as $\overrightarrow{x}=({x}_{1},{x}_{2},\dots ,{x}_{j})$ and $\overrightarrow{y}=({y}_{1},{y}_{2},\dots ,{y}_{j})$. The vector $\overrightarrow{y}$ is said to dominate the vector $\overrightarrow{x}$ (symbolized as $\overrightarrow{y}\prec \overrightarrow{x}$), if and only if:

$$\forall i\in \left\{\mathrm{1,2},\dots ,j\right\}: {f}_{i}\left(\overrightarrow{y}\right)\le {f}_{i}\left(\overrightarrow{x}\right)\wedge \exists i\in \left\{\mathrm{1,2},\dots ,j\right\}:{f}_{i}\left(\overrightarrow{y}\right)<{f}_{i}\left(\overrightarrow{x}\right)$$

(67)

Definition 3:

Pareto Optimal Set (POS).

The collection of all POS is termed the Pareto set and it can be articulated in the subsequent manner.

$${P}_{s}=\{x,y\in X | \exists F\left(\overrightarrow{y}\right)\succ F(\overrightarrow{x})\}$$

(68)

Definition 4:

Pareto Optimal Front (POF).

The POF consists of a set of optimal solutions found within the POS, as indicated in Eq. (69).

$${P}_{f}=\left\{F\left(\overrightarrow{x}\right)|\overrightarrow{x}\in {P}_{s}\right\}$$

(69)

The Pareto optimal frontier characterizes the relationship between the objective and search domains are illustrated in Fig. 9.

Fig. 9

Objective space and search space in MO-Problem

5.2 Multi-objective Snow Ablation Optimizer (MOSAO)

The advanced MOSAO optimization algorithm employs a combination of the crowding distance (CD) technique and the elitist non-dominated sorting (NDS) approach. The steps involved in the NDS process include:

• 1. Identification of non-dominated solutions.

• 2. Implementation of the NDS strategy.

• 3. Calculation of non-dominated ranking (NDR) for all non-dominated candidates.

The NDR process operates between two layers or fronts. Solutions in the primary front are assigned a rank of “0” because no other solutions dominate them. However, at least one solution from this primary front dominates solutions in the subsequent front. The rank assigned to a solution corresponds to the count of solutions that overshadow it. The CD mechanism ensures that the solutions generated span a diverse range. The CD process is elaborated upon in the following section.

$$C{D}_{j}^{i} =\frac{{fobj}_{j}^{i+1}-{fobj}_{j}^{i-1}}{{fobj}_{j}^{{\text{max}}}-{fobj}_{j}^{{\text{min}}}}$$

(70)

where ${fobj}_{j}^{{\text{max}}}$ and ${fobj}_{j}^{{\text{min}}}$ are the max and min values of jth objective function. The diagrammatic illustration of an NDS-based approach is illustrated in Fig. 10.

Fig. 10

The procedure of the NDS approach

The pseudocode for the MOSAO algorithm is illustrated in Algorithm 2. The process initiates by defining essential parameters such as the population size (Np), stopping conditions, the highest number of generational cycles and the ultimate iteration count (Maxit). Subsequently, the objective space vector F for Po is determined by evaluating every objective function using a spontaneously generated initial population Po within the viable search domain S. In the next step, Po undergoes a refinement process via the elitist-based CD and NDS. Subsequently, Po merges with a newly formed population, Pj, to produce a combined population Pi. This combined set Pi is organized based on the information from CD, NDR and elitist non-dominance in Algorithm 2. From this set, the top Np solutions are selected to form the next generation of parent population. This cycle continues until the stipulated stopping conditions are achieved. The time complexity of MOSAO is $O(M*N*N)$. A visual representation of the MOSAO process can be found in Fig. 11.

Fig. 11

Flowchart of MOSAO algorithm

Algorithm 2:

Pseudocode of Multi-objective Snow ablation optimizer (MOSAO)

5.3 Constraint Handling Approach

In real-world engineering design challenges, many problems are MO and possess complex nonlinear constraints. Effectively addressing constrained MOPs necessitates ensuring all constraints remain within their specified limits. While there are multiple constraint handling methods found in literature, the MOSAO algorithm employs a static penalty approach. This method effectively converts a constrained problem into an unconstrained one by imposing a hefty penalty to the corresponding objective function when a constraint is violated. Below, we delve into the intricacies of the static penalty mechanism.

$${f}_{j}\left(X\right)={f}_{j}\left(X\right)+\sum_{i=1}^{p}{P}_{i}\;{\text{max}}\left\{{g}_{i}\left(X\right)\text{,}0\right\}+\sum_{i=p}^{NC}{P}_{i}\;{\text{max}}\left\{\left|{h}_{i}\left(X\right)\right|-\delta \text{,}0\right\},$$

(71)

where ${f}_{j}\left(X\right) \; ,\,j=1, 2\dots n$ isobjective function, $X \,is$ design variables, ${g}_{i}\left(X\right)\, is$ inequality constraints, ${h}_{i}\left(X\right) \,is$ equality constraints and $\delta$ is tolerance.

5.4 Fuzzy Approach for the Multi-objective Problem

The fuzzy membership methodology can be employed in multi-objective scenarios to pinpoint the most balanced solution among all optimal outcomes. By monitoring the lowest ${f}_{i}^{{\text{min}}}$ and highest ${f}_{i}^{{\text{max}}}$ values for each objective, the fuzzy membership function ${\mu }_{{f}_{i}}$ offers insight. Consequently, the membership function for the $i$th objective can be defined as:

$$\mu _{{f_{i} }} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {f_{i} \le f_{i}^{{{\text{min}}}} } \hfill \\ {\frac{{f_{i}^{{{\text{max}}}} - f_{i} }}{{f_{i}^{{{\text{max}}}} - f_{i}^{{{\text{min}}}} }}} \hfill & {{\text{ }}f_{i}^{{{\text{min}}}} < f_{i} < f_{i}^{{{\text{max}}}} } \hfill \\ 0 \hfill & {f_{i} \ge f_{i}^{{{\text{max}}}} } \hfill \\ \end{array} } \right..$$

(72)

Membership function values range between 0 and 1, indicating the extent to which the function ${f}_{i}$ is satisfied. Subsequently, the decision-making function, denoted as ${\mu }^{k}$ can be computed using the given equation:

$${\mu }^{k}=\frac{\sum_{i=1}^{N}{\mu }_{{f}_{i}}^{k}}{\sum_{k=1}^{M}\sum_{i=1}^{N}{\mu }_{{f}_{i}}^{k}}.$$

(73)

The learning procedure of the MOSAO algorithm for addressing the MO-OPF with wind and FACTS device problem.

Step 1: initiate by setting up line data, real and reactive power load levels, and generator specifics. Establish initial values for real power output at PV buses, reactive power generation of shunt capacitors, voltage magnitudes at generator buses, wind data, FACTS device data, and the tap positions of adjustable transformers. Set the number of control variables, population size, generation count, and the control variables' lower and upper limits.

Step 2: formulate an initial random population within the predefined control variable bounds, based on the control parameters.

Step 3: conduct power flow calculations for each candidate solution and evaluate corresponding objective function values. Adjust any updated individual violating design variable constraints to align with those constraints.

Step 4: identify non-dominated (Pareto-optimal) solutions within the initial population using the Pareto optimal method, storing these in the external elitist repository.

Step 5: update the fittest vectors in the best population using the SAO algorithm’s concept.

Step 6: in revised candidate in the best population, recalculating objective function values and incorporating penalties for constraint violations as per Eq. (71).

Step 7: compare each candidate in the previous and modified best population, and keep the better individual.

Step 8: search for non-dominated solutions in the updated best population, amending the external repository and remove the least significant solution based on crowding distance (Sect. 5.2).

Step 9: conclude and illustrate the Pareto optimal front if the end criterion is met; otherwise, revert to Step 5.

Step 10: organize the final Pareto front’s non-dominated solutions in descending order of fuzzy function values (Sect. 5.4). Select the solution with the highest membership value as the optimal solution.

6 Simulation Results, Analysis and Comparative Study

In this segment, we delve into the performance and results yielded by the Multi-objective Snow Ablation Optimizer (MOSAO) algorithm, which has been employed to optimize the Optimal Power Flow (OPF) considering the integration of non-traditional energy sources and FACTS devices, while managing various control variables.

Algorithm Initialization:

Before the MOSAO algorithm is deployed, several basic parameters are established:

• Population Size: determines the number of potential solutions considered in each iteration.

• Archive Size: denotes the capacity to store previous optimal or near-optimal solutions.

• Maximum Iterations: limits the number of times the algorithm will cycle through its optimization routine.

• Boundary Conditions: sets the permissible range for OPF problems to ensure feasible and meaningful solutions. 

With these parameters in place, the MOSAO algorithm embarks on its quest to locate the most advantageous balance amidst multiple objective functions. It commences by establishing an initial stance, from which it evaluates respective objective function values.

• Algorithm Application and Verification:

The prowess of the MOSAO algorithm, once postulated, is put to the test upon OPF scenarios which harmonize non-conventional energy vectors with FACTS devices. As a precursor to its primary application, the algorithm's mettle is gauged against eight standard multi-objective problems that are devoid of constraints.

• Simulation Details

The evaluation harnesses the computational strength of MATLAB, executed on a computer setup boasting 4 GB RAM and clocking a speed of 3.20 GHz. Each of these benchmark tests is executed ten times, ensuring a thorough assessment and mitigating the influence of outliers or anomalous results.

• Algorithm Parameters

During these tests, the MOSAO operates under predefined parameters:

• Population size capped at 30.

• A limit of 500 iterations.

• An archive storage capacity of 30.

• Performance Metrics

The efficacy of the MOSAO is gauged through a suite of performance indicators:

• Generational Distance (GD): measures the proximity between the obtained solutions and the true Pareto front.

• Inversion Generational Distance (IGD): an inversion of GD, it calculates the mean distance from the true PSs to the nearest point in the obtained solution set.

• Spacing Metrics (SP): assesses the spacing between solutions in the Pareto front.

• Diversity Metrics (DM): provides insights into the variance and spread amongst the derived solutions.

• Spread Metrics (SD): evaluates the breadth and dispersal of the solutions across the Pareto front.

With these metrics, the performance of the MOSAO algorithm is thoroughly scrutinized, offering a comprehensive view of its capabilities in addressing intricate multi-objective optimization problems.

6.1 MOSAO Results for Test Benchmark Problems

Before venturing into complex, real-world challenges, the efficacy of the MOSAO is first benchmarked against a set of standard unconstrained test functions, as documented in [21]. We focus our attention on eight renowned benchmark functions—ZDT1 through ZDT6, KURSAVE, SCHAFFER-1 and SCHAFFER-2 to understand the nuances of MOSAO's performance. Parameter Setting: the underpinning settings of any algorithm, often referred to as control parameters, are instrumental in the quality and precision of the optimization outcome. Upon comprehensive evaluation, the optimal parameter combination was set: 30 for population size, 500 for maximum iterations and 30 for archive size, specific to these benchmark functions. Performance Indicators: the MOSAO algorithm was then put under the lens using a suite of performance metrics. These encompassed—Generational Distance (GD): a measure of solution closeness to the ideal; Inversion Generational Distance (IGD): a metric that quantifies the average distance from the true Pareto solutions to the closest derived solution; Spacing Metrics (SP): an indicator of even distribution of solutions; Diversity Metrics (DM): captures the variety amongst derived solutions and Spread Metrics (SD): assesses the overall spread across the Pareto front. Detailed results from Tables 3, 4, 5, 6, and 7 are telling of MOSAO's commendable capabilities. The algorithm exhibited stellar performance across the board, excelling in metrics such as GD, IGD, SP, DM and SD, which span both convergence speed and solution precision. Visual Insights: The final solutions (or archive solutions) for all eight benchmark tests have been vividly portrayed in Fig. 12. A clear observation from these visuals is MOSAO's adeptness in closely approximating the Pareto Front (PF). A side-by-side comparison with other methods underscores the superior performance of the MOSAO in these scenarios. Concluding Observations: given its impressive display on the benchmark tests, it's deduced that the MOSAO algorithm stands out as an ideal candidate for more intricate challenges, especially the stochastic Optimal Power Flow (OPF) problem interlaced with three FACTS devices and wind energy configurations.

Table 3 Results of GD metric on test functions

Table 4 Results of IGD metric on test functions

Table 5 Results of spacing metric on test functions

Table 6 Results of diversity metric on test functions

Table 7 Results of spread metric on test functions

Fig. 12

Optimal PF using the MOSAO optimizer achieved for unconstrained test functions. a ZDT-1, b ZDT-2, c ZDT-3, d ZDT-4, e ZDT-6, f KURSAVE, g SCHAFFER-1 and h SCHAFFER-2

6.2 MO-OPF Problem with Wind Power Plants and Three FACTS Devices

In this research, the focus is on harnessing the Snow Ablation Optimizer (SAO) algorithm to address the complexities of the stochastic OPF problem, particularly when incorporating wind power facilities and three FACTS devices. This investigation goes hand in hand with a comparative study. We’ve also assessed the performance of some of the recent algorithmic innovations such as: Sine–Cosine Algorithm (SCA): a computational technique inspired by trigonometric functions that optimize numerical problems. Dragonfly Algorithm (DA): drawing inspiration from the static swarming behavior of dragonflies. Ant Lion Optimizer (ALO): this is a bio-inspired optimization algorithm. Ion Motion Algorithms (IMO): modeled after the movement of ions, this is yet another fascinating addition to the realm of optimization algorithms. Our experimentation platform employs a modified version of the IEEE-30 bus system, intricately integrated with both wind power generators and FACTS devices. The nuances of this custom IEEE-30 bus setup, including its primary specifications and configurations, are detailed in Table 1. For the sake of clarity and focused evaluation, our research problem is divided into two distinct scenarios: Scenario 1: it's a singular objective OPF problem, contemplating only the integration of wind power plants and FACTS devices. Scenario 2: a more comprehensive approach, this looks into a MO-OPF, factoring in both the wind power facilities and FACTS devices. Table 8 seems to summarize various case studies for an adapted IEEE-30 bus test system, with each case focusing on different objectives. While all cases involve FACTS devices, wind farms, and coal-based plants, each case is defined by its unique objective used for single objective SAO (case 1 to case 6) or combination of objectives for MOSAO (case 7 to case 11), whether it’s minimizing costs, reducing emissions, improving stability, or a combination of these.

Table 8 Summary of case studies for adapted IEEE-30 bus test system

In Table 8, we see a comprehensive array of eleven distinct test scenarios prepared for in-depth evaluation. This segment presents the findings from these case studies, achieved through the application of various metaheuristic approaches. The initial six cases emphasize singular objective optimization. In contrast, the subsequent five extend their purview to multi-objective optimization challenges. These latter tests notably encompass both unconventional energy sources and the strategic sizing and placement of FACTS devices. We've anchored the search agent value at 30, ensuring each algorithm undergoes a thorough 500-iteration assessment. For a more granular breakdown and exploration of the methodologies involved, readers are encouraged to consult the primary research document.

6.2.1 Scenario-1 (Single Objective OPF with Wind Power Plants and FACTS Devices)

Utilizing methods such as SAO, MVO [36], ALO [32], SCA [37] and IMO [31], each objective detailed in the mathematical framework is addressed as individual optimization challenges. All the constraints of control variables, together with the strategic location cases 1 through 6 are organized in Tables 9, 10, and 11. These tables highlight the optimal minimum values achieved using five contemporary techniques. Notably, when utilizing the SAO method, the overall fuel cost, which encompasses both unconventional power sources and the optimal arrangement of FACTS devices, amounts to $806.999/h. This result is superior when juxtaposed with outcomes from other algorithms, as evidenced in Table 9. When comparing the reduction in Total Fuel Cost (TFC) with methods such as SCA, ALO, DA and IMO, savings of $10.386/h, $1.658/hr, $15.305/h and $8.899/h, respectively are observed. This underscores the preeminence of the SAO technique relative to other benchmarked metaheuristic algorithms. The convergence patterns for TFC minimization are graphically represented in Fig. 13.

Table 9 Case 1 and case 2: individual objective simulation outcomes

Table 10 Case 3 and case 4: individual objective simulation outcomes

Table 11 Case 5 and case 6: individual objective simulation outcomes

Fig. 13

Convergence curve of TFC

Likewise, convergence behaviors for Active Power Loss (APL), VD and Voltage Stability Index (VSI) are delineated in Figs. 14, 15 and 16, respectively. Figure 17 offers a comparative analysis of fuel cost reductions across different algorithms. In the second case study, the SAO method achieves a pollutant gas emission rate of 0.138 tons/hr. In the third case, the APL for various transmission lines under the SAO method stands at 1.445 MW. This demonstrates significant reductions in APL when compared to a myriad of other methods, with reductions ranging from 0.287 to 1.738 MW. Voltage divergence from the 1.0 per unit is a pivotal metric for gauging grid reliability. Hence, in the fourth case study, the SAO algorithm yields the smallest voltage deviation, pegged at 0.114 p.u., outclassing the other four optimization strategies, the fifth case study reveals the SAO method attaining a low of 0.092. In the sixth scenario, the cumulative gross fuel cost as determined by the SAO technique is pinpointed at $1099.401/hr.

Fig. 14

Convergence curve of APL

Fig. 15

Convergence curve of VD

Fig. 16

Convergence curve of VSI

Fig. 17

Comparison of the TFC% reduction

6.2.2 Scenario 2 (MO-OPF with Non-conventional Sources)

In this analysis, multiple objectives, specifically two and three, are concurrently optimized using the Multi-objective Snow Ablation Optimizer (MOSAO) method. MO employs the NDS approach to concurrently derive solutions for various objectives. For the IEEE 30-bus system, a set of 30 non-dominant solutions are preserved.

Scenarios corresponding to Cases 7 through 9 are treated as bi-objective optimization situations. Cases 10 and 11, meanwhile, fall under the category of tri-objective optimization tasks. Leveraging the fuzzy decision-making process, the most suitable equilibrium solution is identified from all the Pareto sets. Optimal solutions derived from the proposed MOSAO and other highlighted metaheuristic methods for Cases 7–11 are emphasized in bold within Tables 12, 13, 14, and 15.

Table 12 Analysis of simulation outcomes for individual objectives

Table 13 MO results for case 7 and case 8

Table 14 MO results for case 9 and case 10

Table 15 MO results for case 11

Figure 18 showcases the optimal PFs of TFC and pollution minimization for Case 7 across different metaheuristic strategies. In a manner akin to Case 7, Fig. 19 illustrates the PF for Case 8, which represents the dual-objective optimization of Active Power Loss (APL) and TFC. The PF for Case 9, focusing on TFC in conjunction with carbon tax and VD minimization, is visualized in Fig. 20.

Fig. 18

PF of TFC and emission

Fig. 19

PF of TFC and APL

Fig. 20

PF of TFC and voltage stability index

Moving on to tri-objective scenarios, Fig. 21 portrays the three-fold objective PFs for Case 10, pinpointing minimizations in APL, TFC and toxic emissions. Similarly, Fig. 22 exhibits the PF for Case 11, highlighting the concurrent minimization efforts on voltage fluctuation, APL and TFC across various algorithmic solutions.

Fig. 21

PF of TFC, emission and APL with MO optimizers

Fig. 22

PF of TFC, APL and VSI with MO optimizers

Drawing upon the data consolidated in the tables, it becomes evident that the MOSAO approach emerges as a leading technique. It adeptly addresses the MOOPF challenges, seamlessly integrating wind power generation and judicious FACTS device placements.

7 Conclusion

The presented methodology focuses on the precise placement and sizing of FACTS (Flexible Alternating Current Transmission Systems) devices. It explores single and multi-objective Optimal Power Flow (OPF) considerations, encompassing both coal-based and wind-driven power plants. Distinct probability density functions are employed to capture the variability inherent in unconventional energy sources. The integration strategy for each unit is meticulously detailed.

In the domain of non-traditional energy sources and FACTS devices, individual objectives such as fuel cost, gas emissions, voltage stability, active loss, and Voltage Stability Index (VSI) are methodically addressed. This research also extends into the multi-objective variant of the OPF challenge, reflecting the contemporary dynamics of the power grid. The derived results are compared with those from a newly established optimization technique.

The findings clearly demonstrate that the introduced Multi-objective Snow Ablation Optimizer (MOSAO) exhibits superior convergence characteristics, distinguishing itself from contemporary algorithms. It consistently delivers high-quality, actionable solutions across a variety of OPF scenarios. Collectively, these outcomes highlight the unique advantage of the proposed approach in tackling OPF problems, whether they are single or multi-faceted in nature.

The MOSAO, particularly when combined with a non-dominated sorting mechanism, proves to be highly effective in managing a range of OPF challenges. Despite the effective performance of MOSAO in addressing various challenges, opportunities for enhancement in other domains still exist. As an optimizer, MOSAO is subject to the limitations outlined by the “No Free Lunch” theorem (NFL); moreover, its nature as a stochastic optimizer introduces an inherent element of randomness. This characteristic occasionally leads to stagnation in convergence, particularly when tackling complex optimization issues. The MOSAO source code is available at: https://github.com/kanak02/MOSAO.