Multi-Objective Optimal Power Flow Analysis Incorporating Renewable Energy Sources and FACTS Devices Using Non-Dominated Sorting Kepler Optimization Algorithm

Abstract

In the rapidly evolving landscape of electrical power systems, optimal power flow (OPF) has become a key factor for efficient energy management, especially with the expanding integration of renewable energy sources (RESs) and Flexible AC Transmission System (FACTS) devices. These elements introduce significant challenges in managing OPF in power grids. Their inherent variability and complexity demand advanced optimization methods to determine the optimal settings that maintain efficient and stable power system operation. This paper introduces a multi-objective version of the Kepler optimization algorithm (KOA) based on the non-dominated sorting (NS) principle referred to as NSKOA to deal with the optimal power flow (OPF) optimization in the IEEE 57-bus power system. The methodology incorporates RES integration alongside multiple types of FACTS devices. The model offers flexibility in determining the size and optimal location of the static var compensator (SVC) and thyristor-controlled series capacitor (TCSC), considering the associated investment costs. Further enhancements were observed when combining the integration of FACTS devices and RESs to the network, achieving a reduction of 6.49% of power production cost and 1.31% from the total cost when considering their investment cost. Moreover, there is a reduction of 9.05% in real power losses (RPLs) and 69.5% in voltage deviations (TVD), while enhancing the voltage stability index (VSI) by approximately 26.80%. In addition to network performance improvement, emissions are reduced by 22.76%. Through extensive simulations and comparative analyses, the findings illustrate that the proposed approach effectively enhances system performance across a variety of operational conditions. The results underscore the significance of employing advanced techniques in modern power systems enhance overall grid resilience and stability.

1. Introduction

In the rapidly evolving landscape of electrical power systems, the concept of optimal power flow has emerged as a cornerstone for efficient energy management. Optimal Power flow is a fundamental process in electrical power systems, especially when the operator system (OS) strives to find the most efficient and economical operating conditions while satisfying a variety of system constraints. It involves the determination of power generation levels, voltage levels at different nodes, and the flow of electricity through transmission lines in a way that minimizes operational costs, such as fuel consumption, or maximizes system efficiency, while ensuring network security and reliability. The challenges of OPF optimization are multifaceted and stem from the complexity of power systems, the non-linear nature of power flow equations, and the need to balance generation and demand in real-time.

Research on the optimal power flow (OPF) that mainly focus on thermal power generators has received extensive attention from researchers. Initially, the resolution of these problems employed traditional optimization methods, among them quadratic programming, interior point methods [1,2], and the Newton method as discussed in reference [3], as well as linear and nonlinear programming as highlighted in [4,5]. Despite their initial utility, these classical approaches encounter notable challenges when applied to large-scale and intricate OPF scenarios characterized by significant nonlinearity and diverse optimization challenges. As a result, there is a tendency for these methods to settle on local optima. This constraint not only limits the use of gradient-based methods but also presents an opportunity to explore the utilization of meta-heuristic- evolutionary algorithms for seeking more efficient solutions.

1.1. Literature Review

Literature reveals that meta-heuristic evolutionary algorithms have been successfully utilized to solve the (OPF) problem. As an illustration, researchers have employed techniques such as the Harris Hawks Optimization (HHO) in [6], the modified pigeon-inspired optimization algorithm (MPIOA) [7], Powell’s pattern-search (PPS) method and the most valuable player algorithm (MVPA) [8], an adaptive multiple team’s perturbation-guiding Jaya (AMTPG-Jaya) technique [9], whale optimization algorithm (WOA) [10], adaptive quasi-oppositional migrated-biogeography-based optimization [11], a new swarm intelligence optimization approach [12], the enhanced equilibrium optimizer (EO) referred to as (EEO) [13], and a space transformational invasive weed optimization (ST-IWO) algorithm [14].

The task of optimizing power flow has grown increasingly complex with the integration of variable renewable energy sources (RESs) like wind, solar, and hydropower. These sources, while environmentally beneficial, introduce significant unpredictability and intermittency, which traditional OPF methods cannot adequately address. To model the stochastic nature of these RESs, advanced probabilistic approaches are employed. For wind power, the Weibull probability density function (PDF) is used to capture wind speed variability, while solar PV output is modeled using a lognormal PDF to represent the fluctuations in irradiance. Hydropower, with its unique variability, is modeled using a Gumbel distribution to account for extreme water flow conditions. These models help to accurately represent the uncertainty and intermittency of RESs, allowing for more robust and reliable OPF solutions in dynamic power grids. Recent studies have focused on various optimization algorithms to address the challenges in OPF with renewable energy integration; as an illustration, Abid et al. [15] examine the (OPF) of a thermal–wind–solar power system using an enhanced (KOA) Kepler optimization algorithm referred to as (EKOA). This research, conducted on a comprehensive practical power system, aims to reduce power production costs, power losses, and toxic emissions while considering the impact of a carbon tax. Duman et al. [16], employed the symbiotic organism search (SOS) approach to address the AC optimal power flow (OPF) problem. Their work explored test scenarios incorporating the stochastic characteristics of wind, solar, and tidal energy systems. Solutions were meticulously evaluated for integrating renewable energy sources into the power networks of the IEEE 30- and 118-bus systems. Various configurations of heat-generating units were strategically positioned within these systems to enhance overall performance. Similarly, in reference [17], the jellyfish search algorithm (JSA) was utilized to improve the financial returns within two IEEE power transmission networks, comprising 30 and 118 buses. The study focused on incorporating wind turbines and electrical market dynamics into the optimization process. T. Samakpong et al. [18] proposed a mutation-based particle swarm optimization (MPSO) technique, which was employed to tackle the OPF problem. This approach incorporated a cost model that accounted for the expenses associated with the uncertainty stemming from the utilization of renewable energy sources (RESs). To address uncertainties, the Monte Carlo method was utilized for simulation and management within the experiment. In [19], Modified Rao-2 (MRao-2) method is suggested to optimize the OPF within a hybrid power system integrating solar, wind, hydro, and biomass sources. The primary aim is to minimize fuel costs across different scenarios for both the IEEE 30-bus and 118-bus systems, while considering RESs under contingency conditions. M. Suleiman et al. [20], applied the barnacles mating optimizer (BMO) to solve the security-constrained optimal power flow (SOPF) problem in a hybrid power system integrating thermal, wind, photovoltaic (PV), and small hydro power sources. It takes into account prohibited operating zones and constraints related to thermal generator limits. In [21], a new multi-objective hybrid evolutionary algorithm (MOEA) is presented, which integrates decomposition techniques with the IWO algorithm. This method is designed to address the OPF problem in transmission networks, especially under conditions of uncertainty. These uncertainties stem from wind energy (WE), PV, and PEV systems, transforming the traditional OPF into a stochastic OPF scenario.

The Flexible AC Transmission Systems (FACTS) devices are crucial for enhancing system performance because of their capability to manipulate system parameters such as transmission-line impedance, voltage magnitudes, and phases, as well as power flow through the lines. Several works have explored the use of FACTS devices in power flow optimization, demonstrating their potential in enhancing system performance. For instance, the SVC provides reactive power compensation, which helps regulate voltage levels on the transmission line. The TCSC can adjust the effective reactance of the transmission line, allowing for better control of power flow. These can help to reduce the power losses, increase the transmission capacity of the line, and mitigate voltage instability. They can improve voltage stability and reduce voltage fluctuations. Several works have explored the use of FACTS devices in power flow optimization, demonstrating their potential in enhancing system performance, namely [22,23,24,25,26,27,28,29]. However, the limitation of these studies is the absence of consideration for RESs. Nusair et al. [30] optimized the rating and sizing of various FACTS devices, such as the SVC, TCPS, and TCSC, in the presence of RESs using recent optimization techniques like the jellyfish search (JS) and marine predators algorithm (MPA). Similarly, Biswas et al. applied the history-based adaptive-differential-evolution (SHADE) technique to solve the optimal power flow (OPF) problem, incorporating both renewable energy resources (RERs) and FACTS devices such as the thyristor-controlled-phase shifter (TCPS), TCSC, and SVC. However, a notable limitation of this research is the exclusive consideration of wind energy resources, with no inclusion of solar PV systems [31]. In reference [32], an enhanced iteration of the hunter–prey optimization (HPO) method has been introduced to amplify its search capabilities in addressing the OPF problem. This advancement encompasses the incorporation of (FACTS) devices alongside the integration of wind power energy. However, only wind generators were considered in the last two references. M. Ebeed et al. [33], utilized a modified version of the Runge-Kutta optimizer (MRUN) to address the stochastic optimal power flow (OPF) problem. Their research concentrated on the optimal integration of wind turbines (WTs) and photovoltaic (PV) systems, alongside a TCSC. However, this study exclusively considered a single type of FACTS device.

1.2. The Novelty and Scope of Work

This study aims to investigate the importance of integrating RESs and FACTS technologies within contemporary power systems. A focal point of this investigation involves determining the optimal location and size of the TCSC and SVC installations. This optimization aims to minimize overall system costs while concurrently reducing real power losses (RPL) and total voltage deviation (TVD), and improving the voltage stability index (VSI). Despite potential additional costs associated with FACTS technology implementation, integrating renewable energy sources offers a promising solution. By decreasing the base cost of energy production and mitigating toxic gas emissions, this integrated approach aligns with environmental and economic sustainability objectives, showcasing a comprehensive strategy to tackle the challenges and seize the opportunities within modern power systems.

• This paper introduces a novel non-dominating sorting KOA referred to as NSKOA, to tackle SOPF problems.
• It addresses the OPF problem by incorporating RESs, namely solar PV, wind, and hydro power systems and FACTs devices such as the SVC and TCSC.
• It optimizes the size and location to maximize the benefits of FACTS devices for the power system.
• The paper utilizes lognormal, Weibull, and Gumbel Probability Density Functions (PDFs) to effectively model and characterize the RES uncertainties within the system.
• A statistical analysis is performed to confirm the effectiveness of the proposed NSKOA and to highlight the advantages gained from integrating RES and FACTS devices.

The remainder of this paper is organized as follows: Section 2 presents a comprehensive review of the mathematical formulation of the OPF problem, Section 3 describes the stochastic models for wind, solar PV, and hydropower generation within the context of RES integration. Section 4 describes the methodology, providing an in-depth explanation of the Kepler optimization algorithm. Section 5 details the simulation setup and discusses the results. Lastly, Section 6 summarizes the conclusions and explores potential avenues for future research.

2. Problem Formulation

The problem presented in this paper aims to solve the SOPF for a power system that incorporates thermal generation as well as stochastic wind and solar PV power generations. The main objective is to determine the optimal settings for control variables in various power system components, while maintaining adherence to all equality and inequality constraints. The mathematical formulation of this problem is as follows:

Minimize

$$OF(d,c)=\left\{O{F}_1(d,c),O{F}_2(d,c),O{F}_3(d,c),\dots \dots ,O{F}_{Nobj}(d,c)\right\}$$

(1)

Subject to

$$g\left(d,c\right)\left\{\begin{array}{l}g\left(d,c\right)=0\\ h\left(d,c\right)\le 0\end{array}\right.$$

(2)

Ensuring both system security and optimal outcomes within an electrical network requires strict adherence to constraints placed on control variables. These limits serve as fundamental safeguards, essential for upholding feasibility, guaranteeing system stability.

-

$OF\left(d,c\right)$ denotes the objective function that needs to be minimized.

-

$g\left(d,c\right)$ represents the collection of equality constraints that must be fulfilled.

-

$c$ represents the vector of decision variables, and d represent the vector of state variables.

2.1. Optimization Problem

2.1.1. Cost of Generation for Thermal Units

The fuel cost associated with thermal generator units can be represented as a smooth and convex quadratic function. This mathematical representation is denoted by Equation (3):

$$C_{TH}\left(P_{TG}\right)= {\displaystyle \sum _{i=1}^{N_{TG}}{a}_i+b_i{P}_{TGi}+c_i{P}_{TGi}^2}$$

(3)

The model described above ignores valve point loading; however, when valve point effects are considered, an additional sinusoidal term is incorporated into the equation to capture the oscillations or fluctuations introduced by the valve points. This modification accounts for the non-linearity in the input–output curve of thermal generators and provides a more accurate representation of the fuel cost function in power system optimization. The valve point-effect can be modeled mathematically using Equation (4) [15]:

$$T{C}_{TH}\left(P_{TG}\right)= {\displaystyle \sum _{i=1}^{N_{TG}}{a}_i+b_i{P}_{TGi}+c_i{P}_{TGi}^2}+\left|d_i\times \sin \left(e_i\left(P_{TG}^{\min }-P_{TG}\right)\right)\right|$$

(4)

$a_i, b_i, c_i, d_i, \mathrm{and} e_i$ are the cost coefficients for the i^th thermal generator. The system comprises $N_{ThG}$ conventional generators, and $P_{ThG}^{,in}$ represents the minimum rated power of these generators.

2.1.2. The Investment Cost of FACTS Modeling

• SVC modeling

The static var compensator (SVC) can exhibit two distinct characteristics: inductive or capacitive. In the former, it absorbs reactive power, while in the latter, it injects reactive power. The SVC is composed of a series capacitor bank that is shunted by a thyristor-controlled reactor, as illustrated in Figure 1.

According to [34], the investment cost of static var compensators (SVCs) varies linearly depending on the reactive power of the SVC to be installed. Therefore, the cost at node I is expressed as follows:

$$C_{SV{C}_i}=0.0003\times {Q_{SV{C}_i}}^2-0.3051\times Q_{SV{C}_i}+127.38$$

(5)

The total investment cost is given as follows:

$$C_{SVC}={\displaystyle \sum _{i=1}^{N_{SVC}}(0.0003\times {Q_{SV{C}_i}}^2-0.3051\times Q_{SV{C}_i}+127.38})$$

(6)

where $Q_{SV{C}_i}$ represent the reactive power generated by i^th SVC while $C_{SVCi}$ represents its associated cost. The total investment cost for all SVC devices is referred to as $C_{SVC}$, and $N_{SVC}$ is the total number of SVC devices.

• TCSC modeling

The thyristor-controlled series compensator (TCSC) is a series compensation device comprising a series capacitor bank shunted by a thyristor-controlled reactor, as shown in Figure 2a,b. The primary concept behind power flow control using the TCSC is to adjust the overall effective series transmission impedance of the lines, either decreasing or increasing it by introducing capacitive or inductive reactive components, respectively. The TCSC is represented as a variable impedance, as illustrated in Figure 2.

The total investment cost of TCSCs is a quadratic function of the reactive power to be installed, expressed as follows [35]:

$$C_{TCSC}={\displaystyle \sum _{i=1}^{N_{TCSC}}((0.0015\times {Q^2}_{TCS{C}_i}-0.7130\times \left|Q_{TCS{C}_i}\right|+153.75})\times 1000\times Q_{TCS{C}_i})$$

(7)

where $Q_{TCS{C}_i}$ represents the reactive power generated by the i-^th TCSC. The total investment cost for the all TCSC devices is referred to as $C_{TCSC}$, and $N_{TCSC}$ is the total number of TCSC devices.

2.1.3. Cost Generation for Renewable Sources

The RES generators, not requiring fuel, have their cost functions evaluated according to specific standards, incorporating operational costs like direct, reserve, and penalty fees.

Direct Cost of RES Generators (DCost)

The DCost encompasses the expenses related to generating RESs, comprising equipment, maintenance, and ongoing operational costs. These costs predominantly entail both the initial investment and the continuous operational endeavors associated with RES production. Specifically, for the wind generator, this is referred to as $D\cos t_{W, j}$, and for the solar generator, this is referred to as $D\cos t_{S, k}$. The direct costs can be mathematically expressed as follows:

$$Dcos{t}_{S, k}\left(S{P}_{Sh}\right)=h_k.S{P}_{Sc,k}$$

(8)

$$Dcos{t}_{W, j}\left(W{P}_{SC,j}\right)=g_j.W{P}_{SC,j}$$

(9)

The mathematical expression for the DCost function of the combined solar–hydro generation plants is as follows:

$$Dcos{t}_{Sh, k}\left(Sh{P}_{Sh}\right)=C(S{P}_{sc}+h{P}_{sc})=h_k.S{P}_{Sc,k}+H_{i}h{P}_{sc,i}$$

(10)

with the $k^{th}$ solar power unit and $J^{th}$ wind power plant and $i^{th}$ hydro power plant, respectively, and denoting the scheduled power from the corresponding RES power plants. $h{P}_{SC,i}$ is the scheduled power from the solar–hydro power plant.

The Evaluation of Cost Uncertainties in RES Generators

In uncertain conditions, two scenarios can occur. When the actual power generated from wind or solar sources is less than the estimated amount, known as power overestimation, operator system necessitates to use the reserve of available resources to ensure supply continuity. The expense associated with activating these reserve units to make up for the overestimated power is referred to as the reserve cost [33]. For wind, solar PV, and solar–hydro generators, the reserve cost is mathematically expressed as follows:

Each constraint of operational constraints is treated as follows:

$$\begin{array}{l}RCos{t}_{W,j}\left(W{P}_{Sc,j}-W{P}_{Av,j}\right)=R{K}_{W,j}\left(W{P}_{Sc,j}-W{P}_{Av,j}\right)\\ = R{K}_{W,j}{\displaystyle {\int }_0^{P_{WG,j}}\left(W{p}_{Sc,j}-W{P}_{,j}\right)} f_W\left(W{P}_{,j}\right)d{p}_{W,j}\end{array}$$

(11)

$$\begin{array}{l}RCos{t}_{RS,k}\left(S{P}_{Sc,k}-S{P}_{Av,k}\right)=R{K}_{S,k}\left(S{P}_{Sc,k}-S{P}_{Av,k}\right)\\ = R{K}_{S,k}\ast f_S\left(S{P}_{Av,k}<S{P}_{Sc,k}\right)\ast \left[S{P}_{Sc,k}-E\left(S{P}_{Av,k}<S{P}_{S,k}\right)\right]\end{array}$$

(12)

$$\begin{array}{l}RCos{t}_{sh}(Sh{P}_{sc}-Sh{P}_{Av})=R{K}_{SH,i}(Sh{P}_{sc}-Sh{P}_{Av})\\ =K_{RSH}\ast f_{SH}(Sh{P}_{Av}<Sh{P}_{sc})\ast \left[P_{SHG}-E(Sh{P}_{Av}<Sh{P}_{sc})\right]\end{array}$$

(13)

where $P{K}_{W,j}$, $P{K}_{S,k}$, and $P{K}_{SH}$ are the penalty cost coefficients associated with the j^th, kth, and $i^{th}$ wind, solar, and solar–hydro power generator, respectively, and $W{P}_{r,j}$ denotes the output power from the corresponding wind units. $f_S\left(S{P}_{Av,k}>S{P}_{Sc,k}\right)$ represents the probability of surplus solar power, indicating actual power surpassing the scheduled power ($S{P}_{Sc,k}$), while $E\left(S{P}_{Av,k}>S{P}_{Sc,k}\right)$ signifies the anticipation of solar power exceeding the scheduled power ($S{P}_{Sc,k}$).

$f_{SH}(Sh{P}_{Av}>Sh{P}_{sc})$ represents the likelihood of energy exceeding the scheduled power, while $E(Sh{P}_{Av}>Sh{P}_{sc})$ denotes the forecast of the combined system power surpassing ($Sh{P}_{sc}$).

$$\begin{array}{l}P{C}_{W,j}\left(W{P}_{Av,j}-W{P}_{S,j}\right)=P{K}_{W,j}\left(W{P}_{Av,j}-W{P}_{S,j}\right)\\ = P{K}_{W,j}{\displaystyle {\int }_{W{P}_{S,j}}^{W{P}_{r,j}}\left(W{P}_{,j}-W{P}_{S,j}\right)} f_W\left(W{P}_{,j}\right)d{p}_{W,j}\end{array}$$

(14)

$$\begin{array}{l}P{C}_{S,k}\left(S{P}_{A,k}-S{P}_{Sc,k}\right)=P{K}_{S,k}\left(S{P}_{Av,k}-S{P}_{Sc,k}\right)\\ = P{K}_{S,k}\ast f_S\left(S{P}_{Av,k}>S{P}_{Sc,k}\right)\ast \left[E\left(S{P}_{Av,k}>S{P}_{Sc,k}\right)-S{P}_{Sc,k}\right]\end{array}$$

(15)

$$\begin{array}{l}PCos{t}_{sh}(Sh{P}_{Av}-Sh{P}_{sc})=P{K}_{SH}(Sh{P}_{Av}-Sh{P}_{sc}) \\ =P{K}_{SH}\ast f_{SH}(Sh{P}_{Av}>Sh{P}_{sc})\ast \left[E(Sh{P}_{Av}>Sh{P}_{sc})-Sh{P}_{sc})\right]\end{array}$$

(16)

2.2. Objective Functions

2.2.1. Minimization of Power Production Cost

The first objective function represents the total cost of energy production, incorporating the presence of RESs and all their relevant cost functions. Mathematically, it can be expressed as follows:

$$\begin{array}{l}\mathrm{T}\mathrm{G}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}= T{C}_{Th}\left(P_{TG}\right)+{\displaystyle \sum _{j=1}^{N_{WG}}\left[C_{W,j}\left(W{P}_{Sc,j}\right)+R{C}_{W,j}\left(W{P}_{G,j}-W{P}_{Av,j}\right)+P{C}_{W,j}\left(W{P}_{Av,j}-W{P}_{Sc,j}\right)\right]}\\ +{\displaystyle \sum _{k=1}^{N_{SG}}\left[C_{S,k}\left(S{P}_{Sc,k}\right)+R{C}_{S,k}\left(S{P}_{Sc,k}-S{P}_{Av,k}\right)+P{C}_{S,k}\left(S{P}_{Av,k}-S{P}_{Sc,k}\right)\right]} \\ +{\displaystyle \sum _{i=1}^{N_{ShG}}\left[C_{Sh,i}(Sh{P}_{sc})+R{C}_{SH}(Sh{P}_{sc}-Sh{P}_{Av})+P{C}_{Sh,i}(Sh{P}_{Av}-Sh{P}_{sc})\right]} \end{array}$$

(17)

where $N_{WG}$, $N_{SG}$, and $N_{ShG}$ represent the number of wind, solar PV, and solar–hydro power generators, respectively, in the grid.

2.2.2. Real Power Losses (RPLs)

In the context of the OPF problem, it is important to consider additional network parameters, including the power loss incurred during system transmission. These parameters play a critical role in assessing the efficiency and stability of the system. The calculation of the (TPL) is expressed through the following equation:

$$RPL={\displaystyle \sum _{q=1}^{nl}{G}_{q(ij)}\left[V_i^2+V_j^2-2V_i{V}_j\cos ({\delta }_i-{\delta }_j)\right]}$$

(18)

where $G_{q(ij)}$ is the conductance of the branch, $nl$ is the number of transmission lines, $V_i$ and $V_j$ are the voltages at bus i and j, respectively, and ${\delta }_{ij}={\delta }_i-{\delta }_j$ is the difference in voltage angles between them.

2.2.3. Total Voltage Deviation (TVD)

Voltage deviation is an indicator used to assess the quality of voltage within a power network. It quantifies the total variation between the voltages at all load buses (PQ buses) and the standard nominal value of one per unit (p.u.). This parameter is determined by summing the absolute differences between the voltage at each load bus and the nominal value. The mathematical expression for calculating voltage deviation is as follows:

$$VD=\left({\displaystyle \sum _{p=1}^{NL}\left|V{L}_p-1\right|}\right)$$

(19)

2.2.4. Voltage Stability Index (VSI)

The importance of monitoring and controlling power networks has grown significantly within the operation of contemporary electrical power systems, particularly with regards to enhancing voltage stability amid the increasing integration of renewable energies. To better understand voltage drops, the operational range of index L has been defined as (0, 1) Consequently, the third objective function, aimed at minimizing the voltage stability index within the transmission branches, can be modeled as follows:

$$\mathrm{Min}\left(\mathrm{VSI}\right)=\min \left(\max \right(L_j\left)\right)$$

where ($L_j$) of the j-th bus is calculated using the following equation:

$$L_j=\left|1-{\displaystyle \sum _{i=1}^{Ngb}{F}_{ij}\times \frac{V_i}{V_j}\angle }\left\{{\theta }_{ij}+({\delta }_i-{\delta }_j\right\}\right|j=1,2\dots \dots N_{lb}$$

(20)

$$F_{ij}=\overline{){\theta }_{ij}}, V_i= \left|V_i\right|\overline{){\theta }_i}, V_j= \left|V_j\right|\overline{){\theta }_j}$$

(21)

$$F_{ij}=-{\left[Y_1\right]}^{-1}\times \left[Y_2\right]$$

(22)

$V_i$ and $V_j$ represent the voltage magnitudes at bus i and j. ${\theta }_{ij}$ signifies the voltage angle difference between bus i and j. $Ngb$ and $N_{lb}$ stand for the number of generator and load buses, respectively.

$Y_1$, $Y_2$, $Y_3$, and $Y_4$ represent the sub-matrices of the system Ybus, obtained through the rearrangement of the generator and load bus parameters as shown in Equation (22).

$$\left[\begin{array}{l}{I}_{gb}\\ I_{lb}\end{array}\right]=\left[\begin{array}{l}{Y}_1 Y_2\\ Y_3 Y_4\end{array}\right]\times \left[\begin{array}{l}{V}_{gb}\\ V_{lb}\end{array}\right]$$

(23)

2.3. Constraints

2.3.1. Equality Constraints

In this study, the specific equality constraints represent the fundamental power balance in the system, ensuring that the total power generation equals the total demand plus transmission losses. This is critical for maintaining a stable and balanced power system. Mathematically, this is expressed as the sum of power generated at all buses being equal to the load demand plus losses. These constraints ensure the continuous supply of power and the reliable operation of the network, which is essential in both traditional and RES-integrated systems.

$$P_{TGi}^{\min }\le P_{TGi} \le P_{TGi}^{\max }, i= 1,2\dots \dots N_{TG}$$

(24)

$$P_{wsj}^{\min }\le P_{wsj} \le P_{wsj}^{\max }, j= 1,2\dots \dots N_{WG}$$

(25)

$$P_{SS, k}^{\min }\le P_{SS,k} \le P_{SS, k}^{\max }, k= 1,2\dots \dots N_{SG}$$

(26)

$$Q_{TGi}^{\min }\le Q_{TGi} \le Q_{TGi}^{\max },$$

(27)

$$Q_{wsj}^{\min }\le Q_{wsj} \le Q_{wsj}^{\max }, i= 1,2\dots \dots N$$

(28)

$$Q_{ss,k}^{\min }\le Q_{ss,k}\le Q_{ss,k}^{\max } k\in N_{SG}$$

(29)

$$Q_{Ci}^{\min }\le Q_{Ci}\le Q_{Ci}^{\max } i\in N_C$$

(30)

$$V_{Gi}^{\min }\le V_{Gi} \le V_{Gi}^{\max },$$

(31)

$$V_{Li}^{\min }\le V_{Li}\le V_{Li}^{\max } i\in NL$$

(32)

$$b_{SVC}^{\min }\le b_{SVC} \le b_{SVC}^{\max }, u \in N_{SVC}$$

(33)

$$X_{TCSC}^{\min }\le X_{TCSC} \le X_{TCSC}^{\max }, W \in N_{TCSC}$$

(34)

2.3.2. Security Constraints

The inequality constraints, on the other hand, include operational limits such as generator output limits, voltage limits at buses, and thermal limits of transmission lines. These constraints ensure that the system operates within safe and efficient boundaries.

$$T_k^{\min }\le T_k\le T_k^{\max } k\in NT$$

(35)

$$S_i^{\min }\le S_i\le S_i^{\max } i\in NL$$

(36)

In this context, max and min are the upper and lower boundaries, $N_{TG}$, $N_{WG}$, $N_{SG}$, and $N_{ShG}$ refer to the number of thermal, wind, solar PV, and solar–hydro generators, and $NL$ is the number of load buses the superiority of feasible solutions technique used in this study to ensure solution feasibility. More details are provided in [36].

3. RES Uncertainty Models

Assessing RES units relies on wind speed and solar radiation probability distribution functions (PDFs). Incorporating these PDFs helps estimate energy output and assess project economics, crucial for designing and implementing efficient renewable energy systems.

The Weibull distribution is a common model for wind speed probability distribution in wind generators. It incorporates scale factor (c) and shape factor (k) to describe wind speed probabilities. Historical wind data analysis determines Weibull distribution parameters and constructs the probability density function (PDF). This assists in evaluating energy production potential and wind generator performance at specific locations.

Wind speed, being a stochastic variable, follows a Weibull probability distribution function (PDF) defined by parameters known as the shape factor (k) and scale factor (c). This relationship can be modeled as follows:

$$f_{\upsilon }\left(S\right)=\left({\displaystyle \frac{k}{c}}\right){\left({\displaystyle \frac{S}{c}}\right)}^{\left(k-1\right)}\times \exp -{\left({\displaystyle \frac{S}{c}}\right)}^k for 0<S<\infty$$

(37)

3.1. Wind Power Model

The power generated by a wind turbine is directly related to the wind speed, as expressed by the following equation [37]:

$$P_w\left(\upsilon \right)=\left\{\begin{array}{l}0, \mathrm{for} \nu \langle {\nu }_{in} and \nu \rangle {\nu }_{out}\\ P_{wr}\left({\displaystyle \frac{\nu -{\nu }_{in}}{{\nu }_r-{\nu }_{in}}}\right) \mathrm{for} {\nu }_{in}\le \nu \le {\nu }_r \\ P_{wr} \mathrm{for} {\nu }_r\le \nu \le {\nu }_{out}\end{array}\right.$$

(38)

where ${\nu }_{in}$, ${\nu }_r$, and ${\nu }_{out}$ represent the cut-in, rated, and cut-out wind speeds of the turbine, respectively. $P_{wr}$ denotes the rated power output of the wind turbine.

3.2. Probability of Wind Power at Various Wind Speeds

According to these models, power will cease to be generated if the wind speed $\nu$ drops below threshold ${\nu }_{in}$ or exceeds threshold ${\nu }_{out}$. However, if the wind speed is within the range, ${\nu }_r\le \nu \le {\nu }_{out}$, the turbine will produce a power output of $P_{wr}$. The probabilities for these different wind speed zones can be calculated using the provided equations

$$f_w\left(P_w\right)\left\{P_w=0\right\}=1-\exp \left[-{\left({\displaystyle \frac{{\nu }_{in}}{\alpha }}\right)}^{\beta }\right]+\exp \left[-{\left({\displaystyle \frac{{\nu }_{out}}{\alpha }}\right)}^{\beta }\right]$$

(39)

$$f_w\left(p_w\right)\left\{p_w=p_{wr}\right\}=1-\exp \left[-{\left({\displaystyle \frac{{\nu }_r}{\alpha }}\right)}^{\beta }\right]+\exp \left[-{\left({\displaystyle \frac{{\nu }_{out}}{\alpha }}\right)}^{\beta }\right]$$

(40)

In contrast to the discrete zones, the power output of the wind turbine varies continuously in the bounds of a specific range ${\nu }_{in}\le \nu \le {\nu }_r$. Thus, the probability for this region can be modeled as follows [36]:

$$f_w\left(p_w\right)={\displaystyle \frac{\beta \left({\nu }_r-{\nu }_{in}\right)}{{\alpha }^{\beta }\ast P_{wr}}}{\left[{\nu }_{in}+{\displaystyle \frac{P_w}{P_{wr}}}\left({\nu }_r-{\nu }_{in}\right)\right]}^{\beta -1}\exp \left[-{\left({\displaystyle \frac{{\nu }_{in}+\frac{P_w}{P_{wr}}\left({\nu }_r-{\nu }_{in}\right)}{\alpha }}\right)}^{\beta }\right]$$

(41)

Solar radiation is represented by the lognormal distribution, taking into account parameters such as average radiation, standard deviation (σ), and mean (μ). Analyzing historical data assists in estimating these distribution parameters, which in turn aids in predicting energy production and evaluating the performance of PV units. Additionally, the conversion of solar irradiance to energy for photovoltaic plants can be described as follows:

$$P_s\left(G\right)=\left\{\begin{array}{l}{P}_{sr}\left({\displaystyle \frac{G^2}{G_{std}{R}_c}}\right) \mathrm{for} 0\le G\le R_c\\ P_{sr}\left({\displaystyle \frac{G^2}{G_{std}}}\right) \mathrm{for} G\ge R_c \end{array}\right.$$

(42)

In this context, the solar irradiance is referred to as $G_{std}$, under standard environmental conditions, and the rated power output of the PV power plant is referred to as $P_{sr}$. For a more detailed understanding of the uncertainty model of renewable energy sources (RESs), please consult the provided reference [33].

4. The Proposed Optimization Technique

4.1. Kepler Optimization Algorithm

The Kepler optimization algorithm, introduced by Mohamed Abdel-Basset et al. in March 2023 [38], takes inspiration from Johannes Kepler’s laws of planetary motion. This meta-heuristic approach is applied to solve optimization problems, particularly in evolutionary computation.

In this algorithm, the search space is depicted as the Sun and planets moving in elliptical orbits, mirroring Kepler’s first law. The eccentricity (e = 1), representing the degree of elongation of an ellipse, ranges from a line segment (e = 1) to a circle (e = 0). Different ellipse shapes are demonstrated in Figure 3 and Figure 4.

In this method, candidate solutions are represented by planets that systematically explore and exploit the search space. They occupy varied positions relative to the Sun, symbolizing the best solution, at different stages.

In our study, we replace “iteration” with “time” to align with principles from the solar system and cosmology. In the Kepler optimization algorithm (KOA), we follow the following principles:

-

The orbital period of each planet (candidate solution) is randomly determined from a normal distribution.

-

The eccentricity of the planet’s orbit is randomly chosen within the range of 0 to 1.

-

$r_3$, $r_4$, ${\overrightarrow{r}}_5$, and ${\overrightarrow{r}}_6$ are random numbers chosen within the range of 0 to 1.

-

Solution fitness is evaluated based on the objective function.

-

The Sun, the central star, symbolizes the best solution at each time.

The following section delineates the mathematical model and procedures employed in the (KOA).

• Phase 1: Initialization process

In this phase, a population of N planets is randomly distributed across d dimensions, representing the decision variables, using Equation (43).

$$X_i^j=X_{i,low}^j+rand \left[0, 1\right]\times \left(X_{i,up}^j-X_{low}^j\right),\left\{\begin{array}{c}i=1,2,\dots \dots ,N.\\ j=1,2,\dots \dots ,d.\end{array}\right.$$

(43)

$X_i$ denotes the i-^th planet within the search space, while ${X^j}_{i, up}$ and ${X^j}_{i, low}$ represent the lower and the upper bounds, respectively, for the i-^th decision variable. Additionally, the term rand (0, 1) signifies a randomly generated number ranging from 0 to 1.

The initial value of the orbital eccentricity (e) for each object is determined according to Equation (44):

$$e=rand [0, 1], i=1,\dots \dots ,N$$

(44)

($T_i$) is the orbital period for the ith object and can be determined by Equation (46):

$$T_i=\left|r\right|, i=1,\dots \dots ,N$$

(45)

• Phase 2: Determining the gravitational force (F)

The gravitational interaction between the Sun (denoted as Xs) and any given planet (denoted as Xi) is mathematically expressed by the universal law of gravitation, as follows:

$$F{g}_i(t)=e_i\times \mu (t)\times {\displaystyle \frac{\overline{M}s\times {\overline{m}}_i}{{{\overline{R}}^2}_i+\epsilon }}+r_1$$

(46)

In this equation, Ms and mi represent the normalized masses of Xs and Xi, respectively. The parameter ε denotes a small value, while $e_i$ signifies the eccentricity of the planet’s orbit, which varies between 0 and 1. μ stands for the universal gravitational constant. Ri denotes the normalized distance between Xs and Xi, and its expression is given by the following:

$$R_i(t)={‖Xs(t)-X_i(t)‖}_2=\sqrt{{\displaystyle {\sum }_{i=1}^N{(Xs(t)-X_i(t))}^2}}$$

(47)

The distance between Xs and Xi is denoted as ${‖X_s(t)-X_i(t)‖}_2$. The evaluation of the Sun’s mass and the object (i) at time (t) are calculated as follows:

$$M_s={\displaystyle \frac{fi{t}_s(t)-worst(t)}{{\displaystyle {\sum }_{K+1}^N(fi{t}_s(t)-worst(t))}}}$$

(48)

$$m_i=r_2{\displaystyle \frac{fi{t}_s(t)-worst(t)}{{\displaystyle {\sum }_{K+1}^N(fi{t}_s(t)-worst(t))}}}$$

(49)

$$fi{t}_s(t)=best(t)=k\in {\{}^{\min }, 1, 2, \dots ,N\}fi{t}_k(t)$$

(50)

$$Worst(t)=k\in {\{}^{\max }, 1, 2, \dots ,N\}fi{t}_k(t)$$

(51)

The function μ(t) regulates search precision by exponentially declining over time, and its expression is given by Equation (52):

$$\mu (t)={\mu }_0\times \exp \left(-\gamma {\displaystyle \frac{t}{T_{\max }}}\right)$$

(52)

In the context provided, γ stands for a constant, ${\mu }_0$ denotes an initial value, t represents the current iteration, and $T_{\max }$ symbolizes the maximum number of iterations.

• Phase 3: Calculating of planets velocity

Planet velocity adjusts based on its proximity to the Sun: it increases as it becomes closer in order to counteract the stronger gravitational forces and decreases as it becomes further away due to the weaker pull. This behavior, mathematically captured in Equations (53) and (54), enhances the search strategies in KOA.

$$\begin{gathered} v_i(t)=\mathcal{l}\times (2r_4\overrightarrow{X_i}-\overrightarrow{X_b})+\rho (\overrightarrow{X_a}-\overrightarrow{X_b})+(1-R_{i-norm}(t))\times F\times \overrightarrow{U\times }(\overrightarrow{X_{i,up}}-\overrightarrow{X_{i,low}}) \\ if R_{i-norm}(t)\le 0.5 \end{gathered}$$

(53)

Else

$$v_i(t)=r_4\times L\times (\overrightarrow{X_a}-\overrightarrow{X_i})+(1-R_{i-norm}(t))\times F\times \overrightarrow{U_1}\times (r_3\overrightarrow{X_{i,up}}-\overrightarrow{X_{i,low}})$$

(54)

$$\mathcal{l}=\overrightarrow{U}\times \mathrm{M}\times L,$$

(55)

$$L=[\mu (t)\times (M_S+m_i)(\frac{2}{R_i(t)+\epsilon }-\frac{1}{a_i(t)+\epsilon }){]}^{\frac{1}{2}}$$

(56)

$$\mathrm{M}=(r_3\times (1-r_4)+r_4)$$

(57)

$$\overrightarrow{U}=\left\{\begin{array}{cc}0& \overrightarrow{r_5}\le \overrightarrow{r_6}\\ 1& else\end{array}\right.$$

(58)

$$\mathcal{F}=\left\{\begin{array}{cc}1& if r_4\le 0.5\\ -1& else\end{array}\right.$$

(59)

$$\rho =(1-\overrightarrow{U})\times \overrightarrow{\mathrm{M}}\times L$$

(60)

$$\overrightarrow{\mathrm{M}}=(r_3\times (1-\overrightarrow{r_5})+{\overrightarrow{r}}_5)$$

(61)

$$\overrightarrow{U_1}=\left\{\begin{array}{cc}0& r_5\le r_4\\ 1& else\end{array}\right.$$

(62)

$$\overrightarrow{U_2}=\left\{\begin{array}{cc}0& r_3\le r_4\\ 1& else\end{array}\right.$$

(63)

In these equations, $\overrightarrow{v_i}(t)$ denotes the velocity of object (i) at time (t), where ${\overrightarrow{X}}_i$ represents the object. Additionally, ${\overrightarrow{X}}_a$ and ${\overrightarrow{X}}_b$ represent solutions chosen randomly from the population, while $M_s$ and $m_i$ represent their respective masses. μ(t) represents the gravitational constant.

$a_i$ represents the semi-major axis of the ellipse defined by Kepler’s third law, as mentioned in Equation (64).

$$a_i(t)=r_3\times {\left[T_i^2\times \frac{\mu (t)\times (M_S+m_i)}{4{\pi }^2}\right]}^{\frac{1}{3}}$$

(64)

where $T_i$ is the orbital period of object i and is calculated using Equation (45). $R_{i-nom}(t)$ defines the normalizing Euclidian distance between $X_s$ and $X_i$, and it is remodeled as follows:

$$R_{i-norm}(t)={\displaystyle \frac{R_{i-norm}\left(t\right)-R_i}{\mathit{min}(R\left(t\right))-\mathit{max}(R(t\left)\right)}}$$

(65)

Equation (65) computes the step percentage for adjusting each object’s velocity. When $R_{i-norm}(t)$ ≤ 0.5, it signifies that the object is in close proximity to the Sun and requires an increase in velocity to counteract gravitational drift toward the Sun, driven by its strong gravitational force. Conversely, if the value exceeds 0.5, the object will decelerate.

• Phase 4: Preventing the local optimum

The celestial bodies in the solar system typically orbit the Sun in an anticlockwise direction while rotating on their axes. However, some bodies do orbit clockwise. The proposed algorithm utilizes this phenomenon to avoid local optima. The KOA incorporates a flag F to adjust the search direction, thereby dynamically improving the exploration of the search space.

• Phase 5: Updating objects’ locations

As objects orbit the Sun, they exhibit a cyclic pattern of approaching and retreating. The KOA algorithm imitates this behavior through an exploration–exploitation strategy. In the exploration phase, distant objects are extensively surveyed to uncover new solutions, while those closer to the Sun assist in refining the search for optimal solutions.

$${\overrightarrow{X}}_i(t+1)={\overrightarrow{X}}_i(t)+\mathcal{F}\times {\overrightarrow{v}}_i\times (F{g}_i(t)+\left|r\right|)\times \overrightarrow{U}\times ({\overrightarrow{X}}_s(t)-{\overrightarrow{X}}_i(t))$$

(66)

In this scenario, $X_i(t+1)$ represents the updated position of object i at time t + 1, while ${\nu }_i(t)$ indicates the velocity required for object i to reach its new position. $X_s(t)$ denotes the best-known position of the Sun up to the current time.

• Phase 6: Updating the distance of objects from the Sun

In KOA, optimization focuses on exploiting planets near the Sun and exploring those farther away. Decision-making depends on a gradually evolving regulatory parameter. A high value of (h) prioritizes exploration, promoting orbital distancing from the Sun. Conversely, a small (h) favors exploitation, focusing on areas near the best-known solution.

$${\overrightarrow{X}}_i(t+1)={\overrightarrow{X}}_i(t)\times {\overrightarrow{U}}_1+(1-{\overrightarrow{U}}_1)\times \left(\frac{{\overrightarrow{X}}_i(t)+X_s+{\overrightarrow{X}}_a}{3.0}+h\times \left(\frac{{\overrightarrow{X}}_i(t)+X_s+{\overrightarrow{X}}_a}{3.0}-{\overrightarrow{X}}_b(t)\right)\right)$$

(67)

where (h) serves as an adaptive factor, controlling the distance between the current planet and the Sun at time (t), as elaborated below:

$$h={\displaystyle \frac{1}{e^{\eta r}}}$$

(68)

η is a factor that decreases linearly from 1 to −2, as shown below:

$$\eta =(a_2-1)\times r_4+1$$

(69)

$a_2$ represents a cyclic control parameter that gradually decreases from −1 to −2 over ($\overline{T}$) cycles throughout the entire optimization process, as depicted in Equation (70).

$$a_2=-1-1\times \left({\displaystyle \frac{t\times \frac{T_{\max }}{T}}{\frac{T_{\max }}{T}}}\right)$$

(70)

• Phase 7: Elitism

This step aims to employ an elitist approach by preserving the optimal positions of both the Sun and the planets. This procedure is succinctly represented by Equation (71).

$${\overrightarrow{X}}_{i,new}(t+1)=\left\{\begin{array}{ll}{\overrightarrow{X}}_i(t+1),& if f({\overrightarrow{X}}_i(t+1)\le {\overrightarrow{X}}_i(t))\\ {\overrightarrow{X}}_i(t),& else\end{array}\right.$$

(71)

4.2. Non-Dominated Sorting Kepler Optimization Algorithm (NSKOA)

In the original version of the Kepler optimization algorithm (KOA), the algorithm was designed to solve mono-objective optimization problems, where a single global best solution was determined as the “Sun” (XS) based on fitness values. The algorithm updated the positions and velocities of the objects (planets) accordingly, aiming to improve their fitness relative to this singular objective.

To extend KOA into a multi-objective framework, the key modification was the integration of the non-dominated sorting procedure, allowing the algorithm to handle multiple conflicting objectives simultaneously. The key changes are as follows:

• Non-Dominated Sorting: This mechanism is essential for navigating trade-offs in multi-objective problems by ranking solutions based on Pareto dominance. In the modified version, after initializing the population and calculating the objective functions, non-dominated sorting is performed to assign ranks to solutions. Non-dominated solutions, which are not outperformed in all objectives, receive a rank of 1, while subsequent ranks are assigned iteratively to solutions that are dominated by others. This ranking system helps preserve a diverse set of optimal solutions, ensuring that different trade-offs between objectives are explored.
  • Crowding distance (CD): An important aspect of this approach is the calculation of crowding distance, which measures the proximity of solutions to their neighbors in the objective space. A higher crowding distance indicates a less crowded region, helping to maintain diversity among solutions and favoring those that are well-distributed along the Pareto front.

By evaluating solutions based on both rank and crowding distance, the algorithm ensures a balance between convergence toward the Pareto front and the preservation of solution diversity across the front. This process ensures that, instead of converging to a single global best solution, a set of non-dominated solutions is maintained, offering multiple viable options that address different trade-offs between the conflicting objectives.

2.

Best Compromise Solution (BCS): The best compromise solution (BCS) is introduced in the multi-objective version, calculated using Equation (72). This BCS serves as a reference point to guide the search toward an optimal trade-off among objectives, considering both ranks and crowding distances.

3.

Elitism and Population Combination: After each iteration, the algorithm combines the newly generated population (Npi) with the previous population (Pi) to form an augmented population (Upi). Non-dominated sorting is again applied to this combined population. From this merged pool, the best N elitist objects are selected to form the next generation. This ensures that high-quality solutions from both the new and previous generations are preserved, promoting convergence to the Pareto front. Algorithm 1 below represents the pseudocode of NSKOA. More details on the process how to transform KOA of NSKOA are given at Appendix A.

Algorithm 1: Pseudocode of NSKOA

Step 1: Define input power system data (line data þ bus data) and identify the control variable limits and number of variables. Set KOA parameters N, Tmax, and ${\mu }_0$.

Step 2: Initialize objects population with random position, orbital eccentricities, and orbital periods using Equations (43), (44) and (45), respectively.

Step 3: Run a power flow algorithm based on the Newton Raphson method to calculate the value of the objective functions for the initial population.

Step 4: Perform non-dominated sorting: Calculate ranks (RK) and crowding distance (CD) using the eq() of population using the proposed PFA with the sorting and crowding distance calculation procedure. Calculate the best compromise solution (BCS) using Equation (73).

Step 5: While (t < Tmax): Update ei … i = 1,2, …, N, best(t), worst(t), and μ(t), using Equations (50), (51) and (52), respectively.

Step 6: For i = 1: N Pi = population Calculate the gravitational force between the Sun and the object i using Equation (46). Calculate the Euclidian distance between the Sun and the object i using Equation (48). Calculate the velocity of the object Xi using Equation (53). Generate two random numbers r and r1 between 0 and 1. If r > r1 Update the position of the planet. Update the object position using Equation (66). Else Update the distance between the planet and the Sun. Update the object position using Equation (67). End if

Step 7: Run power flow algorithm based on the Newton Raphson algorithm to calculate the objective functions values for the new population (Npi). Combine new population (Npi) with previous population (Pi) to form Upi   Upi = Npi U Pi.

Step 8: Perform non dominated sorting: Calculate ranks (RK) and crowding distance (CD) of population using the proposed PFA with the sorting and crowding distance calculation procedure.

Step 9: Extract N elitist objects from Upi.

Step 10: Generate the Pareto optimal front and extract the best compromise solution.

End for End while

4.3. The Best Compromise Solution (BCS)

Fuzzy set theory is frequently employed to effectively select a candidate Pareto-optimal solution from numerous options along the Pareto front. Given the inherent irrationality of decision-makers, the i-^th objective function of a solution within the Pareto-optimal set, denoted as $f_i$ [38], is expressed through a membership function ${\mu }_i$, defined as follows:

$${\mu }_i=\left\{\begin{array}{l}1 f_i\le f_i^{\min }\\ {\displaystyle \frac{f_i^{\max }-f_i}{f_i^{\max }-f_i^{\min }}} f_i^{\min }\le f_i\le f_i^{\max }\\ 0 f_i\ge f_i^{\max }\end{array}\right.$$

(72)

where $f_i^{\max }$ and $f_i^{\min }$ represent the maximum and minimum values of the $i-th$ objective function, respectively.

The normalized membership function ${\mu }^k$ is calculated for each non-dominated solution k as follows:

$${\mu }^k={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^{Nof}{{\mu }_i}^k}}{{\displaystyle {\sum }_{j=1}^M{\displaystyle {\sum }_{i=1}^{Nof}{{\mu }_i}^j}}}}$$

(73)

where the number of non-dominated solutions is denoted as M. The best compromise solution is determined as the one with the highest value of ${\mu }^k$. By organizing all solutions in descending order based on their membership function, a priority list of non-dominated solutions is generated. This prioritized list serves as guidance for the decision-maker, aiding in navigating through the current operational circumstances.

5. Results and Discussion

To address stochastic OPF problems, an analysis was conducted on both the conventional (base case) and the modified IEEE 57-bus network. The base case was simulated using the mono-objective KOA to demonstrate the effects of RES and FACTs on the system. To determine the optimal control variables, sizing, and location of SVC and TCSC, 17 tap changers, as well as the power and voltage of generators, were used as control variables. SVC placement was considered among 50 load buses, and TCSC could be situated across 81 branches, totaling 161 control variables as represented in

Table 1. Control variables of the test system.

Elements	IEEE 57-Bus Test System
No of buses	57
No of branches	81
No of generators	07
No of thermal generators	04
No of RES generators	03
No of load buses	50
No of control variables	161
Initial active and reactive load demand	1250.80 MW; 336.40 Mvar

.

The proposed algorithm was implemented using MATLAB software and simulations were conducted on a personal computer with an Intel Core™ i7-8300H 2.22 GHz processor. To determine the optimal population size for the NSKOA algorithm, empirical tests were conducted with different population sizes, taking into account the search space complexity and the number of control variables. Population sizes of 100, 200, and 300 were tested. Although specific test results are not provided here, the best outcomes were achieved with a population size of 200 individuals, which was then used for all simulation cases. For equitable comparison, the control variables of the test system were treated as continuous variables.

5.1. Test-System: Conventional and Modified IEEE 57 BUS Network

The IEEE 57-bus test system consists of seven power plants installed at buses 1, 2, 3, 6, 8, 9, and 12; eighty transmission lines of which 17 are equipped with tap changer transformers, and three parallel compensators are installed at buses 18, 25, and 53, respectively. The complete data are available in [23]. The modified IEEE-57 bus system consists of a combined production of a solar–hydro power generator replacing a thermal power plant at bus 6 and a solar PV generator at bus 9, as well as wind generators at bus 12. The parameters of the probability density function (PDF) and the cost coefficients of the renewable energy sources are detailed in

Table 2. PDF parameters of renewable generators.

Wind-power unit
No. of turbines	Rated power, $P_{wr}$ (MW)		Weibull PDF parameters
25	75		l = 9; p = 2
Photovoltaic plant
Rated power, ${\mathit{P}}_{\mathit{s}\mathit{r}}$ (MW)		Lognormal PDF parameters
50		μ = 5.2 σ = 0.6
Combined solar and small hydro power
Photovoltaic rated power ${\mathit{P}}_{\mathit{s}\mathit{r}}$ (MW)		Lognormal PDF parameters
45		μ = 5.0 σ = 0.6
Small hydro rated power ${\mathit{P}}_{\mathit{h}\mathit{r}}$ (MW)		Gumbel PDF parameters
5		λ = 15 γ = 1.2

[20].

The study was conducted using three scenarios:

• Base Case Scenario: In this case, the conventional IEEE 57 bus was simulated in order to show the impact of RESs and FACTs devices on the four optimization cases (cost, power losses, voltage deviation, and voltage stability index) in the next two scenarios.
• Scenario number 1: This study was conducted on the modified IEEE 57 system after the integration of RES sources.
• Scenario number 2: This study was conducted on the modified IEEE 57 system after the integration of RES sources and FACTs devices.

The findings of the case studies are presented in a tabulated format. Figure 5, Figure 6, Figure 7 and Figure 8 exhibit the characteristics of the probabilities of RESs, while Figure 9, Figure 10, Figure 11 and Figure 12 represent the available real power at the RES units.

5.2. Case 1: Total Generation Cost and the Investment Cost of FACTS Optimization

A comprehensive investigation into the optimal power flow for the IEEE 57-bus system was undertaken, examining three distinct scenarios. Initially, the conventional configuration of the system was analyzed, followed by two additional cases: one integrating renewable energy sources (RESs) and the other involving RES integration alongside Flexible Alternating Current Transmission System (FACTS) devices. The primary objective was to minimize the cost of power production in the conventional IEEE 57 system designed as the base case scenario, and then with the integration of RES designated as scenario 1, while simultaneously considering the investment cost associated with FACTS devices, defined as scenario 2.

In the base case scenario, the production cost was recorded at 5570.956 (USD/h), with an emission rate of 234.75 tons/h. Upon the integration of RESs into the system, the production cost decreased to 5217.635 (USD/h) saving 353.32 (USD/h), approximately 6.34% of the power production cost, alongside a reduction in emissions by approximately 22.76%, to 181.294 tons per hour (scenario 1). The subsequent deployment of FACTS devices further reduced the production cost to 5208.97 (USD/h), saving approximately 8.8 (USD/h) representing a marginal saving of approximately 0.17% compared to the RES-integrated scenario. However, this enhancement incurred an additional cost of 288.973 (USD/h) for FACTS deployment, yielding a total cost of 5497.94 (USD/h), but still saving 73 (USD/h) from the total cost compared to the base case.

Notably, the best compromise solution (BCS) yielded a production cost of 5211.722 (USD/h), and a FACTS cost of 82.5115 (USD/h), with a total cost of 5294.231 (USD/h) reflecting savings of approximately 4.98% compared to the base case.

Figure 13 shows the Pareto front of case 1, where NSKOA provides a well-distributed front, with the BCS solution positioned almost in the center, while the optimal result and control variables of case 1 are represented in

Table 3. Optimal result and control variables of case 1.

Variables	Base Case	Scenario 1	Scenario 2		BCS
Pg1	331.516	552.331	550.564		551.9723
Pg2	99.986	100.000	100.000		100.000
Pg3	76.635	76.6229	76.6240		76.6183
Pg6	99.997	100.000	100.000		99.9995
Pg8	53.916	50.5657	50.8772		50.0453
Pg9	160.173	199.9998	200.000		199.999
Pg12	209.847	210.000	209.999		209.998
Vg1	1.100000	1.09983	1.099983		1.099990
Vg2	1.095161	1.092973	1.093325		1.093092
Vg3	1.079965	1.070328	1.071347		1.070841
Vg6	1.065115	1.043710	1.040434		1.040998
Vg8	1.064082	1.034728	1.029753		1.031650
Vg9	1.050054	1.027750	1.032134		1.028320
Vg12	1.06068	1.03707	1.03637		1.03660
T(4,18)	1.0539	1.0084	1.0507		1.0724
T(4,18)	0.9958	0.9868	0.9727		0.9664
T(21,20)	1.0538	1.0122	1.0010		1.0133
T(24,25)	0.9601	0.939	0.9495		0.9457
T(24,25)	1.0158	0.9375	0.9622		0.9585
T(24,26)	1.0020	0.9820	0.9787		0.983
T(7,29)	1.0006	0.9599	0.9659		0.9634
T(34,32)	0.9717	0.9184	0.9135		0.9107
T(11,41)	0.9220	0.9000	0.9158		0.9123
T(15,45)	0.9835	0.9887	1.0014		0.9959
T(14,46)	0.9744	0.9648	0.9800		0.9744
T(10,51)	0.9978	0.9723	0.9782		0.9775
T(13,49)	0.9503	0.9331	0.9784		0.9771
T(11,43)	0.9997	0.9591	0.9872		0.9776
T(40,56)	1.0168	0.9861	0.9509		0.9432
T(39,57)	0.9744	0.9549	0.9577		0.9617
T(9,55)	0.9922	0.9643	0.97113		0.9679
Cost (USD/h)	5570.956	5217.635	5208.97		5211.722
TFcost (USD/h)	—	—	288.973		82.5115
Emission	234.75	181.294	180.125		181.06
RPL (MW)	43.958	38.7203	37.263		37.833
TVD (pu)	1.4994	1.5821	1.6739		1.6739
VSI (pu)	0.2899	0.2757	0.2533		0.2533
Optimal size and location of SVC–TCSC
			svc(13)	50.000	svc(42)	0.7796
			svc(16)	24.050	svc(44)	6.5088
			svc(21)	7.135	svc(21)	3.5401
			svc(35)	15.038	svc(35)	12.438
			svc(52)	5.825	svc(52)	0.6752
			svc(54)	2.567	svc(54)	1.5976
			Tcsc(21,22)	0.0377	Tcsc(30,31)	0.3588
			Tcsc(30,31)	0.3976	Tcsc(13,49)	0.1184
			Tcsc(13,49)	0.1260	Tcsc(29,52)	0.1430
			Tcsc(29,52)	0.1393	Tcsc(56,42)	0.2421
			Tcsc(56,42)	0.1939	Tcsc(38,49)	0.0138
			Tcsc(38,49)	0.0332

P_Gi (MW), V_gi (p.u.), T_(I,j) (p.u.),Q svc_(i) (MVAr), X Tcsc_(I,j) (p.u.).

. The results show that the placement of SVC devices at buses 13, 16, 21, 35, 52, and 54 played a significant role in managing reactive power and stabilizing voltage profiles. The addition of SVCs helped improve voltage levels, particularly in buses with high reactive power demand, ensuring more efficient power flow across the network. Similarly, the installation of TCSC devices on transmission lines (21–22), (30–31), (13–49), (29–52), (56–42), and (38–49) optimized power transfer capability by reducing line reactance and enhancing reactive power compensation, leading to more stable voltage regulation and an overall reduction in power losses. The combination of SVCs and TCSCs helped balance reactive power more effectively, improving voltage profiles and reducing strain on power generation units by minimizing unnecessary reactive power generation, thereby contributing to cost optimization. Figure 14 and Figure 15 illustrate the voltage profile in the load buses and the generated reactive power, respectively, in the same case. All values are within their limits, indicating that the constraints are completely satisfied. Importantly, the combination of RES integration and FACTS deployment not only reduced power production costs and emissions but also enhanced the voltage stability index, significantly mitigating power losses and voltage deviations. This synergy underscores the potential for designing and operating future power systems that ensure a reliable and sustainable electricity supply.

Table 4. Statistical comparison of case 1 and case 2 for the first scenario.

Case 1			Case 2
Algorithms		Results ($/h)	Algorithms		Results (MW)
BMO	[20]	5300.457	BMO	[20]	20.7850
MFO	[20]	5316.140	MFO	[20]	21.3031
PSO	[20]	5417.538	PSO	[20]	21.3621
GTO	[39]	5260.009	GTO	[39]	19.7703
AEO	[39]	5260.249	AEO	[39]	19.7633
INFO	[39]	5259.204	INFO	[39]	19.7040
NSKOA		5217.635	NSKOA		16.8360

Bold indicates the best solutions found so far.

represents a statistical comparison of the proposed method with six other methods used in case 1 and case 2 for the first scenario, for instance, the barnacles mating optimizer (BMO), moth-–flame optimization algorithm (MFO), particle swarm optimization (PSO), efficient optimization algorithm based on weighted mean of vectors optimization (INFO), and artificial ecosystem-based optimization (AEO).

5.3. Case 2: Real Power Loses and the Investment Cost of FACTS Optimization

In a parallel investigation, the optimization of real power losses (RPLs) within the same system configurations was pursued. In the base case scenario, RPL stood at 18.022 MW, which decreased to 17.729 MW following the integration of RESs, achieving a reduction of approximately 1.62% (Objective 1). Subsequent intervention with FACTS devices further diminished RPL to 16.398 MW, demonstrating an additional reduction of approximately 9.05% compared to the base case and 7.45% compared to the RES-integrated scenario. However, this enhancement incurred an additional cost of 227.86 (USD/h) due to FACTS deployment (scenario 2). Conversely, the base case scenario led to an RPL of 16.836 MW, with a FACTS cost of 68.945 (USD/h).

Figure 16 shows the Pareto front of case 2, where NSKOA provides a well-distributed front, with the BCS solution positioned almost in the center, as shown in

Table 5. Optimal result and control variables for case 2.

Control Variables	Base Case	Scenario 1	Scenario 2		BCS
Pg1	198.3536	301.5189	297.4671		299.6088
Pg2	19.6385	8.424944	7.759523		8.143402
Pg3	136.582	140.000	139.9811		139.962
Pg6	94.5306	99.9969	99.9978		99.9930
Pg8	319.912	308.590	311.9925		309.931
Pg9	199.070	199.999	199.999		200.000
Pg12	209.994	209.998	210.000		209.997
Vg1	1.074194	1.075784	1.073172		1.075249
Vg2	1.067047	1.068157	1.066069		1.067677
Vg3	1.061879	1.064005	1.063692		1.063580
Vg6	1.063364	1.059364	1.059340		1.059306
Vg8	1.075176	1.061550	1.064236		1.064473
Vg9	1.051757	1.043234	1.051782		1.047636
Vg12	1.040481	1.037484	1.037586		1.038832
T(4,18)	1.042820	0.964453	0.955661		0.950783
T(4,18)	0.969822	1.026266	1.035669		1.040006
T(21,20)	1.00078	1.019229	1.006507		1.014671
T(24,25)	0.946581	0.937832	0.960908		0.954601
T(24,25)	0.962140	0.971551	0.975584		0.972877
T(24,26)	1.025195	1.013087	1.009955		1.016856
T(7,29)	0.991647	0.983174	0.990604		0.986726
T(34,32)	0.943693	0.927002	0.929053		0.932515
T(11,41)	0.922844	0.900000	0.930783		0.919565
T(15,45)	0.982334	0.983402	0.993976		0.984972
T(14,46)	0.963354	0.963000	0.984149		0.975964
T(10,51)	0.975757	0.969271	0.986212		0.977325
T(13,49)	0.937087	0.936787	0.991151		0.980100
T(11,43)	0.975505	0.977018	0.988254		0.982823
T(40,56)	1.003564	1.006911	0.958281		0.960962
T(39,57)	0.962797	0.968382	0.957930		0.969004
T(9,55)	0.988405	0.986404	0.991583		0.984642
Optimal size and location of SVC–TCSC
			svc(13)	17.3838	svc(35)	8.28284
			svc(14)	11.2275	svc(38)	6.5553
			svc(35)	12.2902	svc(50)	2.08968
			svc(38)	16.354	svc(53)	1.7785
			svc(50)	11.123	svc(54)	0.6111
			svc(53)	6.240
			Tcsc(9,12)	0.0573	Tcsc(9,12)	0.0269
			Tcsc(1,16)	0.0308	Tcsc(1,16)	0.0042
			Tcsc(30,31)	0.3034	Tcsc(47,48)	0.0039
			Tcsc(47,48)	0.0175	Tcsc(13,49)	0.1356
			Tcsc(13,49)	0.1465	Tcsc(56,42)	0.0228
			Tcsc(38,48)	0.0371	Tcsc(38,48)	0.0347
Cost (USD/h)	10,979.12	10,154.82	10,269.41		10,198.18
RPL (MW)	18.7160	17.7299	16.3981		16.8366
TVD (pu)	1.5782	1.5298	1.9708		1.6799
VSI (pu)	0.2766	0.2759	0.2468		0.2593

P_Gi (MW), V_gi (p.u.), T_(I,j) (p.u.), Q svc_(i) (MVAr), X Tcsc_(I,j) (p.u.).

. The strategic placement of SVCs and TCSCs significantly influenced power loss reduction by optimizing reactive power flow. Notably, buses 50 and 53 were equipped with SVCs, highlighting their critical role in reactive power support and reducing power losses. Similarly, the results show that in the IEEE 57 BUS network, the branches (9–12) and (13–49) are the regions that need TCSC device installation to effectively control line impedance and improve power flow efficiency, which reduces power losses in these regions by reducing line reactance. The presence of TCSCs on these repeated branches indicates their importance in minimizing congestion and reducing transmission losses, making them key points for reactive power compensation and power loss mitigation. These findings underscore the efficacy of FACTS deployment in significantly mitigating real power losses within the system, even with an added expense, presenting a nuanced trade-off between loss reduction and investment expenses. Figure 17 and Figure 18 illustrate the voltage profile in the load buses and the generated reactive power, respectively, in the same case. All values are within their limits, indicating that the constraints are completely satisfied.

5.4. Case 3: Total Voltage Deviation and the Investment Cost of FACTS Optimization

In a comprehensive analysis encompassing the optimization of total voltage deviation (TVD) within the same system configurations, notable improvements were observed. In the base case scenario, TVD registered at 0.7029 per unit (p.u.), which decreased to 0.6842 p.u. following the integration of RESs, marking a reduction of approximately 2.66% for scenario 1 (Objective 1). Subsequent intervention with FACTS devices yielded a significant improvement, reducing TVD to 0.2138 p.u., reflecting a substantial enhancement of approximately 69.5% compared to the conventional configuration of the IEEE 57 bus and 68.79% compared to the RES-integrated scenario. However, this advancement incurred an additional cost of 237.38 (USD/h) due to FACTS deployment. Conversely, the base case scenario led to TVD of 0.3435 p.u., with a FACTS cost of 73 (USD/h).

Figure 19 shows the Pareto front for case 3, where NSKOA provides a well-distributed front, with the best compromise solution (BCS) positioned near the center. Figure 20 and Figure 21 present the voltage profile at the load buses and the generated reactive power, respectively, for the same case. All values remain within permissible limits, confirming that the constraints are fully satisfied. The numerical results presented in

Table 6. Optimal result and control variables for case 3.

Control Variables	Base Case	Scenario 1	Scenario 2		BCS
Pg1	198.3536	545.1541	425.9918		519.1516
Pg2	19.6385	42.6645	39.6999		57.4019
Pg3	136.582	110.9823	118.1849		96.2114
Pg6	94.5306	14.0295	58.4176		0.0000
Pg8	319.912	239.4711	224.5698		262.5789
Pg9	199.070	171.4433	199.999		182.9161
Pg12	209.994	164.5253	209.999		168.6121
Vg1	1.074194	1.04032	1.02062		1.02712
Vg2	1.067047	1.02816	1.01073		1.02244
Vg3	1.061879	1.02422	1.00583		1.02102
Vg6	1.063364	1.00302	1.00025		1.00153
Vg8	1.075176	1.02873	1.02203		1.03107
Vg9	1.051757	1.01095	1.01148		1.01119
Vg12	1.040481	1.01479	1.00556		1.00799
T(4,18)	1.042820	0.97587	0.955661		0.97037
T(4,18)	0.969822	1.04991	1.05236		1.04938
T(21,20)	1.00078	0.96580	0.96297		0.96335
T(24,25)	0.946581	0.96748	0.96883		0.97148
T(24,25)	0.962140	0.96013	0.98919		0.98409
T(24,26)	1.025195	1.03451	1.02960		1.03335
T(7,29)	0.991647	0.95720	0.96383		0.95958
T(34,32)	0.943693	0.92014	0.95538		0.95523
T(11,41)	0.922844	0.900000	0.92422		0.90000
T(15,45)	0.982334	0.93812	1.02301		0.98128
T(14,46)	0.963354	0.98081	1.00294		0.99451
T(10,51)	0.975757	0.99571	0.99637		0.99666
T(13,49)	0.937087	0.90023	0.90764		0.90000
T(11,43)	0.975505	0.97609	1.00385		1.00178
T(40,56)	1.003564	1.02057	0.92593		0.94376
T(39,57)	0.962797	0.90000	0.94573		0.95405
T(9,55)	0.988405	0.98505	1.01831		0.98132
Optimal size and location of SVC–TCSC
			svc(14)	21.585
			svc(21)	12.785
			svc(35)	30.907	svc(35)	28.511
			svc(44)	7.757
			svc(53)	5.156
			svc(54)	9.580
			Tcsc(19,20)	0.3155	Tcsc(19,20)	0.3472
			Tcsc(21,22)	0.0936	Tcsc(21,22)	0.0657
			Tcsc(37,39)	0.0303	Tcsc(37,39)	0.0216
			Tcsc(36,40)	0.0373	Tcsc(36,40)	0.0199
			Tcsc(56,42)	0.1961
			Tcsc(38,49)	0.1416
Cost (USD/h)	10,979.12	8102.0217	7546.581		8671.4377
TFcost($/h)	—	—	237.380		73.0604
Emission	17.7160	37.4701	26.0637		36.0721
RPL (MW)	1.5782	0.6842	0.2138		0.3436
TVD (pu)	0.2766	0.29440	0.2945		0.29601

P_Gi (MW), V_gi (p.u.), T_(I,j) (p.u.),Q svc_(i) (MVAr), X Tcsc_(I,j) (p.u.).

illustrate that the deployment of SVCs and TCSCs effectively reduced voltage deviation, enhancing voltage stability across the network. SVCs were installed at buses 14, 21, 35, 44, 53, and 54, providing crucial reactive power support to improve voltage profiles at these locations. Notably, in the IEEE 57-bus network, bus 35 required a significant SVC installation, which played a key role in reducing voltage deviation by supplying or absorbing reactive power as needed. This capability helps maintain voltage levels within desired limits, especially during fluctuations in load.

The TCSCs installed on branches (19–20), (21–22), (37–39), and (13–49) facilitated better control of line reactance and optimized the flow of reactive power, which mitigated voltage drops along transmission lines, contributing to improved voltage levels and further aiding voltage regulation. The consistent presence of SVC at bus 35 and TCSC on branch (13–49) illustrates their pivotal role in mitigating voltage deviations and ensuring stable reactive power management throughout the IEEE 57-bus network. These findings highlight the efficacy of FACTS deployment in reducing total voltage deviation within the system, underscoring the trade-off between improved voltage stability and FACTS deployment cost.

5.5. Case 4: Voltage Stability Index and the Investment Cost of FACTS Optimization

In this phase, the focus shifted towards optimizing and enhancing the voltage stability index within the examined system configurations. The base case exhibited a voltage stability index of 0.2757 p.u., which notably improved to 0.2018 p.u. following the integration of both RESs and FACTS technologies, representing a significant improvement of approximately 26.80%. However, this improvement came at an additional cost of USD 62.04 per hour due to FACTS deployment. Interestingly, in the BCS, the voltage stability index of 0.2074 p.u. was achieved with a significantly lower FACTS cost of USD 10.38 per hour, showcasing comparable performance with substantial savings of approximately USD 52 per hour, equivalent to approximately 83.78% in investment expenses of FACTS.

Figure 22 shows the Pareto front for case 4, where NSKOA provides a well-distributed front, with the best compromise solution (BCS) positioned near the center.

Table 7. Optimal result and control variables for case 4.

Control Variables	Base Case	Scenario 1	Scenario 2		BCS
Pg1	320.9496	552.3335	402.9396		516.5286
Pg2	44.2970	100.000	32.0975		93.31830
Pg3	70.9802	76.6229	135.201		92.04518
Pg6	49.1657	99.9997	0.45037		0.000000
Pg8	271.1674	50.5656	295.800		196.9922
Pg9	197.0577	199.998	200.000		180.3202
Pg12	152.410	209.999	208.955		210.000
Vg1	1.076947	1.099834	1.037665		1.038911
Vg2	1.064682	1.092973	1.021399		1.029121
Vg3	1.040932	1.070328	1.021756		1.021645
Vg6	1.035024	1.043710	1.008789		1.006930
Vg8	1.050021	1.034728	1.031409		1.029410
Vg9	1.030837	1.027750	1.014471		1.011513
Vg12	1.023650	1.037090	1.018090		1.014264
T(4,18)	0.9726	1.008410	1.022652		1.017152
T(21,20)	1.07613	0.986867	0.901523		0.922051
T(24,25)	0.98528	1.012262	1.069454		1.017742
T(24,25)	0.96546	0.939723	0.998648		1.007304
T(24,26)	1.04711	0.937514	0.992255		0.997993
T(7,29)	0.959181	0.982061	1.054628		1.084371
T(34,32)	0.90985	0.959973	0.939237		0.931607
T(11,41)	0.95951	0.918400	0.900827		0.900019
T(15,45)	0.97504	0.900000	0.900589		0.900000
T(14,46)	0.95079	0.988708	0.955885		0.954432
T(10,51)	0.95684	0.964576	0.935021		0.934472
T(13,49)	0.90762	0.970265	0.941255		0.939803
T(11,43)	0.95171	0.934064	0.900053		0.900661
T(40,56)	0.98339	0.95915	1.053052		0.975194
T(39,57)	1.09573	0.986130	1.083148		1.058040
T(9,55)	1.032468	0.954959	1.012770		1.037772
		0.964379	0.989115		0.988271
Optimal size and location of SVC–TCSC
			svc(16)	0.75793
			svc(21)	0.66050	svc(32)	0.39169
			svc(28)	0.01766	svc(28)	0.00171
			svc(54)	0.0710	svc(54)	0.00202
			Tcsc(18,19)	0.5480	Tcsc(18,19)	0.4040
			Tcsc(24,25)	0.9840	Tcsc(24,25)	0.9840
			Tcsc(25,30)	0.1616	Tcsc(25,30)	0.1616
			Tcsc(30,31)	0.3952	Tcsc(30,31)	0.3964
			Tcsc(37,38)	0.0807	Tcsc(37,38)	0.0807
			Tcsc(11,41)	0.5992	Tcsc(11,41)	0.1752
Cost (USD/h)	9193.49	5217.74	9748.93		6993.24
TFcost(USD/h)	—	—	62.045		10.380
Emission	29.635	38.719	24.645		38.4040
RPL (MW)	1.4720	1.5869	1.4312		1.3522
TVD (pu)	0.2757	0.2757	0.2018		0.2074

shows that the integration of SVCs and TCSCs significantly influenced the voltage stability index. Notably, the SVCs at buses 16, 21, 28, and 54 provided vital reactive power compensation, especially at bus 28, crucial for keeping the voltage profile at the desired levels. This ability reduces the risk of voltage collapse and improves the voltage stability index as shown in Figure 23 and Figure 24, which present the voltage profile at the load buses and the generated reactive power where all values remain within permissible limits.

The results indicate that branches (18–19) and (24–25) required larger TCSCs due to high reactive power demand, which directly impacted voltage stability, by enhancing power transfer capability, helping to stabilize voltage levels, and preventing conditions that could lead to instability. The placement of TCSCs on branches (30–31) and (37–38) further emphasizes their crucial role in optimizing voltage stability across the IEEE 57 network. The branches (18–19) and (24–25) demonstrated the need for larger TCSCs as these lines were under significant reactive power demand, directly impacting voltage stability. TCSCs improve the power transfer capability of lines, which can help stabilize voltage levels and prevent conditions leading to voltage instability. The TCSCs on branches (30–31) and (37–38) underscore their critical role in optimizing voltage stability across the IEEE 57 network. These findings underscore the effectiveness of FACTS deployment in enhancing voltage stability while highlighting the importance of cost considerations in optimizing system performance.

6. Conclusions

This paper presents a comprehensive study on the optimization of power flow through an exploration of a multi-objective optimization strategy that can simultaneously address a range of critical objectives using the non-dominated sorting Kepler optimization algorithm (NSKOA), focusing on the integration of renewable energy sources (RESs) and FACTS devices into the electrical network. The proposed NSKOA has demonstrated its effectiveness in achieving significant improvements in several key performance indicators. Notably, the integration of RESs and FACTS devices resulted in a 6.49% reduction in power production costs, a 9.05% reduction in real power losses (RPLs), a 69.5% decrease in voltage deviations (TVDs), and a 26.80% improvement in the voltage stability index (VSI). The approach also achieved a substantial 22.76% reduction in emissions, contributing to environmental sustainability. These results illustrate the robustness of the NSKOA in optimizing power system performance under various operational conditions and its effectiveness in addressing multi-objective problems. It underscores the potential of these approaches in managing the complexities of power systems. It is true that the regulation of control variables, such as generation voltages and transformer tap ratios, ensures the operability of the power system but remains technically insufficient. This highlights the practical implications of RES and FACTS device integration into the design and operation of future power systems, paving the way for a more efficient, reliable, and sustainable electricity supply.