Multi-objective resilient-constrained generation and transmission expansion planning against natural disasters

In this paper, the resilient-constrained generation and transmission expansion planning (RCGTEP) considering natural disasters such as earthquake and flood is presented. The proposed model contains a two objective problem that the first objective function minimizes the construction and operation costs of resiliency sources (RSs) and, in the second objective function, the minimization of the expected energy not supplied (EENS) due to the outage of the system against the mentioned natural disasters is formulated. Also, the proposed model is constrained to linearized AC optimal power flow (AC-OPF) equations, the planning and operating of RSs constraints such as generation units, hardening transmission lines, parallel and series FACTS devices, resiliency constraints and power system angular stability against earthquakes and floods in the presence of critical and noncritical loads. In order to achieve the best compromise response, the RCGTEP single-objective problem is formulated based on the ε constraint-based Pareto optimization. In addition, this model has the uncertainty of RSs availability against natural disasters, so stochastic programming is used for RCGTEP. Then, the Benders decomposition (BD) method is used to solve the proposed problem and achieve the optimal solution in the shortest possible time. Finally, by implementing the proposed model on an IEEE standard transmission network, the numerical results confirm the capability of this method in improving the economic, operation, angular stability and resiliency indices of the power system simultaneously.


Introduction
One of the most important power system outage causes is natural disasters. Hence, various researches on the destructive effects of natural disasters on the power system have been done in recent years. Their main important purpose is determining the causes of widespread outages and finding methods to strengthen the system against natural disasters [1]. Notably, the characteristic of natural disaster blackouts is different from the characteristic of internal equipment failures in power systems. Therefore, a large area of the power system will out of service which requires a long time for the restoration of this system [2]. But, the blackouts due to internal factors such as equipment failures can cover a smaller area by proper network control and protection. Also, the cause of the error is eliminated immediately [3]. According to these conditions, general techniques and strategies which are used to maintain the stability of the power system against the occurrence of common faults, have many challenges in the face of natural disasters. Therefore, the concept of resiliency in power systems has been considered modeled by the robust method. In [8], a two-objective GTEP is used to minimize the cost of construction and operation of generation units in the first objective function and the network losses in the second one. In [9], the reliability based on N − 1 contingency due to the occurrence of internal faults in transmission equipment is investigated in which to establish adequate reliability the high operating cost is required. The same plan is done in [10], however, it is assumed that the network outage is caused by physical attacks that GTEP is modeled according to the minimum system vulnerability or the maximum system resiliency against these attacks. In [11], the risk-based dynamic GTEP is presented, where it uses post-contingency load-shedding penalty costs in the objective function according to a mixed-integer second-order cone programming (MISOCP) model to penalize high-risk contingencies more dominantly. In [12], PSEP is modeled considering the flexibility of the power system. Note that GTEP modeling in [6][7][8][9][10]12] is based on DC power flow (DC-GTEP) which generally uses the mathematical approaches (MA) such as the simplex method to determine the optimal solution. Noted that in DC-GTEP the changes in voltage, reactive power and network losses are not considered, some network indices such as voltage security can't be evaluated. To solve this issue, some studies such as [13][14] have modeled the GTEP based on AC power flow (AC-GTEP). In [13], the improvement of voltage stability and power losses are considered, and in [14], the voltage stability based on the L-index method has been added to the GTEP. In most studies such as [13][14], evolutionary algorithms (EAs) are used to achieve the optimal solution in AC-GTEP. Another limitation in DC-GTEP is the inability of modeling FACTS devices in the PSEP problem. This issue has been resolved in [15], in which the effects of these devices on economic and operation indices are studied. In [16], AC-GTEP is modeled along with energy storage planning in order to achieve high flexibility and reliability in the transmission network. In [16], a mixed-integer linear programming (MILP) AC-GTEP model is obtained which in the next step is solved by conventional math methods. Generally, in most studies on PSEP such as [6][7][8][9][10][11][12][13][14][15][16], the objective is the improvement of economic and operation indices and in some works, other indices such as flexibility, reliability and security have also been evaluated. Another important index that can be defined for a power system is resiliency against natural disasters like flood and earthquake. Therefore, PSEP can be designed by considering this objective, in this case, the electricity and gas network expansion planning is formulated by taking into account the resiliency index in [17]. In [18] and [19], the resiliency based on graph theory is studied in PSEP and distribution network, respectively. Also, the distribution network resiliency model in the format of the stochastic and hybrid stochastic/robust optimization is respectively presented in [20] and [21], where these papers used the expected energy not supplied (EENS) index to investigate the network resiliency. The reason for this is that shutdown in the network may arise from the outage of network equipment due to natural disasters. Since the improvement of network operation in the condition of natural disasters refers to high network resilience, an index that can be used to assess resiliency in these conditions is the EENS. Because this index refers to the rate of shutdown in the event of N -k contingencies caused by natural disasters. In [22], the problem of transmission line hardening planning as probabilistic power flows considering renewables has been studied. To tackle the issue of incomplete probabilistic data of renewables, a data-driven approach to approximate the renewable uncertainty sets has been developed. Afterwards, it extends the N − 1 security criteria and provides a hardening plan for the worst-case scenarios. A two-stage data-driven stochastic model was formed by taking into account the joint worst-case wind output distribution and transmission line contingencies. Then, it Nomenclature 1) Indices and Sets g, n, j, b, p, s, l, d, m, r, k, h Indices of generation unit, bus, bus, line, parallel FACTS, scenario, load level, load block, feasibility cut, infeasibility cut, polygon side, linearization parts in conventional piecewise linearization method, respectively.  Based on the research background and Table 1, it can be seen that there are the following major research gaps for PSEP: -In most researches such as [6][7][8][9][10][11][12][13][14][15][16], economic and operation indices and in some cases reliability index have been considered. But, some outages in the power system are due to external factors such as natural disasters, so the power system should be developed in such a way to prevent extensive outages in these events. In a few studies such as [10,[17][18][19]22], the resiliency of the transmission network is considered, but they generally consider hardening transmission lines and generation units as RSs while there are different RSs in power systems. -In many studies, economic and operation indices in the PSEP problem have been considered simultaneously. A power system has different economic, operation, reliability and resiliency indices which the improving of one index does not necessarily improve another. For example, cheap generation resources should be used in order to improve the economic index which may not lead to the desired resiliency in the power system. So, to achieve an optimal power system, it is necessary to model indices simultaneously in the PSEP problem. -In available researches for PSEP such as [6][7][8][9][10]12,[17][18]22], the modeling is based on the DC power flow. However, parallel and series FACTS devices can be used as RSs, this model is not applicable for PSEP in this condition. In other words, it is expected that these devices can maintain the voltage/angular stability of the power system with their proper performance in natural disasters. But, these targets are not achievable in the DC power flow model because of the elimination of the reactive power. Therefore, the PSEP based on AC power flow (AC-PSEP) is required. -In most researches such as [13][14], evolutionary algorithms have been used to solve the AC-PSEP problem. Powerful evolutionary algorithms are able to achieve a reliable optimal response by a significant number of iteration to converge which is time-consuming. In addition, they generally have some adjustment parameters that need to be corrected by trial and error when the problem data are changed.
To solve the first and second research gaps, resiliency-constrained generation and transmission expansion planning (RCGTEP) as a twoobjective problem is presented in this paper. The first objective function of the proposed problem formulates the cost of the RSs construction and the cost of generation units operation. In order to model the resiliency of the system, the minimization of the expected energy not supplied (EENS) in the N -k contingency occurring under floods and earthquakes is considered in the second objective function. In this problem, RSs include hardening transmission lines, generation unit and series and parallel FACTS devices. To resolve the third research gap, the proposed RCGTEP is constrained to AC-optimal power flow (AC-OPF), the planning and operating of RSs constraints and the limitation of resiliency and system angular stability against floods and earthquakes. Due to the uncertainty in RSs availability, stochastic modeling is used for RCGTEP. Therefore, a certain number of scenarios for the RSs availability parameters proportional to the Bernoulli probability distribution function (PDF) which depends on the forced outage rate (FOR) of RSs against floods and earthquakes, is generated by the Roulette wheel (RWM) mechanism. Next, the single objective RCGTEP model is derived based on the ε constraint-based Pareto optimization method, and it minimizes the planning and operating costs of RSs. The EENS objective function is modeled in the new problem as a constraint which is limited to a certain constant. The fuzzy decision-making process is also employed to achieve a compromise solution between the mentioned cost and EENS. To compensate the fourth research gap, the Benders decomposition (BD) method is used to obtain the optimal RCGTEP solution in the shortest possible time. This method has two parts: the master problem (MP) and the sub-problem (SP). The MP includes the RSs planning problem, but the RCGTEP operating model is formulated in the SP based on the linearized AC-OPF. Finally, the innovation aspects of the proposed scheme can be summarized as follows: -Modeling the RCGTEP as a multi-objective problem based on the BD method considering linearized AC-OPF. -Simultaneous modeling of economic, operation, stability and resiliency indices in RCGTEP problem. -Analyzing the capability of parallel and series FACTS devices in the issue of power system resiliency by considering linearized AC-OPF in RCGTEP.
The rest of the article is organized as follows: Section 2 presents the multiobjective and single-objective model of RCGTEP as well as modeling the uncertainty of RSs availability. The solving process of the proposed problem based on BD is described in Section 3. Then in Sections 4 and 5, the numerical results and conclusions will be mentioned, respectively.

Multi-objective RCGTEP model
In this section, the two objective RCGTP problem is described. In equations (1) and (2), the objective functions of the proposed RCGTEP problem are expressed, where Eq. (1) indicates the minimization of the total annual cost of the transmission network expansion or the annual cost of the RSs construction and the expected annual operational cost of the generation units [16]. In Eq. (2), the EENS minimization is presented to investigate the network resiliency against natural disasters [23]. The EENS is the sum of the load not supplied in the outage hours (OH) in different areas of the network due to the natural disasters such as floods and earthquakes. The AC-OPF constraints are presented in equations (3)-(9) which represent the active and reactive power balancing at each bus, the active and reactive power flowing through transmission lines and technical network constraints such as transmission lines capacity limit, voltage magnitude and angle limit [24][25][26][27]. It should be mentioned that in case of floods and earthquakes, some network equipment such as transmission lines and generation units will be taken off-line. In other words, the N -k contingency occurs. Therefore, it is expected that the power angle stability of generation units will not be in optimal conditions. So, constraint (9) is added to the proposed RCGTEP to keep the mentioned stability in critical situations by constructing new elements. This limitation can improve network resiliency. Another resiliency constraint is mentioned in (10) which means the amount of load not supplied in a zone (d) is the part of the requested load of that zone which is proportional to the ζ coefficient. Loads are classified into critical and non-critical loads in this paper. Because the outages can damage critical loads severely, therefore, they do not tend to be outage even in natural disasters and ζ is zero for them. But, non-critical loads are not damaged seriously if they go off and ζ can be considered greater than zero for them.
Subject to: ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ( V⩽V n,l,s ⩽V : μ v n,l,s , μ v n,l,s ∀n, l, s − θ⩽θ n,l,s ⩽θ ∀n, l, s RSs constraints are mentioned in equations (11)- (17). The capacity curve of generation units is expressed in Eqs. (11) and (12) which refers to the limitation of their active and reactive power [28]. Also, the allowed reactive power limit of parallel FACTS devices is formulated in Eq. (13) [15]. In addition, it should be noted that in equations (5) and (6), the conductivity and suspension of the transmission line due to the variable construction status of this line (x L ) and the series reactance value of series FACTSs with this line are defined as the variables. These variables can be calculated from Eqs. (14) and (15), respectively. Series FACTS capacitive control limit is modeled in constraint (16). Logical constraint, (17), refers to the fact that a series FACTS can only be connected in series with a transmission line, if the line is constructed in the network. In the proposed RCGTEP, if the value of x L ، x SF ، x PF or x G is equal to 1, it means that the construction of hardening transmission line, series and parallel FACTSs and generation units is cost-effective in terms of resiliency, operation and economic conditions, otherwise, their value is 0. In this paper, the cost of existing elements in the network is zero, therefore their construction status always is 1. In equations (1)- (17), λ and μ are dual variables. The availability of transmission lines, generation unit, series and parallel FACTS in the event of flood and earthquake is defined with the parameters u L , u G ,u SF and u PF . If any of variables is equal to 1, it means that it can be installed on the network. But the series FACTS installation status depends on the availability of series FACTS and transmission line based on the constraint (16). Note that the proposed plan envisages a GTEP that is able to reduce earthquake and flood outages by constructing generation units, transmission lines, and FACTS equipment. This issue is also in line with the goals of network resiliency, that is, at the moment of natural disasters, the network equipment should be used in such a way that blackout is limited. The GTEP also proposes a preventive measure to reduce earthquake and flood damages. In other words, in this plan, it is assumed that in the future there is a possibility of earthquakes and floods in the power system, so the GTEP optimally allocates generation units, transmission lines, and equipment based on the reduction of outages and financial damages caused by natural disasters. Therefore, this plan is limited to resiliency, which by taking preventive decisions can increase the resilience of the system at the time of floods and earthquakes in the future. In addition, in this paper, it is assumed that the outage of equipment can cause a supply-demand imbalance in some parts of the network, so angular instability may occur. To solve this issue, FACTS devices are used to steadily maintain the bus voltage angle in the allowable range in accordance with (9). Also, connection of some equipment in the event of natural disasters can reduce consumer downtime. Therefore, the EENS index is used in this section, mainly because this index is used in the event of natural disasters as a resiliency index [20][21]. Overall, low EENS implies high network resiliency. Also, since the limitation of the network angle stability based on (9) in the proposed plan is always satisfied, this network has the desired stability with the construction of resources, transmission lines and FACTS devices.

Single-objective model
The RCGTEP problem presented in the previous section has a twoobjective function framework. In order to achieve an integrated or single-objective formulation, the Pareto optimization technique is needed. There are several methods to perform this technique, but, the ε-constraint method is used in this paper due to its simplicity [29]. The process of the proposed method for RCGTEP is as follows [29]: The cost objective function in Eq. (1) is considered as the objective function of single-objective RCGTEP as Eq. (18). Also, the second objective function, EENS in Eq. (2), is considered as constraint (19) which limits to ε EENS in the new problem. It is noted that there is no limit for changing the position of the aforementioned objective functions. In other words, the objective function of the new problem can be equal to EENS minimization and the cost function can be defined as the constraint limited to ε Cost . The optimal solution is the same in the two stated problems. Other constraints of single-objective RCGTEP are the constraints (3)-(17) shown as (20).
Subject to: In the above problem, ε EENS varies between the minimum and maximum values of EENS, EENS min and EENS max . The amount of EENS min (EENS max ) is obtained from the problem with the objective function defined in Eq. (2) (Eq. (1)) and constraints (3)- (17). The cost and EENS functions have different directions in terms of increase or decrease so that the decrease of EENS corresponds to the increase in cost. Therefore the minimum (maximum) EENS value of the problem presented in section 2.1 is calculated by considering only the EENS objective function (cost). In addition, due to the consideration of different values for ε EENS , different values are obtained for EENS and cost which by placing in a set, the Pareto front can be evaluated [29]. The next step in Pareto optimization is to achieve the best compromise point between the mentioned objective functions. For this purpose, the fuzzy decision method (FDM) is used in this study in which its details are presented in algorithm 1 [29].

Algorithm 1 Algorithm of the FDM
The best compromise solution of a Pareto front;Pareto optimal solution along with the preferences of the decision-maker; Step 1: Computation of the Fuzzy membership functionCalculation the values of the linear fuzzy membership function (Fi) for each member of the Pareto optimal front: The fuzzy membership function will be equal to 1; else if F i min ≤ F i ≤ F i max The fuzzy membership function will be equal to

Modeling uncertainty of natural disasters
In this paper, it is assumed that there is the possibility of natural disasters such as floods and earthquakes in different areas of the transmission network. These events may cause that the transmission lines, parallel and series FACTS devices and generation units are disconnected from the network. The availability of these elements in natural disasters are displayed by the parameters of u L ، u SF ، u PF and u G , respectively. Due to the uncertainty in the availability of these elements in natural disasters, they will be uncertain parameters. These parameters have two values, 0 and 1. The value of zero for an element means the element inaccessibility (the element disconnection from the network) and the value of one means the element accessibility (the element connectied to the network). The probability of inaccessibility in floods and earthquakes is proportional to the Bernoulli probability distribution function (PDF) [30]. For example, it is expressed as Eq. (21) for a generation unit.
In Eq. (21), π G and χ G represent the outage probability and the forced outage rate (FOR) of the generation unit due to floods and earthquakes. π 0 is the probability, that all the elements are available as shown in Eq. (22).
The stochastic programming is used to model the mentioned uncertain parameters. First, a certain number of scenarios with different values are generated by the Roulette wheel mechanism (RWM) that different values for uncertain parameters, u L ، u SF ، u PF and u G , are considered in each scenario. Then, the probability of each uncertain parameter is calculated from Eq. (21) and the occurrence probability of each scenario is equal to the probability multiplication of u L ، u SF ، u PF and u G . Next, RWM selects scenarios that are more than a certain probability of occurrence, usually 80%.

Solution method
The RCGTEP problem presented in Section 2.2 has a mixed-integer nonlinear programming (MINLP) model. Since the solvers of this model are based on iterative numerical methods such as BONMIN [31], it is expected that the computational time will be high [16]. When the amount of network data is large, using algorithms to solve the problem cannot achieve the optimal situation [23]. It is predicted that the solution methods based on the decomposition technique will be a good solution which this paper, uses the Benders decomposition (BD) method. In this method, the RCGTEP problem is divided into a master problem (MP) and a sub-problem (SP) [32][33]. The RSs planning model including the first three terms of the objective function (1) and constraint (17), series FACTS operation constraint (16), and the equations of the conductivity and suspension of transmission lines, (14) and (15), are considered as MP. Also, transmission network operation limits in the presence of RSs, the last term of the objective function (1) and the constraints (3)-(13) and (19), are modeled in SP. Therefore, the SP problem is as a nonlinear programming (NLP), which will not have the complexity of MINLP problem. Also, by calculating the conductivity and suspension of transmission lines in MP, the degree of complexity of solving relations (5) and (6) in SP will be reduced.
A) Master problem (MP): In MP, the minimization of the RSs construction cost constrained to series FACTS operation and the conductivity and suspension of transmission lines is presented, which its model based on the RCGTEP problem, (18)- (20), is as follows: Subject to: (25) Constraints (14)-(17) The MP objective function is stated in equation (23), which is equal to the constructing cost of RSs based on the constraint (24). The conductivity and suspension of transmission lines, (14) and (15), the series FACTS operation and planning constraints, (16) and (17), are repeated in constraint (25). It is noteworthy that the problem (23)-(25) is known as the initial MP (MP1) [32]. In the following, (26) refers to the Benders feasibility cut [33]. In other words, if SP has an optimal boundary solution, but the convergence point of the BD method is not obtained, the feasibility cut constraint in proportion to the SP results should be applied to MP [32]. In addition, if the SP does not have an optimal boundary solution, then the Benders feasibility cut constraint as (27) in proportion to the SP results should be applied to MP [32]. In these constraints, J SP is equal to the amount of the SP objective function and λ SP and μ SP are dual variables of SP. Noted that the MP model is MINLP due to binary variables of x L , x SF , x PF and x G and non-linear constraints (14) and (15) that calculate the conductance (G L ) and susceptance (B L ) of the transmission line. Since there are no limitations on G L and B L , thus, the MP model as MINLP is not complex. Hence, it is not needed to obtain a simple or linearized model to MP. Finally, it should be said that the MP output variables are the binary variables, x L , x SF , x PF and x G , and continuous variables of G L , B L and X SF which are applied to SP as the parameters [32]. B) Sub-problem (SP): In SP, the problem of resiliency-constrained transmission network operation in the presence of RSs is formulated. Accordingly, the SP model based on the single-objective RCGTEP problem, (18)- (20), is as follows: Subject to: (29) Constraints (3)- (13), (19) The feasibility region of the problem, (28)- (29), depends on the amount of MP output variables which in turn can complicate the solution of SP [33]. Therefore, the dual sub-problem (DSP) model is commonly used to address this issue [32]. But achieving a dual model for the NLP problem has the complexities of a duality gap and equilibrium constraints [34]. Therefore, the best solution is to use the linear approximation model for the NLP problem, (28)- (29), which the dual model can be obtained subsequently.
In SP, constraints (5), (6) and (7) are nonlinear. Since the voltage magnitude, V, in the transmission network must vary between 0.95 and 1.05p.u, so according to [35], the terms V 2 and V n V j can be approximated to the linear equation 2 V − 1 and V n + V j − 1, respectively, while the maximum approximation error is 0.0025p.u. In addition, the functions Ψ 1 and Ψ 2 in constraints (5) and (6) can be expressed as Δδ h by using the conventional piecewise linearization method [25], respectively. In these equations, Δδ represents the deviation of the power angle, and β is the line slope. Also δ is the minimum value of the power angle which is equal to −2θ. Note that the power angle (δ) is equal to θ nθ j and it can be expressed as δ + ∑ Furthermore, Ψ 1 V n V j or Ψ 2 V n V j can be defined as V n V j Ψ 1 Δδ h by accepting a maximum approximate error of 0.005. Therefore, based on the mentioned statements, constraints (5) and (6) can be approximated as linear equations (30) and (31), respectively:

j,l,s,h
: λ pl n,j,b,l,s ∀n, j, b, l, s In this case, the angle stability constraint, (9), can also be rewritten as (32): Constraint (7) represents a circular plane in P and Q coordinates, ⩽S. To achieve a linear equation, this plane can be approximated to a polygon [26]. If the number of sides is high, the approximation error will be very low. It should be said that the equation of each side of the polygon is P.cos(k.Δα) + Q.sin(k.Δα) = Sas a linear equation. Δα indicates the angle deviation equal to 360 / N e . N e is equal to the number of sides of a polygon and k is the number of a regular polygon side that is selected from the set Ξ K = {1,2,...,N e }. Also, the plane corresponding to the side k is as P.cos(k.Δα) + Q.sin(k.Δα)⩽S. Finally, a polygon plane will be obtained by sharing the corresponding pages with all sides. Hence the linear equation of constraint (7) will be as follows: Subject to: (35) Constraints (3), (4), (7), (9)-(13), (19), (30)- (33) Finally, the DSP model based on the dual theory presented in [32] will be written as follows: Equation (36) is equal to the DSP objective function, which is equal to J SP . Constraints (37)-(44) are equal to the dual equations corresponding to the SP problem variables. In each equation, the corresponding SP variable is presented. The range of changes of the dual variables λ and μ is also expressed in equation (45). Finally, there are three ways to solve DSP, which are as follows [32]: 1) DSP has an optimal boundary solution: In this case, the feasibility cut constraint is added to MP as the Eq. (46).
2) The DSP does not have an optimal boundary solution: In this case, a constraint similar to (47) is considered for the dual variables λ and μ.
Then the feasibility cut constraint is added to MP as equation (48) according to the results of (47).
3) DSP does not have an optimal solution: In this case, problem (1)- (17) does not have a feasible solution. In equations (46) and (48), λ SP and μ SP are the optimal value of the variables λ and μ. In addition, the convergence condition for the BD method is available if |Z U − Z L |⩽ε is met [32]. In this equation, ε indicates the convergence tolerance or computational error allowed in the BD method. Z U and Z L also represent the upper and lower bounds of the BD method, where Z L is determined from MP, and Z U is equal to the right-side of (46). Finally, the RCGTEP solution flowchart based on the BD method is as shown in Fig. 1.

Numerical results
This section presents the numerical results of the proposed RCGTEP on the modified IEEE 6-bus and the IEEE 118-bus transmission networks. The problem is coded in GAMS optimization software, which uses the BONMIN solver to solve the MP, and the SP solution will be based on the CPLEX algorithm [31]. To achieve the guaranteed optimal solution, the number of linear pieces for functions Ψ 1 and Ψ 2 is assumed to be 8, and the circular plane is approximated to a 45 side polygon. It is also assumed that the convergence conditions for the BD method can be achieved by accepting tolerance of $ 3 (ε = 3).

Modified IEEE 6-bus network
A) Problem data: A single-line diagram of the modified IEEE 6-bus network is shown in Fig. 2 [16]. The characteristics of existing and candidate lines including their resistance, reactance, construction cost and capacity are presented in [16]. In this paper, it is assumed that the resiliency of existing transmission lines against floods and earthquakes is low and FOR for these existing lines is 10%. While, the resiliency of candidate transmission lines against floods and earthquakes is high and FOR for these existing lines is 0.1%. It is also assumed that a candidate transmission line can have an series FACTS whose capacitive reactance varies between zero and 20% of the candidate line reactance. Its construction cost is equal to 30% of the transmission line cost and FOR is equal to 0.1%. In addition, this network has 4 existing non-renewable generation units and 13 candidate non-renewable generation units, whose specifications such as the construction location, fuel price, construction cost and capacity are presented in [16]. The FOR rate for the existing and candidate generation units is 3% and 40%, respectively. Their minimum generation capacity is zero. The maximum active and reactive generation capacity and the minimum reactive power of each unit are 80%, 50% and −50% of the unit capacity, respectively. In this  network, three parallel FACTSs including static synchronous compensator (STATCOM) and static var compensator (SVC) can be installed in bus 3-5 that their information such as capacity, type and installation cost are stated in [15]. In this paper, the series and parallel FACTS resiliency against floods and earthquakes is high, so FOR is equal to 1%. Note that in Eqs. (20) and (21), the term FOR for any equipment means its forced outage rate against natural disasters such as floods or earthquakes. In other words, in these equations, the FOR value for an element against different natural disasters can be different. However, in this study, the FOR term is assumed to be the same for floods and. However, there is no limit to the choice of different FOR values for different natural disasters. The upper and lower bound for the voltage magnitude is 1.05 and 0.95p.u, respectively [23]. In order to maintain stability angle in the events of floods and earthquakes, the voltage angle must be within the allowable range of [90,-90]. In the mentioned network, buses 3, 4 and 5 have loads with the capacity of 40%, 30% and 30% of the total network load, respectively [30]. Also, the network load has five levels during a year which the load coefficients of levels are 0.5, 0.65, 0.8, 0.9 and 1, respectively. The event time interval of load levels is 1510, 2800, 2720, 1120 and 610 h, respectively [30] and the power factor is often considered to be 0.9. It is also assumed that the load on bus 4 is the critical load, but the load of buses 3 and 5 are the non-critical one. Therefore, the value of ζ for critical and non-critical loads is considered equal to zero and one, respectively. Planning is for 5 years and floods or earthquakes are likely to occur once in these 5 years. Finally, it is assumed that there is a possibility of the flood (earthquake) in buses 1 to 3 (4-6) and the lines between them which in the worst-case scenario, the outage hours (OH) is equal to 720 h. In this section, RWM creates 30 scenarios for modeling the uncertainty of network equipment availability against floods and earthquakes based on Bernoulli PDF which is applied to the proposed RCGTEP problem. B) Evaluation of best compromise solution: The Pareto front results for the proposed RCGTEP problem, (23)- (27) and (36)-(45), for the network peak load (NPL) of 25, 50 and 75 MW are shown in Fig. 3. According to this figure, the minimum value of EENS, i.e. zero, is proportional to the maximum value of the cost function, which is obtained around $ 30, 39 and 60 million when NPL is 25, 50 and 75 MW, respectively. However, the lower values of cost correspond to the increase of EENS, so that the lowest value of cost is at the point on the Pareto front where EENS has the highest value. According to the data presented in the previous section, RSs are located in different locations of the network which their FOR rate is much higher than the network equipment, therefore, a highresilient transmission network will be costly. Therefore, the reduction of EENS in the event of floods and earthquakes requires high investment in network expansion.
The results confirm the second research gap presented in Section 1. Finally, the results of the best optimal solution based on the fuzzy decision method, algorithm 1, for the proposed RCGTEP model are summarized in Table 2 3%. This indicates the ability of the proposed model to achieve the best compromise response between the network expansion costs and EENS which is in accordance with the first and second innovation in Section 1.

C) Evaluation of economic and technical indices:
The results of the network expansion planning including the investment and operating costs and, and computational time for different NLP values and different RCGTEP problem models according to the best compromise solution are shown in Table 3. It is noted that the MINLP model in this section which includes equations (18)- (20) represents the real model of the proposed problem. Also, the AC-RCGTEP model based on the BD method has the same equations and technique used in Section 3. The DC-RCGTEP model based on the BD method is also based on the equations and techniques of Section 3, with the difference that the reactive power variable is zero, the buses' voltage is 1p.u, the resistance and conductivity of the transmission lines is zero, and the phrase sin ( θ n,l,s − θ j,l,s ) is approximated to ( θ n,l,s − θ j,l,s ) . As shown in Table 3, the results obtained for the MINLP and AC-RCGTEP models based on the BD method are close to each other. The constructed elements based on the calculations of the two models are the same, and the amount of computational error for the operation cost of the AC-RCGTEP model based on the BD method has a higher value (0. Although the computational time in DC-RCGTEP based on the BD method has the lowest value compared to other models, due to high computational error is not acceptable for the RCGTEP problem. The AC-RCGTEP model based on BD method, due to its low error and computational time, is the best model for the RCGTEP problem, which indicates the capabilities of the first and third innovations in Section 1. Finally, according to Table 3, it can be stated that the increase in NPL increases the number of constructed RSs which increases the construction cost and the expected operating cost. Also, due to the reduction of the problem feasibility space for the increase in NPL, the computational time of the problem also increases. It should be stated that problem (18)-(20) is a NP-hard problem, because Eqs. (5) and (6) are strictly nonlinear due to the variables G L and B L , and also the problem has binary variables besides continuous variables. Therefore, the solution time will be very high. But the solution method used in this paper is based on BD. The constraints that complicate problem solving, (14)- (17) and the binary variables, are placed in the MP section, which eventually form a simple problem as (23)- (27). In addition, the remainder of problem (18)- (20) used in SP or (36)-(45) is a linearized optimal power flow (OPF) problem, which has a high calculation time. Finally, the model used in MP and SP does not have the same complexity as (18)- (20), so this has caused the calculation time to be very low. This is also confirmed in [25,[32][33].
Results of technical indices for one year (assumed to be proportional to the year in which floods and earthquakes occurred) include energy Table 2 The best compromise solution results.  Table 4. It is noteworthy that MVD can be calculated using the following formula: max (|1 -V|), and MVA is equal to max (|θ|). As can be seen, the results of the operation indices, namely EL and MVD, in the AC-RCGTEP model based on the BD method are very close to the results of the MINLP model which the highest computational error for EL and MVD is 32.9 MWh (13888.7-13921.6) and 0.005p.u (0.0235-0.024), respectively. However, in DC-RCGTEP based on the BD method due to the elimination of line resistance and voltage drop of the buses in this model, the values of EL and MVD are equal to zero. The comparison of this case with the results of the MINLP model shows that the computational error in the DC-RCGTEP model based on the BD method is very high. Therefore, the results shown in Table 4 also confirm the appropriate status of the AC-RCGTEP model based on BD method. In addition, the increase in NPL increases energy losses, MVD, MVA, and EENS. Note that at the highest NPL value, 75 MW, the system is statically stable, because the absolute value of the buses voltage angle is always less than 90 degrees. Also, the bus voltages are in their allowable range, i.e. 0.95 to 1.05p.u. EENS is also around 77% away from its maximum value based on the results of Section 4.1.B. Moreover, the ratio of losses per hour to NPL is less than 3%. Therefore, this issue indicates the ability of the proposed model which is in line with the first and second innovations in Section 1. D) Determining the outage cost in the event of floods and earthquakes: One of the important indices to determine the system resiliency can be the outage cost (OC)     considered as a parameter. The results are shown in Fig. 4. Based on this figure, it can be seen that increasing VOLL decreases EENS and increases cost. Although, increasing it to a certain value such as 3000 $/MWh (3 k $/MWh) for NPL 25 MW will increase the OC, but increasing the VOLL more than this will reduce the outage cost. Finally, the suitable penalty price for achieving high resiliency, i.e. EENS and OC is zero, in the system of Fig. 2 for NPL 25 MW according to Fig. 4 is equal to 10 k $/MWh, These values are 12 and 16 k$/MWh for NPL 50 MW and NPL 75 MW, respectively. In this case, the cost of power system expansion compared to the results of Section 4.2.B is equal to its maximum value, because EENS in this case is zero.

IEEE 118-bus network
The IEEE 118-bus network has 54 generation units, 186 transmission lines and 91 loads buses. The information of this network is presented in [36]. In this paper, it is assumed that in the location of the existing generation unit, another generation unit with the same specifications can be constructed with a construction cost of 40 $/kVA/year. This assumption also holds for transmission lines, with the difference that their construction cost is 5 $/kVA/year. Like the previous section system information, FOR is considered 30% and 40% for existing transmission lines and generation units, respectively, but for candidate lines and units are 0.1% And 3%, respectively. Each candidate line can have a series FACTS whose specifications are the same as the data presented in the previous section. It is also assumed that a parallel FACTS of SVC type with the capacity of 20 kVAr and the construction cost (FOR) of 4 $/kVAr (1%) can be installed in each bus. The peak load characteristics are presented in [36]. The amount of power factor, load factor at different load levels and other cases are similar to the information in Section 4.1. Planning is for 5 years and floods or earthquakes are likely to occur once in 5 years. In this section, it is assumed that buses 1-59 and the lines between them are in the flood zone, and the buses 18-160 and the lines between them are in the earthquake zones. The OH is equal to 720 h, and RWM creates 60 scenarios in this section.
The results of this section for the following case studies are summarized in Tables 5 and 6: -Case I: RCGTEP problem including generation unit as RS -Case II: RCGTEP problem including generation unit and transmission lines as RS -Case III: RCGTEP problem including generation unit, transmission lines and series FACTS as RS -Case IV: RCGTEP problem including generation unit and parallel FACTS as RS -Case V: RCGTEP problem including generation unit, transmission lines, series and parallel FACTS as RS -Case VI: Case V with high resiliency mode, i.e. EENS and OC equal to zero Based on the results of Table 5, the indices of resiliency, operation, static stability and economic can be improved by using different types of RSs, so that EL, MVD, MVA, EENS, cost are reduced. However each type of RSs has the ability to improve these indices, but if different types of RSs are used in a distributed way over a large area of the network, better results can be obtained. Note that due to the increase of candidate RSs in case V, its computational time will be longer than the other cases. In addition, increasing the NPL increases the values of the various indices but in the allowed ranges. Finally, according to Table 6, it can be seen in case V that the VOLL amount at the best compromise solution for NPL 12,000 and 14000 MW is 7.24 and 8.58 k$/MWh, respectively, followed by the outage cost of 1.18 and 2.07 G$, respectively. But VOLL amount in the highly resilient system, i.e. EENS and OC are zero, for NPL 12,000 and 14000 MW are obtained 17.72 and 20.11 k$/MWh, respectively. In this case, the planning cost will increase by 17% ((16.24-5.53)/16.24) and 61.3% ((18.19-7.04)/18.19) in case VI compared to case V, respectively. Therefore, to achieve a system with high resiliency, high planning cost is required.

Conclusion
This paper addresses the RCGTEP problem in the transmission network in the event of floods and earthquakes. In the proposed model, a two-objective optimization problem is expressed that in the first objective function, the total cost of RSs construction and the expected operation cost of generation units are minimized, and in the second objective function, EENS is minimized. This problem is constrained to the AC-OPF equations, resiliency limits, static stability in the presence of critical and non-critical loads, and operation and planning equations of RSs such as generation unit, transmission line, series and parallel FACTS. Stochastic programming is used to model the uncertainty of network equipment availability against floods and earthquakes in this paper. Then, the Pareto optimization technique based on ε-constraint method is used to achieve the single-objective problem. The problem-solving process is based on the BD method, in which a linearized model is used for the operation problem. Finally, based on numerical results, it is observed that the AC-RCGTEP model based on the BD method has a good ability to achieve the optimal solution with the least computational error in the shortest possible computational time. The proposed problem is also able to achieve the best compromise solution between the planning cost and EENS, so that they are approximately 60% away from their maximum value. With the increase in outage penalty price, VOLL, the system is able to achieve maximum resiliency, i.e. EENS is zero, with respect to the results of the best compromise solution at the expense of 8000 $/MWh, where it comes with a high planning cost. In addition, the proposed RCGTEP model with optimal RSs planning and scheduling is able to improve resiliency, static stability, operation and economic indices simultaneously.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.