Towards Avoiding Cascading Failures in Transmission Expansion Planning of Modern Active Power Systems Using Hybrid Snake-Sine Cosine Optimization Algorithm

Abstract

In this paper, a transmission expansion planning (TEP) model is proposed to guarantee the resilience of power systems and mitigate cascading failures’ impacts. The energy storage systems and fault current limiters’ planning models are integrated into the TEP problem to minimize cascading outages and comply with short-circuit current reliability constraints. Most studies in the literature adopt a single strategy to simulate power systems’ cascading failures that may not be enough to guarantee networks’ resilience. This work elaborates on two scenarios for initiating cascading failures to study the impact of various initiating events on the planned system’s strength and the projects required. The TEP problem is formulated as a non-linear, non-convex large-scale problem. To avoid linearization issues and enhance meta-heuristics performance, a hybridization of two meta-heuristic techniques, namely snake optimizer and sine cosine algorithm (SO-SCA), is proposed to solve the problem. Two hybridization strategies are suggested to improve the exploration and exploitation stages. Defining future loads growth is essential for TEP. Hence, a load forecasting technique based on SO-SCA is investigated and compared with some methods reported in the literature. The results obtained proved the efficiency of the proposed approach in predicting load growth. TEP’s calculations were carried out on the Garver and the IEEE 24-bus system. The results demonstrated the superiority of the hybrid SO-SCA in solving the TEP problem. Moreover, the projects required to expand networks differed according to the type of cascading failures’ initiating scenario.

1. Introduction

The transmission expansion planning (TEP) problem corresponds to allocating different power facilities, such as transmission lines, energy storage systems (ESSs), thyristor-controlled series compensators (TCSC), and fault current limiters (FCLs), to satisfy the expected load growth and fulfill operational and reliability constraints [1]. TEP is a well-studied problem in the literature; however, there are still many challenges related to solution methods and offering a mathematical model that can accommodate the power system’s uncertainties and ensure reliability constraints [2].

Several mathematical, heuristic and meta-heuristic methods were reported in the literature to solve TEP models [2]. Some of these works are introduced below as follows. In Hamidpour et al. [3], the AC power flow-based model (ACTEP) was adopted. The non-linear and non-convex constraints were linearized through Taylor’s expansion and polygon approaches. The CPLEX solver solved the linearized model within GAMS. The suggested strategy was conducted on the IEEE 6-bus and 30-bus test systems. The numerical results verified the efficiency of the proposed solution. Esmaili et al. [4] presented a Benders’ decomposition approach that included a successive linearization technique to solve a multi-stage TEP problem and deal with non-linear constraints. The integration of thyristor TCSCs and FCLs’ allocation problem with the TEP problem increased the problem’s complexity and non-linearity. The proposed method was tested on the IEEE 39- and 118-bus systems. Solvers like GUROBI and CONOPT were used to solve the proposed model. The suggested strategy demonstrated its superiority in terms of linearization error and execution time measured. Franken et al. [5] developed an approach based on the DC power flow-based(DCTEP) model for the simultaneous incorporation of TEP and planning problems of TCSCs, phase-shifting transformers, and high voltage direct current systems (HVDC). The problem was translated into a mixed-integer linear programming model (MILP), and the solver CPLEX within the MATLAB platform was used to solve the problem. In Moradi-Sepahvand et al. [6], a MILP was investigated to formulate TEP integrated HVDC, HVACs, and ESSs in the presence of renewable energy sources. The problem was divided into a master problem and sub-problems and was solved using the CPLEX solver in GAMS. To validate the suggested approach, different test cases were carried out on the Garver network, IEEE 24-bus system, and IEEE 118-bus system.

Regarding the implementation of meta-heuristic techniques, sine cosine algorithm (SCA), Levy flight optimizer (LFO), and the linear population size reduction—success-history-based differential evolution with semi-parameter adaptation hybrid—covariance matrix adaptation evolution strategy algorithm (LSHADE-SPACMA) were tested in Refaat et al. [7] to solve the TEP problem and FCLs allocation problem. DCTEP was formulated, and SCA proved its superiority in obtaining high-quality solutions with a reasonable execution time. In Abdi et al. [8], different algorithms, such as a genetic algorithm, an orthogonal crossover-based differential evolution, grey wolf optimizer, moth–flame optimization, exchange market algorithm, SCA, and imperialistic competitive algorithm, were examined to solve the TEP problem. A comparative analysis was done using the 24 and 118 benchmark systems. The numerical results of different scenarios showed that the grey wolf optimizer was recommended for solving the TEP problem. The non-dominated sorting genetic algorithm II-based Pareto approach was investigated by Abbasi et al. [9] to solve a multi-objective DCTEP problem. The obtained results over two test systems proved the efficiency of the applied technique in dealing with non-linearity in the adopted model.

Meta-heuristic algorithms perform differently in different types of problems. One may perform slightly better than the other in a particular problem and worse in other sets of problems [2]. Further, the success of meta-heuristics mainly hinges on the encoding of candidate solutions and thus the search space. Nowadays, many studies encourage the hybridization of several techniques in one algorithm to overcome these issues. Refaat et al. [7] suggested an approach that combined mathematical and meta-heuristic techniques to solve TEP. A stochastic planning model that considered short-circuit constraints and N-1 security was adopted. The problem was translated into a three-level problem. In the first level, a meta-heuristic algorithm was applied to place new circuits and FCLs such that short-circuit current constraint was achieved. Two mathematical techniques were employed in the second and third stages to decide the best location for new generation units to meet power flow constraints during steady-state operation and define the load shedding to satisfy N-1 reliability during contingencies. The results proved the suggested strategy for handling non-linear constraints and improving the performance of meta-heuristic techniques. Ramachandran et al. [10] proposed an algorithm combining the modified grasshopper optimizer and the improved Harris hawks optimization algorithm to balance the global search and global convergence phases. The results ascertained the capability of the proposed algorithm in providing promising results compared with some algorithms presented in the literature. The hybrid Harris hawks optimizer-arithmetic optimization algorithm was proposed by Cetinbas et al. [11] for autonomous micro-grids design and sizing. The results proved the efficiency of the algorithm suggested. A hosting capacity planning model was developed by Almalaq et al. [12] to increase the use of RESs. A hybrid of two meta-heuristic techniques was employed to deal with non-convexities within the search space. The results showed the effectiveness of the suggested approach in enhancing the exploration and exploitation stages and getting high-quality solutions.

A heuristic approach based on fuzzy logic was developed by Sousa et al. [13] for TEP. The simulation results carried out on different test systems confirmed the proposed technique’s effectiveness in finding high-quality solutions. In Armaghani et al. [14], an analytical multi-stage TEP approach combined with a forward pseudo-dynamic planning method was proposed for finding critical transmission lines that may initiate and propagate the cascading failure.

Recently, cascading failures’ impacts during TEP have been taken into account in several studies [15,16,17]. Mehrtash et al. [15] contributed to a risk-based planning model that minimizes cascading failures’ effects. The suggested problem was modeled as a mixed-integer second-order cone programming problem. Risk indices were considered in the developed model to measure the non-identical probability and individual contingencies’ consequences. The results obtained from the IEEE 24-bus system demonstrated that the system planned was more robust against high-risk contingencies that may have resulted in cascading failures. Qorbani et al. [16] developed a resilient TEP model to mitigate cascading outages’ effects and ensure N-1 security. The outage of two-connected lines to each node was the triggering event for cascading failures. Gjorgiev et al. [17] proposed a risk-informed TEP method. An AC-based cascading failure model was developed to simulate post contingencies power systems’ responses. A multi-objective optimization was investigated and solved using a meta-heuristic algorithm. The suggested procedure succeeded in mitigating the effects of cascading failures.

Calculation of future peak loads is essential for TEP; hence the accuracy of the load forecasting tool implemented is necessary to guarantee the accuracy of the planning process. The hyper-parameters of the forecasting models play an essential role in improving the model’s performance [18,19,20,21]. A long-term load prediction-based particle swarm optimization was implemented by AlRashidi et al. [21]. The results demonstrated the superiority of the applied optimizer in estimating forecasting models’ parameters compared to the least error square method. An adaptive neural fuzzy system was recommended by Fathy et al. [22] and Refaat et al. [23] as an efficient tool for load forecasting. However, it had a problem selecting the type and the number of membership functions. Hong et al. [24] applied the support vector machines method to forecast future loads growth; however, it has parameters selection problems. In this work, a load forecasting approach based on a hybrid of two-meta heuristic techniques is proposed and compared with some methods presented in the literature.

Based on all the above discussion, the existing research gaps are illustrated as follows:

• The mathematical methods efficiently solve linear and simple TEP problems. However, they need a high computational burden for solving large-scale non–linear TEP problems with a vast search space. These methods have difficulty finding the optimal solution, mainly if non-convexities exist with the search space. Most of the existing works apply a linearization scheme and a decomposing strategy to simplify the problem, but it may result in linearization errors that affect the accuracy of obtained solutions.
• As known, the meta-heuristics are not transparent algorithms. They are time-consuming and challenging to reach the optimal solution in a single run [2]. Much effort is still needed for improving the performance of meta-heuristic techniques in solving TEP problems.
• Although heuristic methods need a lower computational burden compared to mathematical and meta-heuristic algorithms, they cannot guarantee the optimal solution or high-quality solutions.
• The hybridization of several optimization algorithms is recommended for solving the TEP problems. However, most current studies proposed it for solving the DCTEP problem that is simple compared to the AC models. Testing the performance of the hybrid algorithms in solving the ACTEP problems is rarely discussed in the literature.
• Most existing works simulate a single strategy to avoid cascading failures of power systems that may not be enough to ensure the electrical networks’ resilience and security.
• Most current research adopts a single optimization algorithm to calculate future loads in the long term. However, the hybrid implementation of several meta-heuristic algorithms is rarely examined in the literature.

This work corresponds to introducing a robust TEP model and investigating a long-term load forecasting method to address the mentioned shortcomings. The ACTEP is adopted in this work. It handles short-circuit current constraints and cascading failures impacts. The proposed model is non-linear and non-convex. To avoid linearization, a hybrid of two meta-heuristic algorithms is developed. The main contributions of this study can be outlined as follows:

• The meta-heuristic algorithm, Snake optimizer (SO) [25], is combined with the sine cosine algorithm (SCA) [26] to solve the proposed TEP problem and calculate the long-term load forecasting. Two combination approaches are adopted to enhance both algorithms’ exploration and exploitation phases.
• The use of renewable energy sources (RESs) threatens the modern power systems’ resilience and may be a trigger to cascading failures and rolling blackouts [27,28]. Therefore, the adopted model deals with RESs’ stochastic behavior.
• Two scenarios for initiating cascading failures are suggested to ensure the planned system’s strength and study the impact of various initiating events on the projects required.
• The impact of ESSs and FCLs in minimizing cascading outages and fulfilling short-circuit current constraints is also investigated. In this regard, Sodium sulfur batteries (NaS) and FCLs’ sizing and placing problems are included in the TEP problem.

The remainder of this paper is systematized as follows. The problem’s formulation is presented in Section 2, while the hybrid of SO-SCA is explained in Section 3. Section 4 shows and discusses the results obtained. Finally, in Section 5, the main remarks derived from this work are presented.

2. Problem Formulation

The TEP problem is modeled as an optimization problem to expand power networks, increase networks’ resilience and achieve reliability constraints in the presence of RESs. The proposed objective function targets to minimize the cost of installed circuits, FCLs, generation units and ESSs.

The AC model is adopted in this study to formulate the proposed problem. It is superior to the DC model in terms of the accuracy of the results obtained and the ability to integrate many problems into a single model, such as combining the problem of reactive power units’ allocation and the TEP into a single problem [27]. A clustered data of wind speed and load consumed over the year simulate power systems’ stochastic behavior.

As mentioned before, TEP is always linked to future loads; therefore, implementing an accurate load forecasting technique is essential to improve the planning process. In this regard, the load forecasting problem is formulated as an estimation process, and the hybrid of SO-SCA is applied to minimize the error associated with the estimated model.

2.1. Mathematical Model of TEP

The proposed TEP model has a single objective function and is subjected to normal and security constraints, as illustrated in (1)–(29). Equation (1) provides the proposed cost function. The first term concerns the cost of installed circuits, while the second and third terms cover the cost of FCLs and NaS batteries required to achieve security and short-circuit current constraints. The fourth term covers the cost of new generation units installed. The reminder terms present the operating cost of NaS batteries and thermal units and RESs, respectively.

$$\begin{array}{l}OF = min\{{\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1,j=1,i\ne j}^{i = n,j=n,i\ne j}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_{ij}(N_{ij}{}^h-N_{ij}{}^{h-1})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+ {\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1,j=1,i\ne j}^{i = n,j=n,i\ne j}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_{FCL}(Z_{ij}{}^{FCL,h}-Z_{ij}{}^{FCL,h-1})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{i = n}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_{ESS}^{inv}(P_{ESS,i}{}^h-P_{ESS,i}{}^{h-1})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{i = n}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_G^{inv}(P_{G,i}{}^h-P_{G,i}{}^{h-1})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{i = n}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_{ESS}^{op}{E}_{ESS,i}{}^h\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{h=1}^{h_{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{i = n}}\frac{\text{λ}{\left(1+\text{λ}\right)}^y}{{\left(1+\text{λ}\right)}^y-1}{C}_G^{op}{P}_{G,i}{}^h\}\end{array}$$

(1)

where the recovery factor $\left(\frac{\mathsf{\lambda }{\left(1+\mathsf{\lambda }\right)}^y}{{\left(1+\mathsf{\lambda }\right)}^y-1}\right)$ is inserted in the cost function to consider the money’s time value [29], where $\mathsf{\lambda }$ and $y$ are the discount rate and the lifetime of the project in years, respectively. The problem constraints are provided in (2)–(29). The steady-state constraints are shown in (2)–(24), whereas reliability constraints are given in (25)–(29). Constraints (2) and (3) are related to active and reactive powers injected at each node. $P_{ESS,i}^h$ is negative if an ESS is in the discharging mode, and it is positive if an ESS is in the charging mode.

$$P_{g,i}^h+P_{R,i}^h-P_{d,i}^h-P_{ESS,i}^h=P_i^h$$

(2)

$$Q_{g,i}^h-Q_{d,i}^h=Q_i^h$$

(3)

So that $P_i^h$ and $Q_i^h{}_i$ are calculated as follows [30]:

$$P_i^h=V_i{}^h{\displaystyle {\displaystyle \sum }_{j=1}^N}{V}_j{}^h\left[G_{ij}cos\left({\theta }_i{}^h-{\theta }_j{}^h\right)+B_{ij}sin\left({\theta }_i{}^h-{\theta }_j{}^h\right)\right]$$

(4)

$$Q_i^h=V_i{}^h{\displaystyle {\displaystyle \sum }_{j=1}^N}{V}_j{}^h\left[G_{ij}sin\left({\theta }_i{}^h-{\theta }_j{}^h\right)-B_{ij}cos\left({\theta }_i{}^h-{\theta }_j{}^h\right)\right]$$

(5)

Constraints (6) and (7) ensure that active and reactive powers generated from thermal stations do not violate their limits. For RESs, it is assumed that they generate active power only, and the output power must be within the permissible bounds, as shown in (8). Constraints (9) and (10) guarantee that the buses’ voltage and angle do not surpass their limits.

$$P_{g,i}^{min}\le P_{g,i}^h\le P_{g,i}^{max}$$

(6)

$$Q_{g,i}^{min}\le Q_{g,i}^h\le Q_{g,i}^{max}$$

(7)

$$P_{R,i}^{min}\le P_{R,i}^h\le P_{R,i}^{max,h}$$

(8)

$${\theta }_i^{min}\le {\theta }_i^h\le {\theta }_i^{max}$$

(9)

$$V_i^{min}\le V_i^h\le V_i^{max}$$

(10)

Equation (11) explains that any circuit built in any previous scenario will be present in the following scenarios. It also ensures that the number of circuits between any two nodes should be less than the maximum. Equations (12) and (13) show that the power flows in new and existing circuits do not exceed the thermal limits.

$$N_{ij}^{h-1}\le N_{ij}{}^h\le N_{ij}^{max}$$

(11)

$$S_{ij}^{s,h}\le N_{ij}{}^h S_{ij}^{max}$$

(12)

$$S_{ij}^{r,h}\le N_{ij}{}^h S_{ij}^{max}$$

(13)

where $S_{ij}^{s,h}$ and $S_{ij}^{r,h}$ are calculated as follows:

$$S_{ij}^{s,h}=\sqrt{{\left(P_{ij}{}^{s,h}\right)}^2+{\left(Q_{ij}{}^{s,h}\right)}^2}$$

(14)

$$S_{ij}^{r,h}=\sqrt{{\left(P_{ij}{}^{s,h}\right)}^2+{\left(Q_{ij}{}^{r,h}\right)}^2}$$

(15)

The other parameters are calculated as follows [30]:

$$P_{ij}{}^{s,h}=V_i^{h^2}{G}_{ij}-V_i^h{V}_j^h\left[G_{ij}cos\left({\theta }_i^h-{\theta }_j^h\right)+B_{ij}sin\left({\theta }_i^h-{\theta }_j^h\right)\right]$$

(16)

$$Q_{ij}{}^{r,h}=-V_i^{h^2}\left(B_{ij}+B_{sh}\right)-V_i^h{V}_j^h\left[G_{ij}sin\left({\theta }_i^h-{\theta }_j^h\right)-B_{ij}cos\left({\theta }_i^h-{\theta }_j^h\right)\right]$$

(17)

$$P_{ij}{}^{r,h}=V_i^{h^2}{G}_{ij}-V_i^h{V}_j^h\left[G_{ij}cos\left({\theta }_i^h-{\theta }_j^h\right)-B_{ij}sin\left({\theta }_i^h-{\theta }_j^h\right)\right]$$

(18)

$$Q_{ij}{}^{r,h}=-V_i^{h^2}\left(B_{ij}+B_{sh}\right)+V_i^h{V}_j^h\left[G_{ij}sin\left({\theta }_i^h-{\theta }_j^h\right)+B_{ij}cos\left({\theta }_i^h-{\theta }_j^h\right)\right]$$

(19)

Constraints related to ESSs’ planning are given in (20)–(24). Equation (20) restricts ESSs’ number at each bus, whilst (21)–(23) maintain ESSs’ power and capacity during charging and discharging modes to be lower than nominal values. ESSs’ storage capacity at each hour is calculated using (24).

$$N_{ESS,i}^{h-1}\le N_{ESS,i}^h\le N_{ESS,i}^{max}$$

(20)

$$0\le P_{dch,i}^h\le P_{ESS,i}^{max}$$

(21)

$$0\le P_{ch,i}^S\le P_{ch,i}^{max,S}$$

(22)

$$0\le E_{ESS,i}^S\le E_{ESS,i}^{max}$$

(23)

$$E_{ESS,i}^h=E_{ESS,i}^{h-1}+{\eta }_{ESS}^{ch}{P}_{ch,i}^S-\raisebox{1ex}{$P_{dch,i}^S$}\!\left/ \!\raisebox{-1ex}{${\eta }_{bat}^{dch}$}\right.$$

(24)

Short-circuit current constraints are provided in (25) and (26). While (27) ensures that FCLs’ size required is lower than the maximum value [4,7].

$$0\le I_i^{SC,s}\le I_{max}^{SC}$$

(25)

$$I_i^{SC}=\frac{V_i\left(0\right)}{Z_{ii}}$$

(26)

$$N_{ij}{}^{s-1} Z_{ij}^{FCL,s-1}\le Z_{ij}^{FCL,s}\le N_{ij}{}^{new,s}{Z}_{ij,max}^{FCL}$$

(27)

Two security scanning matrices, L1 and L2, are suggested as triggering events for the initiation of cascading failures. L1 simulates a power system under a single circuit’s outage between two nodes. In comparison, L2 regards the outage of all circuits between any two nodes as the triggering event for cascading failures. $N_n$ represents total circuits installed at route number n. Systems planned must be able to fulfill the constraints (2)–(27) after contingencies L1 and L2.

$$\mathrm{L}1={\left[\begin{array}{cc}\begin{array}{ccc}{N}_1-1& N_2& \cdots \\ N_1& N_2-1& \cdots \\ ⋮& ⋮& \ddots \end{array}& \begin{array}{ccc}\cdots & N_{n-1}& N_n\\ \cdots & N_{n-1}& N_n\\ ⋮& ⋮& ⋮\end{array}\\ \begin{array}{ccc} ⋮& ⋮& ⋮\\ N_1& N_2& \cdots \\ N_1& N_2& \cdots \end{array}& \begin{array}{ccc}\ddots & ⋮& ⋮\\ \cdots & N_{n-1}-1& N_n\\ \cdots & N_{n-1}& N_n-1\end{array}\end{array}\right]}^{n\times n}$$

(28)

$$\mathrm{L}2={\left[\begin{array}{cc}\begin{array}{ccc}0& N_2& \cdots \\ N_1& 0& \cdots \\ ⋮& ⋮& \ddots \end{array}& \begin{array}{ccc}\cdots & N_{n-1}& N_n\\ \cdots & N_{n-1}& N_n\\ ⋮& ⋮& ⋮\end{array}\\ \begin{array}{ccc}⋮& ⋮& ⋮\\ N_1& N_2& \cdots \\ N_1& N_2& \cdots \end{array}& \begin{array}{ccc}\ddots & ⋮& ⋮\\ \cdots & 0& N_n\\ \cdots & N_{n-1}& 0\end{array}\end{array}\right]}^{n\times n}$$

(29)

2.2. Long-Term Load Forecasting Model

The proposed load forecasting problem is modeled as an estimation process, and the hybrid of SO-SCA is applied to find the optimal coefficients to reduce estimation errors. Linear, quadratic models are commonly used in the literature for forecasting load growth in the long term, and they are adopted in this work.

The historical data used in the forecasting process is a matrix (D) that comprises m rows and two columns. The first column represents the input data, whilst the second column represents the output data. This work uses the input data of historical years (y), and peak loads (P) are the output data. For m rows, the discrete state-space form represents the forecasting system is formulated as follows [21]:

$$P\left(y\right)=H\left(y\right) \phi +\epsilon \left(y\right)$$

(30)

where $P\left(y\right)$ is a vector ($m\times 1)$ that represents calculated outputs, while $\phi$ is a vector ($n\times 1)$ that represents estimated parameters. $H\left(y\right)$ is a row vector representing the relation between $P\left(y\right)$ and $\phi$. $\epsilon \left(y\right)$ is the error vector results in the estimation process.

In this study, the linear and the quadratic models are adopted. For the linear model:

$$H\left(y\right)=\left[\begin{array}{cc}Y& 1\end{array}\right], Y=1,2,\dots .m \mathrm{and} \phi ={\left[\begin{array}{cc}\alpha & \beta \end{array}\right]}^T$$

(31)

while for the quadratic models:

$$H\left(y\right)=\left[\begin{array}{ccc}Y{ }^2& Y& 1\end{array}\right], Y=1,2,\dots .m \mathrm{and} \phi ={\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]}^T$$

(32)

3. The Hybrid Snake-Sine & Cosine-Based Optimization Algorithm

The proposed problem is a large-scale non-linear optimization problem, and to solve it, a hybrid of two meta-heuristic techniques, namely SO-SCA, is developed. SO is recommended for solving several complex and challenging engineering problems. It is capable of converging fast to the global optimum. While, SCA is superior in solving unimodal, multi-modal, and composite test functions compared to other algorithms. SO and SCA have some issues finding optimal or high-quality solutions in one run. Moreover, it is challenging to decide on suitable design parameters for all problems because well-defined algorithm parameters for one situation may not be ideal for solving another problem.

The hybrid SO-SCA is employed in this work to overcome these issues. Two combination approaches are proposed to enhance the exploration and exploitation stages and get high-quality solutions in as few runs as possible.

3.1. Snake Optimizer

SO was developed by Hashim et al. [25]. It is based on the idea of imitating snakes’ mating behavior. Snakes commonly fight to get the best partners if they have enough food and a low ambient temperature. Like any meta-heuristic technique, SO initializes with a randomly generated population, and then the updating process is carried out in two phases: exploration and exploitation. During the updating process, it is assumed that the number of males equals the number of females. This process goes on until the maximum number of iterations (T) is reached.

In the exploration phase, it is assumed that there is not enough food, and snakes randomly move toward searching for food. The food’s quantity (Q) is calculated by $0.5 exp\left(\frac{t-T}{T}\right)$.The exploration’s behavior is modeled as follows:

$$X_{i,m}^{t+1}=X_{rand,m}^t\pm C_2\times A_m\times \left(\left(X_{max}-X_{min}\right)\times rand+X_{min}\right)$$

(33)

$$X_{i,f}^{t+1}=X_{rand,f}^t\pm C_2\times A_f\times \left(\left(X_{max}-X_{min}\right)\times rand+X_{min}\right)$$

(34)

where $X_{i,m}^{t+1}$ and $X_{i,f}^{t+1}$ indicate to the ith male and female positions. $X_{rand,m}^t$ and $X_{rand,f}^t$ are the positions of random males and females selected from the population. $t$ is the current iteration. $rand$ is a random value and $C_2$ is a constant that set to 0.05. $A_m$ and $A_f$ measure respectively the male and female abilities to get food, and they are calculated by $exp\left(\frac{-f_{rand,m}}{f_{i,m}}\right)$ and $exp\left(\frac{-f_{rand,f}}{f_{i,f}}\right)$; respectively. $f_{rand,m}$ and $f_{rand,f}$ are the fitness values of $X_{rand,m}^t$ and $X_{rand,f}^t$; respectively. While $f_{i,m}$ and $f_{i,f}$ are the fitness values of ith male and female.

In the exploitation phase, it is assumed that the food is available, and the next positions are defined according to the ambient temperature ($TEMP)$. $TEMP$ is defined by $exp\left(\frac{-t}{T}\right)$. If the temperature is high and exceeds the threshold value, the new position for males and females is located using (35). Where $X_{food}^t$ is the best individuals position and $C_3$ equals 2.

$$X_{i,f,m}^{t+1}=X_{food}^t\pm C_3\times TEMP\times rand\times \left(X_{food}-X_{i,f,m}^t\right)$$

(35)

Otherwise, the snakes are in the fighting or mating modes and are randomly switched between both modes. The new positions are calculated using (36) and (37) for the fighting mode and using (38) and (39) for the mating mode.

$$X_{i,m}^{t+1}=X_{i,m}^t\pm C_3\times FM\times rand\times \left(X_{best,f}-X_{i,m}^t\right)$$

(36)

$$X_{i,f}^{t+1}=X_{i,f}^{t+1}\pm C_3\times FF\times rand\times \left(X_{best,m}-X_{i,f}^{t+1}\right)$$

(37)

$$X_{i,m}^{t+1}=X_{i,m}^t\pm C_3\times M_m\times rand\times \left(Q\times X_{i,f}^t-X_{i,m}^t\right)$$

(38)

$$X_{i,f}^{t+1}=X_{i,f}^t\pm C_3\times M_f\times rand\times \left(Q\times X_{i,m}^t-X_{i,f}^t\right)$$

(39)

where $FM$ and $FF$ are males’ and females’ fighting abilities and calculated by $exp\left(\frac{-f_{best,f}}{f_i}\right)$ and $exp\left(\frac{-f_{best,m}}{f_i}\right)$; respectively. While $M_m$ and $M_f$ are the mating abilities of males and females fighting and calculated by $exp\left(\frac{-f_{i,f}}{f_{i,m}}\right)$ and $exp\left(\frac{-f_{i,m}}{f_{i,f}}\right)$; respectively.

3.2. Sine Cosine Algorithm

SCA was proposed by Mirjalili et al. [26]. It is based on the use of sine and cosine functions to vary new positions outward and toward the best positions in the hope of enhancing the exploitation and explorations phases. SCA starts with a randomly generated population, and then the updating scheme given in (40) is employed to define the new position of each individual in the population. Where $r_2$, $r_3$ and $r_4$ are randomly defined, and $r_1$ is calculated by $C-\frac{C t}{T}$. $C$ is constant.

$$X_i^{t+1}=\{\begin{array}{c} X_i^t+r_{1}sin\left(r_2\right)\left|r_3\times X_{best,i}^t-X_i^t\right|, if r_4<0.5\\ X_i^t+r_{1}cos\left(r_2\right)\left|r_3\times X_{best,i}^t-X_i^t\right|, if r_4\ge 0.5\end{array}$$

(40)

3.3. Proposed Hybrid Scheme

Two combination approaches for SO-SCA are proposed. The pseudo-code of the first approach is presented in Algorithm 1. SO and SCA operate in parallel to calculate new variables in this approach. According to Algorithm 1, candidate solutions are randomly generated, and their corresponding finesses are calculated. Individuals who violate the steady-state and security constraints are penalized by adding additional costs. After that, the best solution and best fitness value are defined. In the beginning, the population is divided into two groups that have the same number of individuals. New solutions for the first group are calculated using SO‘s updating scheme provided in (33)–(39). While new positions of individuals in the second group are obtained using SCA’s updating scheme given in (40). The best solution for both groups is decided. If the SO gives the best solution, more individuals will be assigned for SO in the next iteration, relying on the switching rate ($N^{rate}$). Otherwise, more individuals will be assigned for SCA.

Algorithm 1: Pseudo-code of SO-SCA algorithm (the first approach)

1: Define the problem’s dimension, population size ($N_p$), stop criterion (T), and upper and lower bounds. 2: Set SO and SCA parameters. 3: Define switching rate ($N^{rate}$).

4: Initialize the population (X). 5: Calculate the cost function and check the problem’s constraints. 6: Add a penalty cost to solutions that violate the constraints. 7: Define the best and the worst solutions. 8: Define SO’s population size ($N_p^{SO}$) = 0.5 $N_p$. 9: Define SCA’s population size ($N_p^{SCA}$) = 0.5 $N_p$.

10: Fort = 2:T 11:   Assign SO’s population ($X^{SO}$) = $X\left(1:N_p^{SO}\right)$. 12:   Assign SCA’s population ($X^{SCA}$) = $X\left(1:N_p^{SCA}\right)$. 13:   For $i=1: N_p^{SO}$ 14:      Define the new position of each individual ($X_{i,SO}^{t+1}$) in $X^{SO}$ using (33)–(39). 15:      Calculate the fitness value of each new solution ($f_i^{SO}$). 16:      Check the problem constraints and add a penalty cost if the new solution violates the constraints. 17:   end for 18:   For$i=1: N_p^{SCA}$ 19:      Define the new position of each individual ($X_{i,SCA}^{t+1}$) in $X^{SCA}$ using (40). 20:      Calculate the fitness value of each new solution ($f_i^{SCA}$). 21:      Check the problem constraints and add a penalty cost if the new solution violates the constraints. 22:   end for 23:   Define the minimum of in the vector $f^{SO}$ ($f_{min}^{SO}$) and its corresponding solution ($X_{Best}^{SO}$). 24:   Define the minimum of in the vector $f^{SCA}$ ($f_{min}^{SCA}$) and its corresponding solution ($X_{Best}^{SCA}$). 25:   If $f_{min}^{SO}\le f_{min}^{SCA}$ 26:     $N_p^{SO}=N^{rate}\ast N_p$ 27:     $N_p^{SCA}=\left(1-N^{rate}\right)\ast N_p$ 28:     $f_{min}=f_{min}^{SO}$ 29:     $X_{min}=X_{Best}^{SO}$ 30:    else 31:     $N_p^{SCA}=N^{rate}\ast N_p$ 32:     $N_p^{SO}=\left(1-N^{rate}\right)\ast N_p$ 33:     $f_{min}=f_{min}^{SCA}$ 34:     $X_{min}=X_{Best}^{SCA}$ 35:   end if 36:   If $f_{min}\le f_{Best}^t$ 37:     $f_{Best}^{t+1}=f_{min}$ 38:     $X_{Best}^{t+1}=X_{min}$ 39:    else 40:     $f_{Best}^{t+1}=f_{Best}^t$ 41:     $X_{Best}^{t+1}=X_{Best}^t$ 42:   end if 43: end for

In the second approach, a series combination is applied, as given in Algorithm 2. The population is first updated using SCA’s updating scheme. Then, the position of the best solution and its fitness is updated regarding the output of SCA. After that, the SO updating scheme is applied to calculate the position of new variables, and the best solution’s position is defined based on the output of SO.

Algorithm 2: Pseudo-code of SO-SCA algorithm (the second approach)

1: Define the problem’s dimension, population size ($N_p$), stop criterion (T), and upper and lower bounds. 2: Set SO and SCA parameters.

3: Initialize the population (X). 4: Calculate the cost function and check the problem’s constraints. 5: Add a penalty cost to solutions that violate the constraints. 6: Define the best and the worst solutions.

7: Fort = 2:T 8:   For$i=1: N_p$ 9:      Define the new position of each individual ($X_{i,SCA}^{t+1}$) in $X^{SCA}$ using (40). 10:    Calculate the fitness value of each new solution ($f_i^{SCA}$). 11:    Check the problem constraints and add a penalty cost if the new solution violates the constraints. 12:   end for 13:   Find the best and the worst solutions obtained using SCA ($X_{Best}^{SCA}$,$f_{Best}^{SCA}$). 14:   Update the positions of the candidate solution ($X_i^t=X_{i,SCA}^{t+1}$); 15:   Update the best solution and the best fitness. 16:    If $f_{Best}^{SCA}\le f_{Best}^t$ 17:        $f_{Best}^t=f_{Best}^{SCA}$ 18:        $X_{Best}^t=X_{Best}^{SCA}$ 19:    else 20:        $f_{Best}^t=f_{Best}^t$ 21:        $X_{Best}^t=X_{Best}^t$ 22:    end if 23:   For i = 1:$N_p$ 24:    Define the new position of each individual ($X_{i,SO}^{t+1}$) in $X^{SO}$ using (33)–(39). 25:    Calculate the fitness value of each new solution ($f_i^{SO}$). 26:    Check the problem constraints and add a penalty cost if the new solution violates the constraints. 27:   end for 28:   Find the best solution and the best fitness. 29:   If $f_{Best}^{SO}\le f_{Best}^t$ 30:     $f_{Best}^{t+1}=f_{Best}^{SO}$ 31:     $X_{Best}^{t+1}=X_{Best}^{SO}$ 32:   else 33:     $f_{Best}^{t+1}=f_{Best}^t$ 34:     $X_{Best}^{t+1}=X_{Best}^t$ 35:   end if 36: end for

4. Numerical Results

4.1. Test Systems

The proposed model and solution method were examined on the Garver network and the IEEE 24-bus system. A data set introduced in [21] was used to investigate the adopted load forecasting algorithm. The wind speed and load consumed depicted in Figure 1 were used to validate the proposed TEP model. Investment and operating factors of different facilities used are given in

Table 1. Investment and operating factors of components used.

Component	Capital Cost	Operating Cost
		Fixed (×106 $/MW)	Variable (×106 $/MWh)
Wind unit [31]	0.139 × 106 $/MW	0.232	0
Thermal units [31]	0.536 × 106 $/MW	0.123	0.144
FCL	0.5 × 106 $/p.u	NA	NA
NaS [29]	0.446 × 106 $/MW	NA	0.2981

. The results were obtained using the MATLAB r2017a platform.

The Garver network is a benchmark system and is widely used for investigating the efficiency of TEP’s models. It has six buses, three-generation units, and five load buses [22], [32]. Bus 5 was suggested as a candidate position for new wind units, and buses 2, 5, and 6 were suggested for new ESSs. The network’s initial configuration is shown in Figure 2. New candidate routes are presented in dash lines.

The IEEE-24 bus system has 24 buses, eleven generation units, and seventeen load buses. Some modifications to the system were made in [32] to make congestions. The candidate locations for new wind units were buses 3, 10, and 19. While the candidate locations for ESSs were buses 2, 6, 8, 14, 20, 21, 22, and 24. The initial configuration of the system is depicted in Figure 3.

4.2. Testing the Performance of SO-SCA

Figure 4, Figure 5, Figure 6 and Figure 7 show the convergence curves of the hybrid algorithm using the two approaches, SO and SCA, over ten runs. Noteworthy, all the simulations in this section were carried out without embedding L1 and L2 securities. The population size was 30 for SO-SCA #2 and 40 for SO-SCA #1, SO, and SCA. It is shown that the series approach effectively solved the Garver network in terms of standard deviation and quality of solutions compared to the parallel approach, as shown in

Table 2. Simulation results of SO-SCA #1, SO-SCA #2, SO, and SCA for the Garver network.

Optimizer	Best Cost (×106 $)	Worst Cost (×106 $)	Mean Value (×106 $)	Standard Deviation	Execution Time (S)
SO-SCA 1 *	58.27	66.31	60.64	2.99	227.15
SO-SCA 2 *	58.27	58.62	58.58	0.02	426.25
SO	58.27	66.03	62.05	2.93	251.08
SCA	58.27	70.23	60.96	4.08	144.18

* SO-SCA #1: the hybrid of SO and SCA using the first approach, SO-SCA #2: the hybrid of SO and SCA using the second approach.

. For the IEEE 24-bus system, SO-SCA #1 gave a lower standard deviation than SO-SCA #2, as illustrated in

Table 3. Simulation results of SO-SCA #1, SO-SCA #2, SO, and SCA for the 24-bus system.

Optimizer	Best Cost (×106 $)	Worst Cost (×106 $)	Mean Value (×106 $)	Standard Deviation	Execution Time (S)
SO-SCA 1	73.32	87.32	78.89	4.84	348.81
SO-SCA 2	73.32	89.32	79.56	6.18	551.45
SO	77.79	106.98	84.67	11.15	370.12
SCA	77.79	84.47	81.09	2.01	261.59

. Regarding SO and SCA, SCA achieved better performance than SO for the 24-bus system, but SO was superior to SCA in solving the Garver network.

Table 2 and Table 3 explain that SCA converged faster than other algorithms. The time required by SCA to execute one run was about 144.18 s for the Garver network and 261.59 for the IEEE 24-bus system. While SO-SCA #1 ranked second in the convergence speed. It consumed 227.15 s to do one run on the Garver network and 348.81 s on the 24-bus system. SO-SCA #2 was the slowest. It needed 426.25 s to solve the Garver network and 551.45 s to solve the 24-bus system.

Based on all results, it can be concluded that the hybrid algorithm is efficient at solving TEP problems compared to the individual algorithms like SO and SCA. Both hybridization frameworks moved SO and SCA in each iteration to a better search area exploiting the best solution obtained by one of them. The results revealed that the series approach is recommended for solving small-scale systems, while the parallel process is preferred in solving larger systems.

4.3. Testing the Performance of SO-SCA-Based Load Forecasting Technique

The performance of SO-SCA-the based load forecasting technique was compared with PSO-, SCA-, SO, and least error squares (LES) estimation-based methods. The estimated parameters are given in

Table 4. The linear model estimated parameters based on different techniques.

Coefficients	Linear Model
	SO-SCA #1,2	SCA	SO	PSO [21]	LES [21]
α	383.00	377.875	377.841	377.841	363.16
β	1683.00	1713.750	1714.226	1714.226	1874.8
γ	-	-	-	-	-
Average error (%)	3.803074	3.832615	3.8336	3.8336	4.2103

and

Table 5. The quadratic model estimated parameters based on different techniques.

Coefficients	Quadratic Model
	SO-SCA #1,2	SCA	SO	PSO [21]	LES [21]
α	−7.3939	−8.283216	−8.2818	−8.2818	−7.1169
β	494.0909	507.43006	508.4816	508.4816	491.27
γ	1490.3939	1467.2727	1455.493	1455.493	1469.1
Average error (%)	3.323598	3.3375	3.3488	3.3488	3.3873

. The predicted load using all algorithms is depicted in Figure 8.

Table 4 and Table 5 show that SO-SCA#1 and SO-SCA#2 gave the least average error of 3.803074% for the linear model and 3.323598% for the quadratic model. The results also illustrated that SCA was superior to SO, PSO and LES in terms of average error. The linear model was 3.832615% for SCA, 3.8336% for SO, and 3.8336% and 4.2103 for PSO-and LES-based techniques. While for the quadratic model, it was 3.3375%, 3.3488%, 3.3488%, and 3.3873% for SCA-, SO-, PSO-, LES-based techniques.

The results demonstrated that the quadratic model was more suitable for fitting the testing data than the linear model because the minimum average error was 3.803074% for the linear model and 3.323598% for the quadratic model.

4.4. Testing Results of Garver Network

Table 6. TEP’s projects for Garver network.

New Projects	L1 Security	L2 Security
TEP	3-5 (1), 2-6 (3), 4-6 (2), 5-6 (1)	2-3 (1), 1-6 (2), 3-4 (1), 3-6 (2), 4-6 (2), 5-6 (2)
Total ESS (MWh)	198.53	192.79
Total FCLs (p.u)	1.24 + 12.4i	0.59 + 5.9i
Total cost (×106 $)	149.72	175.78

details the results obtained using the proposed planning scheme. The results showed a significant change for projects required to achieve the L1 and L2 security constraints. The maximum number of circuits built in each route did not exceed 2 when the triggering event L2 was applied. In contrast, it did not surpass 3 for the triggering event L1. The total number of circuits installed was 7 for the event L1 and 10 for the event L2. The results emphasized that to achieve L1 security, there was a need for a new circuit at routes 3-5 and 5-6, two circuits at route 4-6, and three circuits at route 2-6. While meeting L2 security required new circuits at routes 2-3, 1-6, 3-4, 3-6, 4-6, and 5-6.

The results showed that the use of ESSs was essential to maintain power continuity and avoid cascading failures. The results indicated that NaS’s capacities required for the events L2 and L1 were 198.35 and 192.79, respectively. The best locations for the ESSs, to achieve L1 security were at nodes 2 and 5. Whilst the locations of ESSs needed to achieve L2 security were at buses 2, 5, and 6.

Figure 9 shows the short circuit current at scenario number24 as an example to verify that it did not exceed the maximum limit (6.5 p.u) and that the short-circuit constraints were achieved in each scenario. It is worth mentioning that the total size of FCLs required ranged between 0.59 + 5.9i p.u and 1.24 + 12.4i p.u. The results also illustrated that the buses’ voltage did not violate their limits (0.9–1.1 p.u). Figure 10 points to the voltage at each bus for scenario number 24 as an example.

The previous results imply that the triggering event L2 was the worst, and the probability of occurring complete blackouts was high. The cost of projects needed for fulfilling L2 security exceeded 17.41 % compared to the cost of L1 projects.

4.5. Testing Results of the IEEE 24-Bus System

Similar to the Garver network, simulations on the 24-bus system confirmed the gravity of the L2 event compared to the L1 event. This is clearly shown in the size of ESSs required and the planning cost. There was an increase of about 108% in the size of ESSs to meet L2 security. The results pointed to the size of NaS batteries required for ensuring L1 security was MWh 29.43 at bus 6. Whilst 60.123 and 1.172 MWh of NaS batteries at buses 6 and 21 were needed to fulfil L2 security.

According to

Table 7. TEP’s projects for the 24-bus system.

New Projects	L1 Security	L2 Security
TEP	1-2 (1), 1-3 (1),2-4 (1), 2-6 (1), 7-2 (1), 7-5 (1), 7-8 (1)	1-2 (1), 1-3 (1), 2-4 (1), 2-6 (1), 7-2 (1), 7-5 (1), 7-8 (1)
Total ESS (MWh)	29.43	61.29
Total FCLs (p.u)	4.1 + 23.89i	5.36 + 28.85i
Total cost (×106 $)	575.36	579.12

, new circuits were installed at the routes 1-2, 1-3, 2-4, 2-6, 7-2, 7-5, and 7-8 to guarantee L1 and L2 securities. It is worth emphasizing that the reliance on the location of new circuits was not enough to achieve the short circuit-current constraints. The use of FCLs was necessary. The short-circuit current at scenario number 24, as an example, is depicted in Figure 11. The results indicated that the planning model was able to maintain the buses’ voltage within their limits. They did not violate the upper and lower bounds (0.9–1.1 p.u), as shown in Figure 12.

Based on all previous findings, it can be concluded that relying on a single strategy to avoid the cascading failures of the electrical networks does not guarantee their security and resilience. Two simple strategies were applied in this work, which significantly affected the size and number of projects required and hence the cost. For the Garver network, the impact was clearly shown in the number and the location of new circuits. While for the 24-bus system, the size of ESSs was affected by the type of cascading failures

The TEP problem is complex, and the hybrid of two meta-heuristic algorithms effectively solved TEP models compared to individual algorithms. The two hybridization frameworks were superior in terms of standard deviation and quality of solutions. Further, the hybrid algorithm was more accurate in forecasting the electrical loads in the long term than some methods introduced in the literature.

5. Conclusions

In this work, a resilience TEP model was suggested to ensure networks’ reliability and supply the future electrical demands. The short-circuit current and the cascading failure constraints were considered to achieve the resilience requirements. NaSs and FCLs’ planning models were incorporated into the proposed TEP model to achieve these targets. Two simple scanning matrices, L1 and L2, were modeled as triggering events for cascading failures. L1 simulated a power system under a single circuit’s outage between two nodes. L2 focused on the outage of all circuits between two nodes as the triggering event for cascading failures. The results showed that simulating a single strategy for cascading failures of power systems does not guarantee their security and resilience. The triggering event L2 was the worst. The probability of occurring complete blackouts was high. The impact was clearly shown in the number and the location of new circuits for the Garver network, while it affected the size of ESSs installed for the IEEE 24-bus system.

The proposed TEP problem was a large-scale non-linear optimization problem, and a hybrid of two meta-heuristic techniques, namely SO and SCA, was adopted. Two hybridization approaches were applied to enhance the exploration and exploitation stages. SO and SCA operated parallel in the first approach to calculate new variables in each iteration. While in the second approach, SO and SCA operated in a series mechanism. The results revealed that the second approach effectively solved small-scale systems, whereas the first approach efficiently solved larger systems.

Further, the results proved that both approaches were superior in terms of standard deviation and quality of solutions compared to individual algorithms. A long-term load forecasting technique based on SO-SCA was proposed. It was shown that the suggested method gave accurate results compared to PSO, SO, SCA, and LES.

This work is only concerned with investigating two simple cascading failure events, but more events that pose a greater threat to the network should be considered to get a secure network. More studies are required to evaluate the impact of including the transient stability constraints on enhancing power systems’ securities and resilience. In addition, the planning model needs to be extended to study the effects of the optimal mix of transmission lines, distributed generators, FCLs, ESSs, reactive source compensators, demand response, and thyristor-controlled series compensators on improving power systems’ securities and avoiding cascading failures. The performance of the hybrid SO-SCA should also be tested on larger systems with an extensive search space.