Optimal Power Flow Analysis Based on Hybrid Gradient-Based Optimizer with Moth–Flame Optimization Algorithm Considering Optimal Placement and Sizing of FACTS/Wind Power

Abstract

Optimal power flow (OPF) is one of the most significant electric power network control and management issues. Adding unreliable and intermittent renewable energy sources to the electrical grid increase and complicates the OPF issue, which calls for using modern optimization techniques to solve this issue. This work presents the optimal location and size of some FACTS devices in a hybrid power system containing stochastic wind and traditional thermal power plants considering OPF. The FACTS devices used are thyristor-controlled series compensator (TCSC), thyristor-controlled phase shifter (TCPS), and static var compensator (SVC). This optimal location and size of FACTS devices was determined by introducing a multi-objective function containing reserve costs for overestimation and penalty costs for underestimating intermittent renewable sources besides active power losses. The uncertainty in the wind power output is predicted using Weibull probability density functions. This multi-objective function is optimized using a hybrid technique, gradient-based optimizer (GBO), and moth–flame optimization algorithm (MFO).

1. Introduction

Electrical utilities always strive to find the best solutions for scheduling generation to lower the cost of production and satisfy safe and reliable operating restrictions and power transmission restrictions. Optimal Power Flow (OPF) is a non-linear, complicated optimization technique. The main objective of the OPF problem is to reach the best scheduling of control variables to decrease power losses and fuel cost, taking into account the achievement of the permissible limits in the system [1,2,3]. The OPF issue was developed utilizing a multi-objective function subject to set of constraints. The goal of OPF issue is to piecewise optimize a function with one or several objectives, such as gasoline cost, fuel cost reduction due to the valve-point impact, voltage profile enhancement, decrease of emissions, and stability of voltage, while taking limitations on equality and inequality into account. Controls variables in the OPF issue are transformer tap ratios, active or reactive power of generators, and bus voltages, whereas state variables are reactive power of the output generator and voltage of the load bus [4,5]. In optimization problems, different methods were created to solve OPF. Initially, a variety of mathematical methodologies were used, including interior point [6], linear programming [7], and Newton Raphson [8], to address the OPF problem. However, the results of these approaches were far from ideal, particularly for complex and large systems, and high-dimensional problems are tough to solve. These approaches have some drawbacks, such as a tendency to converge to local optimums and a lengthy computation time. Furthermore, because of their non-linear properties, these approaches are ineffective for issues involving valve-point effects. As a result, heuristic techniques have largely superseded traditional mathematical methods for solving electric power system issues. Heuristic methods for issue resolution can deliver more effective and fast answers when compared to traditional techniques. Heuristic approaches are frequently influenced by natural laws, biological characteristics, or observed collective behavior. Various heuristic methods have recently been developed, such as krill herd algorithm (KHA) [9], moth swarm algorithm (MSA) [10], grey wolf optimization with pattern search (GWO-PS) [11], salp swarm optimizer (SSO) [12], modified rao-2 algorithm [13], a joint self-adaptive particle swarm optimization (SPSO), and differential evolution algorithms [14], and employed in the OPF problem’s solution. In 2002 [15], the author employed the PSO to minimize overall cost, voltage enhancement, as well as voltage stability in IEEE 30-bus system. In 2003 [16], the author used the simulated annealing (SA) approach to reduce fuel costs in IEEE 30-bus systems. In 2004, the GA [17] was used to reduce the cost of gasoline in an IEEE 30-bus system. In 2012, the author used the GSA to minimize fuel cost, voltage profile improvement, and voltage stability in IEEE 30- and 57-bus systems. In 2015 [18], the author used the chaotic krill swarm algorithm (CKSA), while in 2017, different objectives were decreased, such as fuel cost, power loss, and emissions, using moth swarm algorithm (MSA) [10]. Fossil fuels utilized for the production of electricity in electric power networks have now reached the end of their useful life. Furthermore, hydrocarbons are to blame in order to reduce greenhouse gas emissions, which are created by CO₂, NO₂, and SO₂, which contribute to global warming issues. A global majority has agreed to reduce emissions because of rising environmental consciousness and economic concerns [19]. Wind energy is a popular renewable energy source because it is long-lasting, ecologically benign, and cost-effective [20,21]. It has significant benefits in terms of lowering fuel costs and reducing pollutants. However, in this analysis, noise pollution and investment expenses were not taken into account. Wind power has now been included in the answer to the OPF problem by researchers. In 2015 [22], the authors solved the OPF challenge by employing the modified bacteria foraging algorithm (MBFA) to combine wind and fossil fuel generation units [20], while in 2016 [23], the author integrated wind units into an IEEE 14-bus system, and the OPF issues were handled using quadratic programming. Despite all of this development, the wind speed is unpredictably variable is a major issue. The weibull probability distribution function (PDF) is widely employed in long-term distribution problems for wind speed calculation [19]. In the operation of electric power systems, a lack of reactive power produces issues, such as voltage drop, variations in voltage, and voltage collapse, and voltage instability in severe circumstances. There are numerous options for dealing with the aforementioned situations; to mention a few, employing reactive power compensation and load shedding is used, which is less compatible than reactive power injection methods. Furthermore, generators and capacitor banks generate reactive power that is extremely sluggish to adjust to unexpected variations in load level. FACTS devices technology may be an appropriate option in this case. FACTS devices can help electric networks in a variety of ways, including lowering power losses, voltage stability enhancement, and improving the security of the network [24,25]. FACTS devices are extremely beneficial to a power grid; including them in OPF issues increases the difficulty of finding a good solution. To this purpose, the majority of methodologies in the papers to resolve the OPF issue using FACTS devices are based on heuristic and metaheuristic algorithms [26,27], with a quick overview of this set of algorithms given below: hybrid particle swarm optimization, particle swarm optimization (PSO) and adaptive GSA [28], krill herd algorithm (KHA) [29], fuzzy harmony search algorithm (FHSA) [30], bAT search algorithm [31], self-adaptive DE (SADE) [32], and symbiotic organisms search (SOS) algorithm [33]. Many attempts have been made in recent decades to identify effective techniques due to the significant nonlinearity of optimum power flow (OPF), and the inclusion of unpredictable and intermittent renewable sources into the electrical network increases the problem’s complexity; in this paper, the authors suggest machine learning (ML), improved salp swarm algorithm (ISSA), deep belief network (DBN), and combining chaotic particle swarm optimizer and gravitational search algorithm to solve the OPF problem [34,35,36,37], respectively. In practice, power plants face a number of practical limits, two of which include the valve-point impact and various fuels option; however, in order to arrive at a more realistic solution, certain constraints must be considered, and each is discussed further below [38]:

• To begin, power plants can be powered by a number of fuels, including fossil, natural gas, and so on. As a result, it is possible to utilize various cost coefficients for different fuels.
• In thermal power plants, a large number of valves are employed to regulate steam flow and unit output power. It is worth noting that opening the steam-admission valves can cause rapid variations in active power losses. In addition, it adds additional ripples to the cost function of generators as a result of the abrupt increase in active power losses.

However, we can summarize the main contribution of the paper as follows:

• Presenting and employing a novel proposed meta-heuristic methodologies for a transmission system with unexpected wind power and FACTS devices.
• Thermal power and wind power cost models are presented in this paper in detail.
• The models of three facts devices (TCSC, TCPS, and SVC) are presented in this paper in detail.
• This paper proposes a meta-heuristic optimization technique known as hybrid gradient-based optimizer and moth-flame optimization algorithm (GBO-MFO) technique to minimize the generation cost, reduce the power losses, minimize the cost and power losses, and compare with three other techniques (GBO, MFO, SMA, and CFA).
• Four cases are studied in this paper with the following aims: minimizing the generation cost, reducing the power losses, minimizing the cost and power losses, and minimizing the cost and power losses with uncertain load demand.
• Three facts devices are incorporated in the IEEE 30-bus (TCSC, TCPS, and SVC), and the location and rating of three types of FACTS devices are optimized in case studies with objectives of minimizing cost and system real power loss.

In this paper, our study explores an appropriate model of wind power and optimizes the placement of different FACTS devices (SVC, TCSC, and TCPS). The rating and position of these devices are adjusted to reduce the cost of production in addition to minimizing the active power loss. hybrid gradient-based optimizer (GBO) and moth-flame optimization algorithm (MFO) technique are used to carry out the optimization.

The remainder of the paper is laid out as follows: Section 2 includes cost models for both thermal and wind generation. The problem of the placement of the FACTS devices in a wind-powered network is formulated in Section 3. Meanwhile, in Section 4, the objective of the optimization is studied. The proposed GBO-MFO algorithms are introduced in Section 5. Case studies, simulation results, and detailed comparisons are given in Section 6. Section 7 shows the discussion of the results. Finally, the conclusion of this study is presented in Section 8.

2. Thermal Unit Fuel or Generating Costs

For fossil-fuel-burning thermal units, a quadratic cost model is typically used. The thermal unit’s generation cost in dollars per hour ($/h) is presented as follows:

$$C_{T1j=\left(P_{TGj}\right)}=a_{j+ }{b}_j{P}_{TGj}+C_{j }{P}_{TGj}{}^2$$

(1)

where the cost coefficients $a_j$; $b_j$ and $C_{j }$ refer to the thermal generator that produces $P_{TGj}$, which is provided in

Table 1. Thermal power unit cost coefficients [41].

Generator	Location	a ($/h)	b ($/MWh)	c ($/MW2h)	d ($/h)	e (rad/MW)
(TG1)	1	0	2	0.00375	18	0.0375
(TG2)	2	0	1.75	0.0175	16	0.038
(TG8)	8	0	3.25	0.00834	12	0.045
(TG13)	13	0	3	0.025	13.5	0.041

. Several nozzle groups control the steam, which is pumped into the turbines of thermal power plants, and a substantial increase in fuel costs is noticed due to the influence of wire pulling. This is explained by adding ripples to the cost curve, taking into account the valve-point loading effect; the generation cost of thermal power plants is:

$$C_{T2j}\left(P_{TGj}\right)=a_{j+ }{b}_j{P}_{TGj}+C_{j }{P}_{TGj}{}^2+\left|d_{j }\times \mathrm{Sin}\left(e_{j }\times \left(P_{TGj}{}^{Min}-P_{TGj}\right)\right)\right|$$

(2)

where the coefficients $d_{j }$ and $e_{j }$ account for valve point loading as shown in Table 1. The Weibull probability density function (PDF) is used to simulate the wind speed (v in m/s), and it can be defined as follows:

$$f\left(v\right)=\left(\frac{\partial }{\varnothing }\right) {\left(\frac{v}{\varnothing }\right)}^{\left(\partial -1\right)}{e}^{-\left(-v/\varnothing \right)}$$

(3)

where $\varnothing$ and $\partial$ represent the scale and shape parameter of PDF, respectively.

Note that the rating of the turbine is 3 MW. Bus 5 transports the production of a wind farm with 25 turbines, whereas bus 11 transports the amount of energy produced by 20 windmills. The output power of the turbine is calculated as follows:

$$P_{wind}=\left\{\begin{array}{c}0\\ P_r\left(\frac{v-v_{in}}{v_r-v_{in}}\right)\\ P_r\end{array}\right\} \begin{array}{c}for v <v_{in} and v>v_{out}\\ for v_{in}\le v<v_r \\ for v_r\le v<v_{out}\end{array}$$

(4)

where $v_{in}$ represents the turbine’s cut-in wind speed, which is evaluated as 3 m/s; $v_r$ represents the rated wind speed, which is evaluated as 16 m/s; and $v_{out}$ represents the turbine’s cut-out wind speed, which is evaluated at 25 m/s, and $P_r$ is wind turbine’s rated output power. A wind turbine’s output power is discontinuous at specific wind speeds, as shown by (4) [39,40], which calculates the probability of wind power in discrete zones:

$$f_w(P_{wind})\left\{(P_{wind}=0\right\}=1-exp\left[- {\left(\frac{v_{in}}{\varnothing }\right)}^{(\partial )}\right]+exp\left[- {\left(\frac{v_{out}}{\varnothing }\right)}^{(\partial )}\right]$$

(5)

$$f_w(P_{wind})\left\{(P_{wind}=P_r\right\}=exp\left[- {\left(\frac{v_r}{\varnothing }\right)}^{(\partial )}\right]-exp\left[- {\left(\frac{v_{out}}{\varnothing }\right)}^{(\partial )}\right]$$

(6)

In addition, [39,40] calculates the likelihood for the continuous zone between $v_r$ and $v_{in}$ as follow:

$$\begin{array}{cc}\hfill f_w(P_{wind})=& \frac{\partial \left(v_r-v_{in}\right)}{{\varnothing }^{\partial }\ast P_r}{\left[v_{in}+\frac{P_{wind}}{P_r}\left(v_r-v_{in}\right)\right]}^{\partial -1}\hfill \\ & \ast exp\left[-{\left(\frac{v_{in}+\frac{P_{wind}}{P_r}\left(v_r-v_{in}\right)}{\varnothing }\right)}^{\partial }\right]\hfill \end{array}$$

(7)

2.1. Wind Energy Cost Estimation

2.1.1. Wind Energy’s Direct Cost

There is no need for fuel in wind turbines. A wind farm is usually operated by a personal company that sells a certain quantity of energy to the grid operator, often known as the independent system operator (ISO). As a result, if the ISO owns a wind farm, it is possible that there is no direct cost element unless it wants to inflict some punishment for the first configuration, that is, alternatively. This might be classified as a cost of upkeep [40]. If, on the other hand, the wind farm is owned by a private entity, ISO pays a fixed fee that is in relation jointly agreed upon by the ISO and the wind farm owner to the amount of planned power. Therefore, wind energy’s direct costs can be written as

$${\mathrm{C}}_{\mathrm{di}}\left({\mathrm{P}}_{\mathrm{wsci}}\right)={\mathrm{g}}_{\mathrm{wi}}\left({\mathrm{P}}_{\mathrm{sci}}\right)$$

(8)

where ${\mathrm{g}}_{\mathrm{wi}}$ is the coefficient for a wind farm.

2.1.2. Cost Analysis of Unreliable Wind Power

Due to the unpredictable nature of wind power, the real power generated by wind farms can vary. The actual power may be lower than expected. For example, when the potential of wind energy is overstated, to alleviate the requirement, the ISO must maintain a reserve. The expense of setting aside money in case of overvalued wind energy is classified as a reserve expense [39,40]; the reserve cost can be estimated as

$$\begin{array}{cc}\hfill {\mathrm{C}}_{\mathrm{Rewi}}\left({\mathrm{P}}_{\mathrm{wsci}}-{\mathrm{P}}_{\mathrm{wavi}}\right)& ={\mathrm{M}}_{\mathrm{Rewi}}\left({\mathrm{P}}_{\mathrm{wsci}}-{\mathrm{P}}_{\mathrm{wavi}}\right)\hfill \\ & ={\mathrm{M}}_{\mathrm{Rewi}}{{\displaystyle \int }}_0^{{\mathrm{P}}_{\mathrm{wsci}}}\left({\mathrm{P}}_{\mathrm{wsci}}-{\mathrm{P}}_{\mathrm{wind}}\right){\mathrm{f}}_{\mathrm{windi}}({\mathrm{P}}_{\mathrm{wind}}){\mathrm{dp}}_{\mathrm{wind}}\hfill \end{array}$$

(9)

where ${\mathrm{M}}_{\mathrm{Rewi}}$ is assigned to the coefficient of reserve cost, and ${\mathrm{P}}_{\mathrm{wavi}}$ is assigned to real power generated from the plant. The reserve cost, r.h.s., of the previous equation can be increased as follows:

$$\begin{array}{cc}\hfill {\mathrm{M}}_{\mathrm{Rewi}}{{\displaystyle \int }}_0^{{\mathrm{P}}_{\mathrm{wsci}}}& \{\left({\mathrm{P}}_{\mathrm{wsci}}-{\mathrm{P}}_{\mathrm{wind}}\right)\frac{\partial \left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)}{{\varnothing }^{\partial }{\ast \mathrm{P}}_{\mathrm{ri}}}{\left[{\mathrm{v}}_{\mathrm{in}}+\frac{{\mathrm{P}}_{\mathrm{wind}}}{{\mathrm{P}}_{\mathrm{ri}}}\left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)\right]}^{\partial -1}\hfill \\ & \ast \exp \left[-{\left(\frac{{\mathrm{v}}_{\mathrm{in}}+\frac{{\mathrm{P}}_{\mathrm{wind}}}{{\mathrm{P}}_{\mathrm{ri}}}\left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)}{\varnothing }\right)}^{\partial }\right]\}{\mathrm{dp}}_{\mathrm{wind}}\hfill \\ & +{\mathrm{M}}_{\mathrm{Rewi}}\left({\mathrm{P}}_{\mathrm{wsci}}-0\right){\ast \mathrm{f}}_{\mathrm{windi}}({\mathrm{P}}_{\mathrm{wind}})\left\{({\mathrm{P}}_{\mathrm{wind}}=0\right\}\hfill \end{array}$$

(10)

where ${\mathrm{P}}_{\mathrm{ri}}$ is the rated power of the plan. Quite the opposite of overvalued, there are possibilities when the electricity output of wind farms is measured as greater than the legally agreed-upon scheduled power. If the energy cannot be used by decreasing the output power of the thermal generators, it is squandered. The ISO pays a penalty for excess energy.

The penalty cost element for underestimating wind energy is calculated as follows [40]:

$$\begin{array}{cc}\hfill C_{Pewi}\left(P_{wavi}-P_{wsci}\right)& =M_{Pewi}\left(P_{wavi}-P_{wsci}\right)\hfill \\ & =M_{Pewi}{{\displaystyle \int }}_{P_{wsci}}^{P_{ri}}\left(P_{wind}-P_{wsci}\right)f_{windi}(P_{wind})d{p}_{wind}\hfill \end{array}$$

(11)

where $M_{Rewi}$ is assigned as the coefficient of reserve cost; the penalty cost, r.h.s., of the previous equation can be increased as follows:

$$\begin{array}{c}{\mathrm{M}}_{\mathrm{Pewi}}{{\displaystyle \int }}_{{\mathrm{P}}_{\mathrm{wsci}}}^{{\mathrm{P}}_{\mathrm{ri}}}\left\{\left({\mathrm{P}}_{\mathrm{wind}}-{\mathrm{P}}_{\mathrm{wsci}}\right)\frac{\partial \left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)}{{\varnothing }^{\partial }{\ast \mathrm{P}}_{\mathrm{ri}}}{\left[{\mathrm{v}}_{\mathrm{in}}+\frac{{\mathrm{P}}_{\mathrm{wind}}}{{\mathrm{P}}_{\mathrm{ri}}}\left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)\right]}^{\partial -1}\times \exp \left[-{\left(\frac{{\mathrm{v}}_{\mathrm{in}}+\frac{{\mathrm{P}}_{\mathrm{wind}}}{{\mathrm{P}}_{\mathrm{ri}}}\left({\mathrm{v}}_{\mathrm{r}}-{\mathrm{v}}_{\mathrm{in}}\right)}{\varnothing }\right)}^{\partial }\right]\right\}{\mathrm{dp}}_{\mathrm{wind}}\\ +{\mathrm{M}}_{\mathrm{Rewi}}\left({\mathrm{P}}_{\mathrm{ri}}-{\mathrm{P}}_{\mathrm{wsci}}\right)\times {\mathrm{f}}_{\mathrm{windi}}({\mathrm{P}}_{\mathrm{wind}})\left\{({\mathrm{P}}_{\mathrm{wind}}={\mathrm{P}}_{\mathrm{ri}}\right\}\end{array}$$

(12)

3. Modeling of FACTS Devices

In this paper, we use three FACTS devices, including shunts and series; the series compensation devices are the TCSC and TCPS, and they are utilized to improve the line’s loading capacity, whereas the shunt compensation device is SVC, which supplies reactive power to the network.

3.1. Model of Thyristor-Controlled Series Compensator Phase Shifter (TCSC)

TCSC is one of the most effective and extensively utilized FACTS devices. In comparison to standard control devices, this device is mounted in series with a transmission line and provides smooth and flexible line-impedance control with significantly faster response; it has shown to be a good method for managing power flow particularly across lengthy transmission lines, where power angle control at either end is limited. Furthermore, the TCSC can be used to improve the electric network’s transient stability. Figure 1 depicts the TCSC’s basic circuit structure; it comprises of a thyristor-controlled reactor (XL) in parallel with a fixed series capacitor (XC). The TCSC is treated as a variable capacitive reactance by taking into account the reactance ${\mathrm{X}}_{\mathrm{C}}<{\mathrm{X}}_{\mathrm{L}}$. Controlling the firing angle ($\mathsf{\gamma }$) of the thyristors can change the inductive reactance. As a result, the effective reactance of TCSC can be stated as follows:

$${\mathrm{X}}_{\mathrm{TCSC}}\left(\mathsf{\gamma }\right)=\frac{{\mathrm{X}}_{\mathrm{C}}{\mathrm{X}}_{\mathrm{L}}\left(\mathsf{\gamma }\right)}{{\mathrm{X}}_{\mathrm{L}}\left(\mathsf{\gamma }\right)-{\mathrm{X}}_{\mathrm{C}}}$$

(13)

Figure 2 depicts the TCSC static model fitted between buses a and b after adding TCSC. The transmission line’s modified reactance (${\mathrm{X}}_{\mathrm{eq}})$ is written as follows:

$${\mathrm{X}}_{\mathrm{eq}}={\mathrm{X}}_{\mathrm{ab}}-{\mathrm{X}}_{\mathrm{TCSC}}=\left(1-\mathsf{\tau }\right){\mathrm{X}}_{\mathrm{ab}}$$

(14)

The degree of series compensation is referred to as follows:

$$\mathsf{\tau }=\frac{{\mathrm{X}}_{\mathrm{TCSC}}}{{\mathrm{X}}_{\mathrm{ab}}}$$

(15)

where ${\mathrm{X}}_{\mathrm{ab}}$ is the line-inductive reactance. The power flow equations for the transmission line that includes TCSC are expressed as follows:

$${\mathrm{P}}_{\mathrm{ab}}={\mathrm{V}}^2{}_{\mathrm{a}}{\mathrm{G}}_{\mathrm{ab}}-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)$$

(16)

$${\mathrm{Q}}_{\mathrm{ab}}=-{\mathrm{V}}^2{}_{\mathrm{a}}{\mathrm{B}}_{\mathrm{ab}}-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)+{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)$$

(17)

$${\mathrm{P}}_{\mathrm{ba}}={\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)+{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)$$

(18)

$${\mathrm{Q}}_{\mathrm{ba}}=-{\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}+{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)+{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\theta }}_{\mathrm{a}}-{\mathsf{\theta }}_{\mathrm{b}}\right)$$

(19)

$${\mathrm{G}}_{\mathrm{ab}}=\frac{{\mathrm{R}}_{\mathrm{ab}}}{{\mathrm{R}}^2{}_{\mathrm{ab}}+{\left[\left(1-\mathsf{\tau }\right){\mathrm{X}}_{\mathrm{mn}}\right]}^2}$$

(20)

$${\mathrm{B}}_{\mathrm{ab}}=\frac{\left(1-\mathsf{\tau }\right){\mathrm{X}}_{\mathrm{ab}}}{{\mathrm{R}}^2{}_{\mathrm{ab}}+{\left[\left(1-\mathsf{\tau }\right){\mathrm{X}}_{\mathrm{mn}}\right]}^2}$$

(21)

3.2. Model of Thyristor-Controlled Phase Shifter (TCPS)

Figure 3 shows the TCPS model fitted between the transmission lines that connect buses a and b, and the equations for the line’s power flow are as follows:

$${\mathrm{P}}_{\mathrm{ab}}=\frac{{\mathrm{V}}^2{}_{\mathrm{a}}{\mathrm{G}}_{\mathrm{ab}}}{{\mathrm{Cos}}^2{}_{\mathsf{\alpha }}}-\frac{{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}}{\mathrm{Cos}\mathsf{\alpha }}\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)\right]$$

(22)

$${\mathrm{Q}}_{\mathrm{ab}}=-\frac{{\mathrm{V}}^2{}_{\mathrm{a}}{\mathrm{B}}_{\mathrm{ab}}}{{\mathrm{Cos}}^2{}_{\mathsf{\alpha }}}-\frac{{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}}{\mathrm{Cos}\mathsf{\alpha }}\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)\right]$$

(23)

$${\mathrm{P}}_{\mathrm{ba}}=-{\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{G}}_{\mathrm{ab}}-\frac{{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}}{\mathrm{Cos}\mathsf{\alpha }}\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)\right]$$

(24)

$${\mathrm{Q}}_{\mathrm{ba}}=-{\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{B}}_{\mathrm{ab}}-\frac{{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}}{\mathrm{Cos}\mathsf{\alpha }}\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}+\mathsf{\alpha }\right)\right]$$

(25)

At bus a and bus b, TCPS injects actual and reactive power, expressed as follows:

$${\mathrm{P}}_{\mathrm{as}}=-{\mathrm{G}}_{\mathrm{ab}}{\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{Tan}}^2\mathsf{\alpha }-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}\mathrm{Tan}\mathsf{\alpha }\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)-{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)\right]$$

(26)

$${\mathrm{Q}}_{\mathrm{as}}={\mathrm{B}}_{\mathrm{ab}}{\mathrm{V}}^2{}_{\mathrm{b}}{\mathrm{Tan}}^2\mathsf{\alpha }+{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}\mathrm{Tan}\mathsf{\alpha }\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)\right]$$

(27)

$${\mathrm{P}}_{\mathrm{bs}}=-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}\mathrm{Tan}\mathsf{\alpha }\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)+{\mathrm{B}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)\right]$$

(28)

$${\mathrm{Q}}_{\mathrm{bs}}=-{\mathrm{V}}_{\mathrm{a}}{\mathrm{V}}_{\mathrm{b}}\mathrm{Tan}\mathsf{\alpha }\left[{\mathrm{G}}_{\mathrm{ab}}\mathrm{Cos}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)-{\mathrm{B}}_{\mathrm{ab}}\mathrm{Sin}\left({\mathsf{\sigma }}_{\mathrm{a}}-{\mathsf{\sigma }}_{\mathrm{b}}\right)\right]$$

(29)

3.3. The Model of Static var Compensator (SVC)

In electric power systems, SVC is utilized to provide quick operational reactive power control; they control electric variables affecting the power grid quickly using a mix of capacitors and reactors for a variety of reasons. They are used in electric power systems, and one of the most important is quick voltage management at network weak points. SVC is coupled in a thyristor-controlled reactor ($X_{Lcereectivelyd in a Thyristor-}={\omega }_C$) with fixed capacitor $X_{Ccereectivelyd in a Thyristor-}=\frac{1}{{\omega }_C}$ (TCR-FC). Set up in its most basic version, SVC can be used for inductive as well as capacitive correction. Figure 4 and Figure 5 depict the TCR-FC standard diagram and equivalent circuit, respectively.

The equivalent susceptibility is determined by the following formula:

$$B_{SVC}=B_C+B_L\left(\mathsf{\gamma }\right)$$

(30)

where

$$B_C=wC and B_L\left(\mathsf{\gamma }\right)=\frac{1}{wL}\left(1-\frac{2\gamma }{\pi }-\frac{Sin2\gamma }{\pi }\right)$$

(31)

SVC provides reactive power, which may be represented as follows in terms of power flow:

$$Q_{SVC}=-V^2{}_m\ast B_{SVC}$$

(32)

4. Objective of Optimization

In this paper, the researcher proposed the modified electric network IEEE 30-bus as a test case that contains 30 bus; 41 branches [42]; 4 thermal generators at buses 1, 2, 8, and 13; and the 2 thermal generators at buses 5 and 11 are replaced by wind generator. Moreover, three FACTS devices are preferred and placed at an optimal location with appropriate ratings, namely SVC, TCSC, and TCPS.

• Total costs of generation

The fundamental goal of optimum power flow is to be formulated containing all the costs in the modified system, which is illustrated as follows:

$$C_{Tot\left(g\right)}={\displaystyle \sum }_{j=1}^{N_T}{C}_{T2j}\left(P_{TGj}\right)+{\displaystyle \sum }_i^{N_w}\left[C_{di}\left(P_{wsci}\right)+C_{Rewi}\left(P_{wsci}-P_{wavi}\right)+C_{Pewi}\left(P_{wavi}-P_{wsci}\right)\right]$$

(33)

where $N_T$ is the thermal generation number, and $N_w$ is the wind generation number.

• Active power losses

In the OPF issue, the network parameters, such as voltage deviation and active power loss in the transmission line of the system, are very important. The network active power loss is estimated as

$$P_L={\displaystyle \sum }_{d=1}^{Nb}{G}_{d\left(ab\right)}\left[V^2{}_a+V^2{}_b-2V_a{V}_b Cos\left({\sigma }_{ab}\right)\right]$$

(34)

where ${\sigma }_{ab}={\sigma }_a-{\sigma }_b$ is the difference in voltage angles between buses a and b, and ${\sigma }_{ab}$ is transferred conductance of the line between a and b; $Nb$ is the number of transmission lines.

• Voltage deviation

In addition, by decreasing the voltage deviation of all load buses, we improve the voltages at all load busses by considering the constraint control, which will help to make the network more stable and reliable. The voltage deviation is estimated mathematically as follows:

$$VD={\displaystyle \sum }_{a=1}^{NL}(\left|V_{La}-1\right|)$$

(35)

where $NL$ is the number of loaded buses.

4.1. Equality and Inequality Constraints

4.1.1. Operational Equality Constraints

Equality constraints are the power balance equations where the generated active and reactive power in the electrical system is equal to the loads, and power losses in cases without FACTS devices in the network can be defined as follows:

$$P_{ga}-P_{da}-V_a{\displaystyle \sum }_{N=1}^{Nb}{V}_b{\gamma }_{ab}Cos\left({\partial }_{ab}+{\sigma }_a-{\sigma }_b\right)=0$$

(36)

$$Q_{ga}-Q_{da}-V_a{\displaystyle \sum }_{N=1}^{Nb}{V}_b{\gamma }_{ab}Sin\left({\partial }_{ab}+{\sigma }_a-{\sigma }_b\right)=0$$

(37)

where $Nb$ is the number of busses, $P_{ga}$ and $Q_{ga}$ are the generated active and reactive power, while $P_{da}$ and $Q_{da}$ are the active and reactive load at bus a, and ${\partial }_{ab}$ is the difference in the angle of the admittance matrix.

In addition, the power balance in case of integrating FACTS devices can be defined as

$$P_{ga}+P_{as}-P_{da}-V_a{\displaystyle \sum }_{N=1}^{Nb}{V}_b{\gamma }_{ab}Cos\left({\partial }_{ab}+{\sigma }_a-{\sigma }_b\right)=0$$

(38)

$$Q_{ga}+Q_{as}-Q_{da}-Q_{SVCa}-V_a{\displaystyle \sum }_{N=1}^{Nb}{V}_b{\gamma }_{ab}Sin\left({\partial }_{ab}+{\sigma }_a-{\sigma }_b\right)=0$$

(39)

where $P_{as}$ is the generated real power, $Q_{as}$ is the reactive power incorporated by the TCPS, and $Q_{SVCa}$ is the reactive power incorporated by SVC.

4.1.2. Operational Inequality Constraints

The inequality constraints include the three restrictions listed below:

• Generator constraints:

Active, reactive, and voltage of the generator lay between upper and lower boundaries as follows:

$$P_{gj}{}^{min} \le P_{gj}\le P_{gj}{}^{max} \forall j\in N\left(gen\right)$$

(40)

$$Q_{gj}{}^{min} \le Q_{gj}\le Q_{gj}{}^{max} \forall j\in N\left(gen\right)$$

(41)

$${\mathrm{V}}_{\mathrm{gj}}{}^{\min }\le {\mathrm{V}}_{\mathrm{gj}}\le {\mathrm{V}}_{\mathrm{gj}}{}^{\max }\forall \mathrm{j}\in \mathrm{N}\left(\mathrm{gen}\right)$$

(42)

• Transformer tap setting constraints:

Transformer tap setting (T) must be between lower and upper boundaries as shown below:

$$T_s{}^{min} \le T_s\le T_s{}^{max} \forall s\in N\left(Tr\right)$$

(43)

where $N(Tr$) is the number of transformers.

• FACTS devices constraints:

  $$TCSC : {\tau }_{TCSCa}{}^{min} \le {\tau }_{TCSCa}\le {\tau }_{TCSCa}{}^{max} \forall a\in N\left(TCSC\right)$$

  (44)

  $$TCPS : {\tau }_{TCPSa}{}^{min} \le {\tau }_{TCPSa}\le {\tau }_{TCPSa}{}^{max} \forall a\in N\left(TCPS\right)$$

  (45)

  $$SVC : Q_{SVCi}{}^{min} \le Q_{SVCi}\le Q_{SVCi}{}^{max} \forall i\in N\left(SVC\right)$$

  (46)

  where $\mathrm{N}\left(TCSC\right),$$\mathrm{N}\left(TCPS\right)$, and $\mathrm{N}\left(SVC\right)$ are the numbers of FACTS devices.

5. Modified Moth–Flame Optimization and Gradient-Based Optimizer (GBO-MFO)

5.1. Moth–Flame Optimization (MFO) Algorithm

MFO was suggested by Mirjalili in [43]. MFO begins by producing moths at random in the solution space, as shown in the flow chart in Figure 6, and then estimating each moth’s fitness values and labeling the optimal position with flame; the moths’ positions are then updated using a spiral movement function to improve the flame-marked spots and to update the new best individual places and continuing the preceding procedures (the locations of the moths are being updated, and new ones are being created) until the termination requirements are reached.

There are three steps in the MFO technique corresponding to the measurements to take, which are as follows below.

5.1.1. Creating the Initial Moth Population

Each moth, according to Mirjalili in [43], can fly in one of four dimensions; moths can be expressed as a group as follows:

$$S=\left[\begin{array}{cc}\begin{array}{cc}{S}_{1,1}& S_{1,2}\\ S_{2,1}& S_{22}\end{array}& \begin{array}{cc}\cdots & S_{1,m}\\ \cdots & S_{2,\mathrm{m}}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ S_{n,1}& S_{n,2}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & S_{n,m}\end{array}\end{array}\right]$$

(47)

where n denotes the number of moths in the search space, and m denotes the number of dimensions. In addition, all moth fitness values are memorized as

$$OS=\left[\begin{array}{c}O{S}_1\\ O{S}_2\\ \begin{array}{c}⋮\\ O{S}_n\end{array}\end{array}\right]$$

(48)

Flames make up the remaining components of the MFO algorithm. In three-dimensional space, the flames along with their fitness function vector are shown in the following matrix below:

$$A=\left[\begin{array}{cc}\begin{array}{cc}{A}_{1,1}& A_{1,2}\\ A_{2,1}& A_{22}\end{array}& \begin{array}{cc}\cdots & A_{1,m}\\ \cdots & A_{2,m}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ A_{n,1}& A_{n,2}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \dots & A_{n,m}\end{array}\end{array}\right]$$

(49)

$$OA=\left[\begin{array}{c}O{A}_1\\ O{A}_2\\ \begin{array}{c}⋮\\ O{A}_n\end{array}\end{array}\right]$$

(50)

5.1.2. Positions of the Moths Are Being Updated

To find the global best solution for the optimization issues, MFO uses three separate functions. The following are the definitions for these functions:

$$MFO=\left(I,P,T\right)$$

(51)

where I denotes the moths’ initial random locations, P denotes the moths’ mobility in the search space, and T denotes the end of the search process. The I function, which is utilized to implement the random distribution, is represented by using the following formula:

$$S\left(i,j\right)=\left(UP\left(i\right)-LP(j\right)\ast rand\left(i\right))\ast +LP\left(i\right)$$

(52)

where the lower and upper boundaries of variables are denoted by UP and LP, respectively.

The moths use a transverse orientation to fly through the search space, as previously stated. When working with a logarithmic spiral, there are three things to keep in mind:

• The moth should be the initial point on the spiral.
• The flame’s location should be the spiral’s termination point.
• The range of the spiral should not fluctuate beyond the search space.

As a result, the MFO algorithm’s logarithmic spiral is

$$W\left(S_i,A_j\right)=C_i\ast e^{bt}\ast Cos\left(2\pi t\right)+A_j$$

(53)

$C_i$ denotes the distance between the i-th moth and the j-th flame, b provides a fix to determine the logarithmic spiral’s shape, and t is a random number between (−1,1). The spiral motion of the moth towards the flame in the search zone in MFO ensures a balance between exploitation and exploration.

5.1.3. Updating the Number of Flames

This section focuses on enhancing the exploitation of the MFO method (for example, updating the moths’ positions in n distinct locations in the search space may minimize the risk of exploiting the most promising solutions). As a result, reducing the number of flames aids in the resolution of this problem as shown by the equation below:

$$No.flame=round\left(N-l\ast \frac{N-l}{T}\right)$$

(54)

where N represents the maximum number of flames, $l$ represents the current number of iterations, and T represents the maximum number of iterations.

5.2. Gradient-Based Optimizer (GBO)

The suggested GBO [44] combines gradient and population-based approaches, and it utilizes the Newton’s technique, which specifies the search direction to investigate the search domain with the use of a collection of vectors and two primary operators, namely gradient search rule and local escaping operators.

5.2.1. Initialization Process

$$\mathsf{\beta }={\mathsf{\beta }}_{\min }+\left({\mathsf{\beta }}_{\max }-{\mathsf{\beta }}_{\min }\right)\ast {\left(1-{(\frac{\mathrm{h}}{\mathrm{H}})}^3\right)}^2$$

(55)

A set of choice variables, restrictions, and a goal-oriented function are all part of an optimization issue, and the probability rate and the parameters for transition between exploration and exploitation are among the GBO’s control parameters. Each person of the population is known as a “vector” in the suggested algorithm. As a result, in a D-dimensional search space, the GBO contains NB vectors. As a result, a vector may be written in the following way:

$${\mathrm{K}}_{\mathrm{m},\mathrm{d}}=\left[{\mathrm{K}}_{\mathrm{m},1},{\mathrm{K}}_{\mathrm{m},2},\dots {\mathrm{K}}_{\mathrm{m},\mathrm{D}}\right]\mathrm{m}=1,2,\dots \mathrm{N}\mathrm{d}=1,2,\dots \mathrm{D}$$

(56)

$${\mathrm{K}}_{\mathrm{m}}={\mathrm{K}}_{\min }+\mathrm{rand}\left(0,1\right)\ast \left({\mathrm{K}}_{\max }-{\mathrm{K}}_{\min }\right)$$

(57)

where ${\mathrm{K}}_{\max },{\mathrm{K}}_{\min }$ are the decision variable K limits, while rand (0, 1) is a random number in the range (0, 1).

5.2.2. Gradient Search Rule (GSR) Process

A significant factor called $\mathsf{\delta }$ is used in the GBO algorithm to achieve a well-balanced search of a large area while yet attaining near-optimal and global results. The following is how the $\mathsf{\delta }$ is used:

$${\mathsf{\delta }}_1=2\ast \mathrm{rand}\ast \mathsf{\alpha }-\mathsf{\alpha }$$

(58)

$$\mathsf{\alpha }=\left|\mathsf{\beta }\ast \mathrm{Sin}\left(\frac{2\mathsf{\pi }}{2}+\sin \left(\mathsf{\beta }\ast \frac{2\mathsf{\pi }}{2}\right)\right)\right|$$

(59)

where ${\mathsf{\beta }}_{\min }$ is a constant of 0.2, and ${\mathsf{\beta }}_{\max }$ is a constant of 1.2; $\mathrm{h}$ is the current iteration number, and $\mathrm{H}$ is the total number of iterations; the ${\mathsf{\delta }}_1$ parameter is in charge of balancing exploration and exploitation. Within a specific range, the parameter value rises via defined iterations. This broadens the range of possible answers and allows the algorithm to investigate different options. The GSR can be determined in the following way:

$$\mathrm{GSR}={\mathrm{rand}\mathrm{m}\ast \mathsf{\delta }}_1\ast \frac{2\mathsf{\Delta }{\mathrm{k}\ast \mathrm{k}}_{\mathrm{m}}}{\left({\mathrm{k}}_{\mathrm{worst}}-{\mathrm{k}}_{\mathrm{best}}+\mathsf{\sigma }\right)}$$

(60)

In Equation (60), the random offset is defined as the distinction between the best and second-best solution (${\mathrm{k}}_{\mathrm{best}}$) and a randomly chosen solution (${\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}$). The variable’s meaning changes as it is iterated $\mathsf{\Delta }\mathrm{k}$, as seen by the Equation below (63). A second random number (rand n) is also given to enable further research as follows:

$$\mathsf{\Delta }\mathrm{K}=\mathrm{rand}\left(1:\mathrm{N}\right)\ast \left|\mathrm{step}\right|$$

(61)

$$\mathrm{step}=\frac{\left({\mathrm{k}}_{\mathrm{best}}-{\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}\right)+\mathsf{\vartheta }}{2}$$

(62)

$$\mathsf{\vartheta }=2\ast \mathrm{rand}\ast \left(\left|\frac{{\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}+{\mathrm{k}}_{\mathrm{s}2}^{\mathrm{t}}+{\mathrm{k}}_{\mathrm{s}3}^{\mathrm{t}}+{\mathrm{k}}_{\mathrm{s}4}^{\mathrm{t}}}{4}-{\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}\right|\right)$$

(63)

where rand(1: N) is a random vector of N elements in the range of (0, 1). The four integers chosen at random are $\mathrm{s}1,\mathrm{s}2,\mathrm{s}3,\mathrm{and}\mathrm{s}4$ such that $\mathrm{s}1\ne \mathrm{s}2\ne \mathrm{s}3\ne \mathrm{s}4\ne \mathrm{m}$ step represents a phase scale, which is represented by ${\mathrm{k}}_{\mathrm{best}}$ and ${\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}$. The term DM is calculated by taking given a set of potential vectors, the optimal vector, and transferring the current vector (${\mathrm{k}}_{\mathrm{m}}$) in the direction of the best vector (${\mathrm{k}}_{\mathrm{best}}-{\mathrm{k}}_{\mathrm{m}})$ as follows:

$${\mathrm{k}}_{\mathrm{m}+1}={\mathrm{k}}_{\mathrm{m}}-\mathrm{GSR}$$

(64)

$$\mathrm{DM}={\mathrm{rand}\ast \mathsf{\delta }}_2\left({\mathrm{k}}_{\mathrm{best}}-{\mathrm{k}}_{\mathrm{m}}\right)$$

(65)

${\mathsf{\delta }}_2$ is a random parameter that is used to change the size of each vector agent’s phase, and it can be calculated as follows:

$${\mathsf{\delta }}_2=2\ast \mathrm{rand}\ast \mathsf{\alpha }-\mathsf{\alpha }$$

(66)

Finally, DM and GSR can be modified as

$$\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}=\mathrm{k}{1}_{\mathrm{m}}^{\mathrm{t}}-\mathrm{GSR}+\mathrm{DM}$$

(67)

where $\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}$ denotes the modified vector as a result of changing $\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}$, and the transformation $\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}$ can be known as follows according to Equations (55) and (65):

$$\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}={\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}-{\mathrm{rand}\mathrm{m}\ast \mathsf{\delta }}_1\ast \frac{2\mathsf{\Delta }{\mathrm{k}\ast \mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}}{({\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}}-{\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}}+\mathsf{\sigma })}+{\mathrm{rand}\ast \mathsf{\delta }}_2\ast \left({\mathrm{k}}_{\mathrm{best}}-{\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}\right)$$

(68)

where ${\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}},{\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}}$ are identical to ${\mathrm{y}}_{\mathrm{n}}+\mathsf{\Delta }\mathrm{k}$, and is equal to the average of two elements, ${\mathrm{A}}_{\mathrm{m}+1}$,${\mathrm{k}}_{\mathrm{m}}$.

$${\mathrm{A}}_{\mathrm{m}+1}={\mathrm{k}}_{\mathrm{m}}-\mathrm{rand}\mathrm{m}\ast \frac{2\mathsf{\Delta }{\mathrm{k}\ast \mathrm{k}}_{\mathrm{m}}}{({\mathrm{k}}_{\mathrm{worst}}-{\mathrm{k}}_{\mathrm{best}}+\mathsf{\sigma })}$$

(69)

Despite the fact that ${\mathrm{k}}_{\mathrm{m}}$ represents the current solution, we get the present solution element $\mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}$ by using the preceding procedure and replacing the present solution element ${\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}$n with the new solution element ${\mathrm{k}}_{\mathrm{best}}$.

$$\mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}={\mathrm{k}}_{\mathrm{best}}-{\mathrm{rand}\mathrm{m}\ast \mathsf{\delta }}_1\ast \frac{2\mathsf{\Delta }{\mathrm{k}\ast \mathrm{k}}_{\mathrm{m}}^{\mathrm{t}}}{({\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}}-{\mathrm{yP}}_{\mathrm{m}}^{\mathrm{t}}+\mathsf{\sigma })}+{\mathrm{rand}\ast \mathsf{\delta }}_2\ast \left({\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}-{\mathrm{k}}_{\mathrm{s}2}^{\mathrm{t}}\right)$$

(70)

Using Equation (68), the GBO method seeks to increase discovery and exploitation by making them more transportable than the last Equation (70), which is utilized to aid the local search’s exploitation phase. Finally, a revised version of the response available has been discovered, which is as follows:

$${\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}+1}={\mathrm{s}}_1\ast \left({\mathrm{s}}_2\ast \mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}+\left(1-{\mathrm{s}}_2\right)\ast \mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}\right)+\left(1-{\mathrm{s}}_1\right)\ast \mathrm{K}{3}_{\mathrm{m}}^{\mathrm{t}})$$

(71)

where ${\mathrm{s}}_1{\mathrm{and}\mathrm{s}}_2$ are random numbers, and $\mathrm{K}{3}_{\mathrm{m}}^{\mathrm{t}}$ is presented as

$$\mathrm{K}{3}_{\mathrm{m}}^{\mathrm{t}}={\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}+1}-{\mathsf{\delta }}_1\ast (\mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}-\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}})$$

(72)

5.2.3. The Local Escaping Operator (LEO) Process

The LEO is a tool that may be used to improve an optimization’s performance method by assisting in the solution of complex engineering issues. The LEO operator aids the algorithm’s convergence by allowing it to easily move out of local optima positions. The LEO operator sets out to create a new solution that is a more efficient solution (${\mathrm{K}}_{\mathrm{LEO}}^{\mathrm{t}})$ by different solutions (${\mathrm{k}}_{\mathrm{best}}$); $\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}$ solutions are chosen randomly selected from the population, and ${\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}},{\mathrm{k}}_{\mathrm{s}2}^{\mathrm{t}}$ solutions are created at random. It efficiently upgrades current solutions and follows the following procedure:

If rand < P,

$${\mathrm{k}}_{\mathrm{LEO}}^{\mathrm{t}}\{\begin{array}{c}{\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}+1}+{\mathrm{f}}_1\left({\mathrm{h}}_1{\mathrm{k}}_{\mathrm{best}}-{\mathrm{h}}_2{\mathrm{K}}_{\mathrm{e}}^{\mathrm{t}}\right)\\ +{\mathrm{f}}_2{\delta }_1\left({\mathrm{h}}_3\left(\mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}-\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}\right)\right)\\ {\mathrm{h}}_3\left({\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}-{\mathrm{k}}_{\mathrm{s}2}^{\mathrm{t}}\right)/2\text{ if }\mathrm{rand}<0.5\\ {\mathrm{k}}_{\mathrm{m}}^{\mathrm{t}+1}+{\mathrm{f}}_1\left({\mathrm{h}}_1{\mathrm{k}}_{best}-{\mathrm{h}}_2{\mathrm{K}}_{\mathrm{e}}^{\mathrm{t}}\right)\\ +{\mathrm{f}}_2{\text{δ}}_1\left({\mathrm{h}}_3\left(\mathrm{K}{2}_{\mathrm{m}}^{\mathrm{t}}-\mathrm{K}{1}_{\mathrm{m}}^{\mathrm{t}}\right)\right)\\ +{\mathrm{h}}_3\left({\mathrm{k}}_{\mathrm{s}1}^{\mathrm{t}}-{\mathrm{k}}_{\mathrm{s}2}^{\mathrm{t}}\right)/2,\mathrm{otherwise}\end{array}$$

(73)

P is the value of probability, and P = 0.5 and ${\mathrm{h}}_1,{\mathrm{h}}_2,{\mathrm{and}\mathrm{h}}_3$ are random created values as follows:

$${\mathrm{h}}_1=\{\begin{array}{c}2\ast \mathrm{rand} \mathrm{if} {\mathsf{\omega }}_1<0.5\\ 1 \mathrm{otherwise}\end{array}$$

(74)

$${\mathrm{h}}_2=\{\begin{array}{c}\mathrm{rand} \mathrm{if} {\mathsf{\omega }}_1<0.5\\ 1 \mathrm{otherwise}\end{array}$$

(75)

$${\mathrm{h}}_3=\{\begin{array}{c}\mathrm{rand} \mathrm{if} {\mathsf{\omega }}_1<0.5\\ 1 \mathrm{otherwise}\end{array}$$

(76)

where rand is random number between (0, 1). The above values of ${\mathrm{h}}_1,{\mathrm{h}}_2,{\mathrm{and}\mathrm{h}}_3$ can be reorganized as

$${\mathrm{h}}_1={\mathrm{D}}_1\ast 2\ast \mathrm{rand}+\left(1-{\mathrm{D}}_1\right)$$

(77)

$${\mathrm{h}}_2={\mathrm{D}}_1\ast \mathrm{rand}+\left(1-{\mathrm{D}}_1\right)$$

(78)

$${\mathrm{h}}_3={\mathrm{D}}_1\ast \mathrm{rand}+\left(1-{\mathrm{D}}_1\right)$$

(79)

where ${\mathrm{D}}_1$ is a binary vector takes a value of (0, 1); as an example, if parameter $\mathsf{\omega }$ < 0:5, then, the value is ${\mathrm{D}}_1$; otherwise, ${\mathrm{D}}_1=0$, where the answer ${\mathrm{k}}_{\mathrm{e}}^{\mathrm{t}}$ is created as follows:

$${\mathrm{k}}_{\mathrm{e}}^{\mathrm{t}}=\{\begin{array}{c}{\mathrm{k}}_{\mathrm{rand}} \mathrm{if} {\mathsf{\omega }}_2=0.5 \\ {\mathrm{k}}_{\mathrm{P}}^{\mathrm{t}} \mathrm{otherwise}\end{array}$$

(80)

${\mathrm{K}}_{\mathrm{m}}$ is a random answer that is created using the formula below:

$${\mathrm{K}}_{\mathrm{m}}={\mathrm{K}}_{\min }+\mathrm{rand}\left(0,1\right)\ast \left({\mathrm{K}}_{\max }-{\mathrm{K}}_{\min }\right)$$

(81)

5.3. Proposed GBO-MFO Algorithm

The original GBO algorithm suffers from the local optima when handling non-linear objective functions. Additionally, the convergence speed and precision are low because of the randomness. To overcome all these drawbacks and improve the solution efficiency as well as the rate of convergence, we proposed the GBO-MFO technique [45]. Figure 7 presents the flowchart of the suggested GBO-MFO algorithm. As shown in this Figure, the GBO is represented as the fundamental algorithm in the suggested hybrid algorithm. The GSR and LEO are used to update the place in the search space regarding the optimal solution in the original GBO. On the other hand, the original MFO takes use of the spiral movement ability to achieve the exploitation more than exploring the solution space [46]. The proposed algorithm is the hybridization of the GBO and MFO that owns the advantage of both algorithms.

6. Simulation Results and Discussion of the IEEE 30-Bus Test Network

As an example, the IEEE 30-bus electric network was utilized. The modified single-line diagram of IEEE 30-bus is presented in Figure 8. It has four thermal generators, two wind generators, and 24 load nodes. The major data and operational conditions of the test system are investigated in [42]. The electric network consists of four thermal generators at buses 1, 2, 8, and 13; two wind generators at buses 5 and 11; and four tap-changing transformers in branches 11, 12, 15, and 36. The voltage values of the PV nodes are within the boundaries of (0.95–1.1 p.u.). The tap transformers settings are considered within the range of (0.9–1.1 p.u.). In order to verify the strength of the proposed algorithm, GBO-MFO was implemented to solve the OPF issue with TCSC, SCPS, and SVC and practical limitations. The proposed technique is advanced in MATLAB 2014 programming circumference, and all sections of simulations have been coded in a PC with 2.40 GHZ frequency and the installed memory (RAM) is 4.0 GB. In this paper, the power flow is estimated by MATPOWER version 7 [47]. The electrical network loads are 283.4 MW and 126.2 MVar, and the main simulation results are presented in this section as shown in

Table 2. The results of GBO-MFO with fixed loads for the modified IEEE 30-bus system.

Control Variable	Minimum	Maximum	Case A	Case B	Case C
${\mathrm{P}}_{\mathrm{TG}2}\left(\mathrm{MW}\right)$	20	80	40.36711	25.17703	39.96882
${\mathrm{P}}_{\mathrm{WG}5}\left(\mathrm{MW}\right)$	0	75	49.60228	74.99606	74.99998
${\mathrm{P}}_{\mathrm{TG}8}\left(\mathrm{MW}\right)$	10	35	10.00001	34.99992	34.99953
${\mathrm{P}}_{\mathrm{WG}11}\left(\mathrm{MW}\right)$	0	60	42.08583	59.99999	60
${\mathrm{P}}_{\mathrm{TG}13}(\mathrm{MW}$	12	40	12.00002	39.99076	25.3085
${\mathrm{V}}_1\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.077394	1.055951	1.05933
${\mathrm{V}}_2\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.062138	1.050173	1.054368
${\mathrm{V}}_5\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.039988	1.039544	1.045164
${\mathrm{V}}_{11}\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.040376	1.045106	1.0472
${\mathrm{V}}_{12}\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.094419	1.095728	1.096266
${\mathrm{V}}_{13}\left(\mathrm{p}.\mathrm{u}\right)$	0.95	1.10	1.054076	1.080656	1.068085
${\mathrm{T}}_{11}\left(\mathrm{p}.\mathrm{u}\right)$	0.90	1.10	0.998881	1.02985	1.010771
${\mathrm{T}}_{12}\left(\mathrm{p}.\mathrm{u}\right)$	0.90	1.10	0.974578	0.942647	0.989468
${\mathrm{T}}_{15}\left(\mathrm{p}.\mathrm{u}\right)$	0.90	1.10	0.970168	1.011504	1.006795
${\mathrm{T}}_{36}\left(\mathrm{p}.\mathrm{u}\right)$	0.90	1.10	0.98209	0.978315	0.987358
FACTS ratings
${\mathsf{\tau }}_{\mathrm{TCSC}1}(\%)$	0	50%	0.261798	0.499378	0.260031
${\mathsf{\tau }}_{\mathrm{TCSC}2}(\%)$	0	50%	0.256903	0.180707	0.499985
${\varnothing }_{\mathrm{TCSC}1}\left(\mathrm{deg}\right)$	−5	5	2.891543	4.606817	2.906457
${\varnothing }_{\mathrm{TCSC}2}\left(\mathrm{deg}\right)$	−5	5	−0.07969	−2.38928	−0.99951
${\mathrm{Q}}_{\mathrm{SVC}1}\left(\mathrm{MVAr}\right)$	−10	10	9.998593	5.468421	9.441854
${\mathrm{Q}}_{\mathrm{SVC}2}\left(\mathrm{MVAr}\right)$	−10	10	9.967787	9.980042	9.999946
FACTS locations
TCSC1 Branch	1	40	5	34	2
TCSC2 Branch	1	41	2	41	9
TCPS1 Branch	1	40	14	35	33
TCPS2 Branch	1	41	39	14	5
SVC1 Bus	3	29	7	19	24
SVC2 Bus	3	30	24	24	21
Parameters
${\mathrm{P}}_{\mathrm{TG}1}\left(\mathrm{MW}\right)$	50	200	134.9094224	50.0016021	49.9999977
${\mathrm{Q}}_{\mathrm{TG}1}\left(\mathrm{MVAr}\right)$	20 -	150	5.7579908	−3.37055433	−2.3974411
${\mathrm{Q}}_{\mathrm{TG}2}\left(\mathrm{MVAr}\right)$	20 -	60	19.0386182	9.216845035	9.67332152
${\mathrm{Q}}_{\mathrm{WG}5}\left(\mathrm{MVAr}\right)$	30 -	35	19.5321975	21.38442996	22.2010301
${\mathrm{Q}}_{\mathrm{TG}8}\left(\mathrm{MVAr}\right)$	−15	48.7	34.73924937	31.42417937	31.0560877
${\mathrm{Q}}_{\mathrm{WG}11}\left(\mathrm{MVAr}\right)$	−25	30	24.97840417	27.42443353	27.6472454
${\mathrm{Q}}_{\mathrm{TG}13}\left(\mathrm{MVAr}\right)$	−15	44.7	9.24650880	25.99634823	17.2964094
Objective function
${\mathrm{C}}_{\mathrm{gen}}\left(\$/\mathrm{h}\right)$			807.120060	939.3285458	916.651707
${\mathrm{P}}_{\mathrm{loss}}\left(\mathrm{Mw}\right)$			5.56466263	1.76537072	1.87682289
${\mathrm{C}}_{\mathrm{gross}}\left(\$/\mathrm{h}\right)$			1363.586323	1115.86561	1104.33400
$\mathrm{VD}\left(\mathrm{p}.\mathrm{u}.\right)$			0.6622278	0.862291816	0.83473561
Emission			0.213568466	0.141615940	0.14196922

. The obtained results of the proposed hybrid GBO-MFO algorithm are compared with the results of four recent algorithms, including the gradient-based optimizer (GBO) [44], moth–flame optimization algorithm [43], slime mould algorithm (SMA) [48], and coulomb–franklin’s algorithm (CFA) [49]. The contol parameters of the proposed and others algorithms are illustrated in details in

Table 3. The control parameters of the proposed and other techniques.

Algorithms	The Control Parameter
Common parameters	Number of population size = 200
	Iterations number = 500
	Dimensions number = 27
	Number of runs = 20
GBO-MFO	b = 1 and a decreases linearly from −1 to −2 (Default), pr = 0.5
MFO	b = 1 and a decreases linearly from −1 to −2 (Default)
GBO	pr = 0.5
CFA	Pc and the value equal to 0.5
SMA	(vb) is a parameter with a range of (−a, a) and gradually approaches zero as the iterations increase. The value of $\overrightarrow{vc}$ oscillates between (−1, 1) and tends to zero eventually.

. This comparison with other modern meta-heuristic techniques is presented in

Table 4. Simple comparison between the proposed and other techniques for the modified IEEE 30-bus system.

Technique		SMA	CFA	GBO-MFO	GBO	MFO
Case A	${\mathrm{C}}_{\mathrm{gen}}\left(\$/\mathrm{h}\right)$	807.277	807.4699193	807.12006	807.2502	807.4733
	${\mathrm{P}}_{\mathrm{loss}}\left(\mathrm{Mw}\right)$	5.5798	5.628556222	5.56466263	5.6002	5.6304
	${\mathrm{C}}_{\mathrm{gross}}\left(\$/\mathrm{h}\right)$	1.370 × 103	1370.325541	1363.586323	1367.3	1.370 × 103
	$\mathrm{VD}\left(\mathrm{p}.\mathrm{u}.\right)$	0.8747	0.75741542	0.6622278	0.8514	0.6534
	Emission	0.2136	0.21360518	0.213568466	0.2136	0.2136
Case B	${\mathrm{C}}_{\mathrm{gen}}\left(\$/\mathrm{h}\right)$	936.6358	939.30095	939.3285458	938.649	1.120 × 103
	${\mathrm{P}}_{\mathrm{loss}}\left(\mathrm{Mw}\right)$	1.8149	1.796401	1.76537072	1.7717	1.8102
	${\mathrm{C}}_{\mathrm{gross}}\left(\$/\mathrm{h}\right)$	1.120 × 103	1118.941	1115.86561	1115.8	1.120 × 103
	$\mathrm{VD}\left(\mathrm{p}.\mathrm{u}.\right)$	0.8391	0.82594	0.862291816	0.8927	0.8551
	Emission	0.1414	0.14157	0.141615939	0.1416	0.1416
Case C	${\mathrm{C}}_{\mathrm{gen}}\left(\$/\mathrm{h}\right)$	918.7776	916.62101	916.651707	918.4229	918.0247
	${\mathrm{P}}_{\mathrm{loss}}\left(\mathrm{Mw}\right)$	1.8602	1.90008	1.876822891	1.8608	1.9179
	${\mathrm{C}}_{\mathrm{gross}}\left(\$/\mathrm{h}\right)$	1104.8	1106.65064	1104.334003	1104.5	1109.8
	$\mathrm{VD}\left(\mathrm{p}.\mathrm{u}.\right)$	0.9109	0.85175	0.834735606	0.9166	0.8610
	Emission	0.1416	0.1416933	0.141969216	0.1417	0.1418

. In the our optimization technique, the number of control variables are 27 as defined in detail in Table 4, which includs the locations and rating for each FACTS device, and while integrating the FACTS devices, some criteria must be considered, such as (1) only one FACTS device is put in the location, (2) SVC is not added on the PV bus, and (3) in case of line having tap transformer, TCSC and TCPS are not integrated into this lines. TCSC can compensate up to 50% of the branch reactance where it is placed, while SVC can provide or absorb up to 10 MVAr of reactive power. Phase shifter’s angle (TCPS) can be changed from −5 and 5. The case studies carried out under this section assume that the electric network is fully loaded; for the given full-loaded situation, three case studies are carried out, which are illustrated as follows:

• Case A: Minimizing the cost generation.
• Case B: Minimizing power losses.
• Case C: Minimizing cost and power losses.
• Case D: Load demand uncertainty

Case study findings using the GBO-MFO method are given in this part along with explanations. In each case study of optimization, each scenario was performed 30 times to determine the goal function’s most optimal value and the accompanying control variable settings.

• Case A: Minimizing the Cost Generation

In Case A, it was observed that the electrical network total generation cost is greatly decreased to 807.12 ($/h) using GBO-MFO technique, which is less than the generation cost in Case B and Case C. Additionally, in this case, the minimum value of the voltage deviation was obtained. Because wind generators are less expensive, their scheduled power outputs are frequently higher than those of thermal units. However, because adequate wind energy is not available, a wind units’ output cannot be used to its most tremendous potential, as this will increase the reserve cost to provide scheduled power for extended periods of time.

The placements of FACTS devices as well as their ratings were adjusted in this article to reduce generating costs. In a genuine system or a variant of a well-known system, FACTS devices are installed without any research, which may result in wasteful and ineffective electric system resource consumption. In this case, buses 7 and 24 are considered to be the best sites for the SVCs, and TCSC was installed in branch 5 and branch 2, and TCPS was installed in branch 14 and branch 39.

In this paper, in all cases, we ran the technique multiple times after adjusting the optimal places to improve ratings while reducing costs. Due to the system’s heavy inductive load, SVC’s reactive power correction was nearly at its maximum. FACTS devices were placed in networks with the goal of improving loading capacity. The electrical network’s actual power loss was kept as low as possible in such a way as to optimize the capacity of the system’s loads. Furthermore, the suggested algorithm’s capacity was used to keep bus voltage amplitudes within a reasonable range. The state variables in the IEEE 30-bus test system ere kept within the allowed limits.

• Case B: Minimizing power losses

The major goal of Case B was to lower active power loss to the absolute minimum. Many transmission lines in power networks result in significant active power loss. As a result, for large-scale systems, minimizing active power loss is critical. The GBO-MFO in Case B had the lowest active power loss.

In this case, the power losses achieved a minimal loss of 1.7467 MW using the control variables listed in Table 2. This value of power losses is better than Case A and Case C, and also the voltage deviation obtained in this case was 0.8622, and the emissions from traditional thermal power plants reached to 0.1416159. In this case, the SVCs were placed at buses 19 and 24, while the TCSC were placed at branches 34 and 41, and the two TCPS were installed at branches 35 and 14 given that reactive power in power networks is constantly distributed based on voltage magnitude differences between buses, resulting in a rise in power loss. This is a critical achievement, and it is also where the suggested algorithm’s effect is scrutinized. As an example, one of the most visible aims of the OPF issue is to keep the numerical value of voltage values of the buses between their borders and around typical operating circumstances.

Constraints on load bus voltage are particularly important in the OPF issue since the voltages at which load buses operate are frequently found to be near to their limitations, which must be kept between 0.95 and 1.05 p.u. in our investigation. Furthermore, the limits on reactive power of generators are also met in all situations.

• Case C: Minimizing cost and power losses

Through the comparison presented in Table 2, we found that the cost of generation in Case A is less than the cost of generation in Case B, while the reduction in power losses in Case B is better than in Case A. This statistic emphasizes the need to have a goal that includes both of them, which forced us to resort to the use of Case C. Developing a cost model in which the loss is translated into an equal price of energy is a straightforward method to account for both objectives. The cost of electricity used in this research is $0.10 per kWh. The cost of power was $0.10 per kWh in this research as presented in Equation (82).

$${\mathrm{C}}_{\mathrm{Gross}}=C_{Tot\left(g\right)}+{\mathrm{P}}_{\mathrm{L}}\ast {10}^3\ast 0.10)$$

(82)

In Case C, the technique’s minimal goal value is 1104.0771 $/h, and the active power loss is 1.87682 MW, total generation cost is 916.65170 $/h, voltage deviation is 0.8347 p.u., and the emission is 0.14196. The two SVCs were installed at busses 24 and 21 and the two TCSC at branches 2 and 9, while the two TCPS WERE placed at branches 33 and 5. In this case, the best cost and loss values WERE determined by the coefficient of price chosen for generators, including thermal and wind. However, taking into account both goals results in the lowest total cost.

• Case D: Load demand uncertainty

In OPF research studies, determination of FACTS devices placement and rating considering uncertainties in load demands is a well-developed research topic. The research scenario discussed in this section considers a genuine incident of changing load. The probability density function (PDF) [45] WAS used to depict the unpredictablity of demands. It Is worth noting that this technique also optimizes the placement of the FACTS, and also changing the placement of FACTS during various loading circumstances is not possible as a result, so the FACTS were carefully located for a particular scenario with the greatest chance of occurring. The probability of Scenario 3 and Scenario 4 occurring is the greatest. We improved the positions of the FACTS devices for Scenario 3, and its same optimal placement was used to optimize the FACTS for other scenarios. In addition, the values of the FACTS devices used were adjusted to account for various system scenarios; the research is useful since it optimizes a variety of electrical system aims by executing OPF at regular intervals. The normal distribution with mean ${\mathsf{\mu }}_{\mathrm{d}}$ and standard deviation ${\mathsf{\alpha }}_{\mathrm{d}}$ is commonly used to model the load demand uncertainty and are 70 and 10, respectively. In this work, all loads are assumed to obtain the same PF, mean, and standard deviation. As a result, rather than using a normal distribution to every load separately, ${\mathsf{\Delta }}_{\mathrm{S},\mathrm{j}}$ shows the probability of each load level that is computed by Equation (83); Equation (84) can also be used to calculate the related value of every load level, which is indicated by $\overline{{\mathrm{P}}_{\mathrm{d},\mathrm{i}}}$.

These formulas [50] are used to compute the mean load and probability of an event for a given case.

$${\mathsf{\Delta }}_{S,j}=\underset{P_{d,j}^{Min}}{\overset{P_{d,j}^{Max}}{{\displaystyle \int }}}\frac{1}{{\alpha }_d\surd 2\pi } \exp [-\frac{{(P_d-{\mu }_d)}^2}{2{\alpha }_{d^2}}]d{P}_d$$

(83)

where $P_{d,j}^{Max}$ and $P_{d,j}^{Min}$ are the maximum and minimum boundaries of ith scenarios of loading, and the mean of ith scenarios of load demand is

$$\overline{P_{d,i}}=\frac{1}{{\mathsf{\Delta }}_{S,j}}\underset{P_{d,j}^{Min}}{\overset{P_{d,j}^{Max}}{{\displaystyle \int }}}\left(P_d\ast \frac{1}{{\alpha }_d\sqrt{2\pi }}\exp \left[-\frac{{(P_d-{\mu }_d)}^2}{2{\alpha }_{d^2}}\right]d{P}_d\right)$$

(84)

Additionally, for the premeditated medians and probability for any and all possible load-carry situations, the percent load in a situation means that the nodes needs are multiplied by the % of that scenario. As a result, the technique optimized the performance index of gross cost as written in (82) under various conditions, and the active power generated from all generators was used at full capacity in every scenario.

The variable’ limits, such as the limits presented in Table 2 and Table 4, present effective system parameters for the four scenarios, which show the best amplitude of gross cost. For all scenarios, the expected gross cost is determined as follows:

$$\mathrm{EGRC}={\displaystyle \sum }_{\mathrm{SC}=1}^{{\mathrm{NO}}_{\mathrm{SC}}}{\mathsf{\Delta }}_{\mathrm{SC}}{\ast \mathrm{C}}_{\mathrm{gross},\mathrm{sc}}$$

(85)

where ${\mathrm{NO}}_{\mathrm{SC}}$ is the scenario’s number. ${\mathsf{\Delta }}_{\mathrm{SC}}$ is the scenario’s probability, which is estimated in

Table 5. Scenarios for loading and probabilities of each.

Loading Scenario	$\%\mathbf{Loading},\overline{P_{d,i}}$	$\mathbf{Scenario}\mathbf{Probability},{\mathsf{\Delta }}_{SC}$
SCA	54.749	0.15866
SCB	0.15866	0.34134
SCC	74.599	0.34134
SCD	85.251	0.15866

. EGC is the expected generation cost, and EPL is expected power loss, which is estimated in [51] and computed by

$$EGC={\displaystyle \sum }_{SC=1}^{N{O}_{SC}}{\mathsf{\Delta }}_{SC}\ast C_{Gen, sc}$$

(86)

$$EPL={\displaystyle \sum }_{SC=1}^{N{O}_{SC}}{\mathsf{\Delta }}_{SC}\ast P_{Loss,sc}$$

(87)

7. Discussion

The numerical results show that, when compared to previous techniques in the literature, such as GBO, MFO, SMA, and CFA, the suggested methodology provides a better and more acceptable solution for handling various versions of the OPF issue while taking into account realistic limitations and FACTS devices. The proposed technique GBO-MFO became stuck in local optima quite fast, and it may converge to a superior solution in less iterations as presented in Figure 9 for Case A, Figure 10 for Case B, and Figure 11 for Case C, unlike the other techniques investigated in this study. For more knowledge, consider Case 1, a well-known case study in the field of OPF problem-solving. The most talented generation cost achieved by the suggested method when FACTS devices are taken into account is 807.12 ($/h). These numbers are lower than those of the other techniques evaluated and those found in the paper, which are 807.2502 ($/h) in GBO and 807.469 in CFA, 807.277 in SMA, and 807.4733 in MFO. It is worth noting that a successful optimization method not only converges to a reduced total generation cost while solving the OPF issue, but it also reduces active power losses. Additionally, the power loss is 5.56 MW in the proposed method, which is better than 5.60, 5.62, 5.57, and 5.6304 in GBO, CFA, SMA, and MFO, respectively.

Figure 12 shows a bar chart depicting the breakdown of different costs; because of the tiny cost coefficient for penalties, the cost of not using obtainable wind energy is minimal in all circumstances. Generator units linked to buses 5, 8, and 11 are necessary in Cases B and C to decrease power losses because system loading in the area is extremely high; therefore, it should run at or near its rated capacity. The reserve cost in Cases B and C increases as the planned power from the wind turbines increases. Because the direct cost is related to its anticipated output, it rises as the scheduled power grows. As a result of the reduced scheduled power of the thermal units, the cost of these units decreases in Cases B and C compared to Case A. The cost of loss is calculated using the energy unit cost. When the power loss is included as a goal in the optimization, the cost is properly reduced.

Likewise, in the Cases B and C, the ideal settings of the generator units’ active power are shown below in Figure 13 for Cases A to C. Figure 14 shows the voltage of load buses of the case studies for the modified system when it is fully loaded. For Cases A to C, the best settings of voltages on the generator bus are presented in Figure 15, and Figure 16 also explains the setting of the transformer tap ratio for Case A to Case C This demonstrates the effectiveness of the suggested approach, which is capable of achieving the best global value with fewer rounds.

Table 6. Wilcoxon signed rank test results.

Algortihms	Case 1				Case 2				Case 3
	R+	R−	p-Value	H0	R+	R−	p-Value	H0	R+	R−	p-Value	H0
GBO-MFO vs. GBO	141	69	1.79 × 10−1	Yes	74	136	2.47 × 10−1	Yes	133	77	2.96 × 10−1	Yes
GBO-MFO vs. MFO	179	31	5.73 × 10−3	No	190	20	1.51 × 10−3	No	206	4	1.63 × 10−4	No
GBO-MFO vs. CFA	166	44	2.28 × 10−2	No	159	51	4.38 × 10−2	No	183	27	3.59 × 10−3	No
GBO-MFO vs. SMA	186	24	2.49 × 10−3	No	183	27	3.59 × 10−3	No	92	118	6.27 × 10−1	Yes

contains the obtained results of using the Wilcoxon signed rank test. The proposed GBO-MFO algorithm is compared with four well-known algorithms in the literature, including GBO, MFO, CFA, SMA, respectively, for the Cases A, B, and C.

In case of loading scenarios, the system costs and active power losses are significantly reduced, as presented in

Table 7. The results of MFO-GBO with different scenarios of loads in Case D.

Control Variable	Scenario A	Scenario B	Scenario C	Scenario D
PTG 2 (MW)	20	20	20.22	22.756
PWG 5 (MW)	36.9422	52.126	62.24	73.61
PTG 8 (MW)	10	10.9434	20.83	28.769
PWG 11 (MW)	27.19	41.31	47.29	52.04
PTG 13 (MW)	12	12	12	15.81
V1 (p. u)	1.057	1.056	1.0594	1.058
V2 (p. u)	1.052	1.0507	1.0539	1.053
V5 (p. u)	1.044	1.044	1.0470	1.047
V11 (p. u)	1.045	1.042	1.0473	1.047
V12 (p. u)	1.0629	1.07	1.0641	1.095
V13 (p. u)	1.061168	1.062836	1.065	1.067
T11 (p. u)	1.042204	0.990643	1.0085	1.055
T12 (p. u)	0.949792	0.969312	0.9837	0.908
T15 (p. u)	0.992094	1.000278	1.0099	0.9944
T36 (p. u)	0.987191	0.984462	0.98495	0.986
FACTS ratings
${\mathsf{\tau }}_{\mathrm{TCSC}1}(\%)$	0.219086	0.219927	0.136841	0.844894
${\mathsf{\tau }}_{\mathrm{TCSC}2}(\%)$	0.063736	0.222954	0.446708	2.262013
${\varnothing }_{\mathrm{TCSC}1}\left(\mathrm{deg}\right)$	1.174577	1.763584	2.946295	0.132565
${\varnothing }_{\mathrm{TCSC}2}\left(\mathrm{deg}\right)$	−4.85577	1.061893	0.170417	0.485028
${\mathrm{Q}}_{\mathrm{SVC}1}\left(\mathrm{MVAr}\right)$	2.837804	7.110695	9.916534	7.510003
${\mathrm{Q}}_{\mathrm{SVC}2}\left(\mathrm{MVAr}\right)$	8.046143	6.707114	6.721043	9.988798
FACTS locations
TCSC1 Branch	14	14	14	11
TCSC2 Branch	23	23	23	33
TCPS1 Branch	24	24	24	22
TCPS2 Branch	41	41	41	39
SVC1 Bus	21	21	21	7
SVC2 Bus	24	24	24	24
Parameters
PTG 1 (MW)	50	50	50	50
QTG 1 (MVAr)	−2.57	−2.04	−1.66	−2.04
QTG 2 (MVAr)	2.76	5.51	7.38	5.65
QWG 5 (MVAr)	10.76	11.40	15.95	15.52
QTG 8 (MVAr)	14.22	19.46	25.62	24.52
QWG 11 (MVAr)	12.68	12.56	9.34	29.89
QTG 13 (MVAr)	8.56	11.96	12.37	14.62
Objective function
${\mathrm{C}}_{\mathrm{gen}}\left(\$/\mathrm{h}\right)$	418.1275	520.6703	624.6884	749.3056
${\mathrm{P}}_{\mathrm{loss}}\left(\mathrm{Mw}\right)$	0.9748	1.0649	1.1673	1.3674
${\mathrm{C}}_{\mathrm{gross}}\left(\$/\mathrm{h}\right)$	515.6028	627.1640	741.4172	886.0478 886.0478
$\mathrm{VD}\left(\mathrm{p}.\mathrm{u}.\right)$	1.0190	0.9441	0.9822	0.9400
Emission	0.1547	0.1 544	0.1525	0.1491

, because of their permanent placements, which are optimized depending on Scenario 3, the FACTS compensate in a major way across all scenarios. In addition, as shown in Case C, the optimal planned power in buses 5, 8, and 11 rise as system loading rises. Furthermore, with the variation in loads, the power of generator units is nearly unaffected. The requirement to reduce to reduce system active power loss is the primary reason behind this. Figure 17 presents the convergence characteristics for different loading scenarios in Case D, and Figure 18 depicts the voltages for each scenario and the upper and lower limits.

8. Conclusions

This study provides a technique to achieve optimal power flow. In the IEEE 30-bus network, there is an issue involving stochastic wind and FACTS devices (OPF), and the suggested MFO-GBO was successfully applied to solve the OPF problem in electric power network. In the three case studies, the placement and rating of the FACTS devices (SVC, TCSC, and TCPS) were optimized with the goal of decreasing generation cost and active power loss in the electric network system. The necessity of considering both factors together demonstrated a composite target that includes generation cost and active power loss. The electric network power loss is translated to an energy cost as well as the expense of production. The MFO-GBO was proven to provide a rapid, accurate, and optimal solution to the OPF issue. The MFO-GBO was contrasted and analyzed for calculating the best set of control variables on IEEE 30-bus to examine the efficacy and superiority of the proposed algorithms to keep generation costs down, decrease power losses to minimum value, and enhance the voltage profile. In addition, MFO-GBO was provided to solve the OPF issue considering uncertain load demand. The collected findings demonstrate the proposed approach’s advantage in obtaining the best solution with the shortest time to convergence. These findings support its efficacy in resolving large-scale OPF issues. The suggested technique outperforms existing in terms of well-known optimization techniques computing performance.

9. Future Recommendation

For future OPF study, the authors recommend incorporating tiny hydro-generators and storage in the form of battery or hydro power in a large system. FACTS devices can be used to create a realistic model of doubly-fed induction generators for wind generators.