Critical Propagation Path Identification for Cascading Overload Failures With Multi-Stage MILP

Cascading overload failures of transmission lines are one of the dominant factors that induce catastrophic blackouts. The reliability risk of cascading overload failures is significantly underestimated in traditional adequacy assessment studies. This paper proposes an evaluation model for the assessment of <inline-formula> <tex-math notation="LaTeX">$N$ </tex-math></inline-formula>-<inline-formula> <tex-math notation="LaTeX">$k$ </tex-math></inline-formula> contingencies as cascading overloading failures. The interactions between outaged lines in the propagation path are revealed with a multi-stage optimization problem. The objective function is to maximize the probability of the propagation path of cascading overload failures. Several linearization techniques are proposed to convert the nonlinear evaluation model to mixed-integer linear programming. A rolling search algorithm is adopted to improve the capability of the proposed model to evaluate high-order contingencies. The case studies on the IEEE 118-bus systems show that the proposed model enables the identification of critical propagation paths with probabilities that are a few orders of magnitude larger than those with the traditional method.


I. INTRODUCTION
Power supply is a fundamental resource in modern society. The blackouts in past decades have shown that a power outage can lead to catastrophic consequences in human society. Many reports and studies of past blackouts reveal that the outage of multiple transmission lines is one of the dominant factors in the propagation process of blackouts [1], [2]. The evaluation of N -k contingency is widely required in the reliability criterion of power system operation and planning. For example, the NERC requires that the transmission system planner should investigate the interaction potential of cascading outages due to overloading or local protection [3], [4].
As reliability describes the ability of power supply from the perspectives of adequacy and security, the risk of N -k contingency has been studied with different evaluation frameworks. For adequacy evaluation, many advanced sampling methods are proposed to improve the search efficiency of high-risk events. [5] adopted the cross-entropy-based importance sampling to increase the sampling frequency of critical events that cause load shedding. [6] proposed a pseudoanalytical mix sampling strategy to avoid the local-scatting effect of samples. [7] used the intelligent searching algorithm to identify the critical N -k contingencies in the system state space. Recently, the adequacy risk of multiple line outages is also widely studied from the perspective of power system resilience under extreme events [8].
The adequacy evaluation is often studied with the assumption that the multiple lines fail independently. As a result, the risk of adequacy evaluation is often dominated by contingencies in the lower order and may underestimate the risk of high order contingency [9], [10]. However, past blackouts have proved that N -k contingencies do occur and consequently lead to catastrophic impact. Thus, it is vital to develop the evaluation framework of N -k contingency with the consideration of security risk.
The security assessment of contingencies cares more about the post-contingency violation of power system operation limits. For transient stability assessment, the security performance index is widely used to quantify the severity of contingency [11]. However, the stability analysis is mostly performed for the N -1 contingency as it deals with the shorttime scale mechanism of cascading failure [12]. Another critical issue in the security assessment is the violation of capacity limits for transmission lines [13], [14]. [15] proposed a fast N -2 contingency screening algorithm to analyze the overload risk with line outage distribution factors. [16] ranked the N -1 contingency with the remedial cost required to correct the violations of overload post-contingency limits.
As the cascading overload failure of transmission lines is one of the dominant factors of blackouts, in most cascading failure models, for example, the well-known OPA model [13] and its variants [17], the propagation of line outages is analyzed by simulating the redistribution of power flow following the outage of critical transmission lines. When the high-order contingencies are considered in the security assessment of cascading overload failure, massive power flow calculations are required. To reduce the complexity of the identification of critical propagation paths, many advanced searching techniques have been developed. [18] constructed a random vector functional-link neural network to estimate the cascading failure risk with simple matrix calculations. Similarly, [19] developed a deep convolutional neural network to perform the AC optimal power flow and screened the cascading failure chains with the depth-first search algorithm. [20] proposed a fast power flow calculation with line outage distribution factor. [21] used multi-parametric linear programming to accelerate the repeated DC optimal power flow calculations. Reference [22] and [23] introduce the random chemistry algorithm to search for the critical initial N -k contingencies in cascading overload failures. Reference [24] designed the reinforcement learning framework to identify the cascading overload failure under multi-stage terrorist attacks. Recently, the complex interaction network theory is also investigated to extract the propagation pattern of cascading failure from simulation results and utility data [25], [26], [27]. These methods greatly release the computational burden. However, the sampling methods or machine learning algorithms may have a robustness limit for identification of the most severe propagation path, which is specially cared by the power system operators.
Recently, mixed-integer linear programming (MILP) is widely applied to the problems of power system reliability assessment due to its robustness. Reference [28] constructed the MILP model to evaluate the adequacy shortage risk of power system under the terrorist-attack induced N -k contingencies. [29] developed an MILP model that can deal with the N -k contingency under the AC power flow constraints. Reference [30] proposed an MILP model in which the objective function evaluates the probability of N -k contingency so that the load shedding risk can be assessed. Reference [31] screened the N -k contingency with high incremental risk so that the masking effect of the risk index can be avoided. Reference [32] proposed a space pruning algorithm with root events to significantly accelerate the contingency evaluation with MILP. However, these studies focus on the adequacy risk of N -k contingencies. The evaluation framework of adequacy risk assumes that the outages of transmission lines occur independently. Independent line outages are usually modeled with a fixed probability of terrorist attacks or random outages. The probability of N -k contingency is calculated as the product of that of independent line outage events in these studies. On the other hand, the cascading overload failure considers a different line outage model, in which the outage probability of a transmission line increases with the power flowing through it. As line outages can lead to power flow redistribution and may cause more outages, there are interactions between line outages for the probability evaluation of N-k contingency as cascading overload failure [10]. To reveal the interdependence of outaged lines in the N -k contingency, [10] developed an MILP model to evaluate the possible propagation path of cascading overload failure. This method does not assess the propagation probability directly due to the nonlinear complexity limitation. Instead, the power flow increment is selected in [10] as an alternative so that the evaluation model can be constructed as MILP. Similar assumptions can be found in the related contingency screening models and risk mitigation models [33], [34], [35]. However, as the probability is not explicitly modeled in [10], the result may miss the actual propagation path of cascading failure with the highest probability. This is more likely to happen in large power systems as the number of possible propagation paths for an N -k contingency is k!, which is enormous when k is large. Table 1 shows a summary of the previous studies on N -k contingency evaluation. These studies are compared from the perspectives of the modeling of failure patterns, contingency analysis, evaluation objectives, and solution methodology.
To address the above problem, this paper develops a novel evaluation model of N -k contingency, in which the interdependence of failure probability of transmission lines is explicitly modeled. Several linearization techniques are adopted to convert the proposed model to a multi-stage MILP, which can be efficiently solved with commercial software such as Gurobi and Cplex. The proposed method provides a new choice for the system operators to analyze the propagation path with the highest risk, which can be complementary to the existing search models. The main contribution of this paper is to improve the evaluation model of N -k contingency in [10], [33], [34], and [35] by avoiding the approximation of the increment of the outage probability with that of power flow. The detailed contributions include: 1) An evaluation model is proposed to assess the propagation paths of N -k contingency. The outage probability of overloaded transmission lines at different stages is directly modeled. The proposed model can identify the propagation path with the largest probability, instead of the largest power flow increment, as that with the highest risk. Thus, it can enhance the capability to reveal the interdependency of outaged lines in the N -k contingency more explicitly.
2) To address nonlinearity from the evaluation of outage propagation probability, several linearization techniques, including the logarithmic transformation and piecewise linearization, are adopted to convert the proposed evaluation model to a multi-stage MILP. The propagation path of cascading overload failures with the maximum probability can be effectively identified.
3) To accelerate the identification of propagation paths, a rolling search algorithm is proposed in which the critical lines are recursively selected with the proposed evaluation model.
The rest of this paper is organized as follows. Section II proposes the probability evaluation model for propagation VOLUME 10, 2022 paths of N -k contingency. Section III provides the solution methodology of the evaluation model with multi-stage MILP. Section IV tests the effectiveness of the proposed model on the IEEE 118-bus system. Section V concludes the paper.

II. PROBABILITY EVALUATION MODEL FOR N-K CONTINGENCY A. OuTAGE PROBABILITY MODEL FOR OVERLOADED TRANSMISSION LINES
In most reliability evaluation studies, the unavailability of transmission lines captures the random outage probability due to equipment aging or severe weather. Generally, as the random outage happens with a low probability, multiple line outage events rarely occur. However, the cascading tripping events of transmission lines are observed in many blackouts. For example, more than 400 lines tripped in the US-Canada blackout on August 14, 2003 [36]. Many blackouts indicate that the outage probability of transmission lines may increase with the power flowing on them. Studies show that the overloading operation of a transmission line could increase its sag and thus increase the probability of tree-contact fault [14]. Moreover, the overloaded line could be tripped by the temperature relay or overcurrent relay [37].
In this paper, the outage probability of transmission line l in cascading failures is modeled with (1) as [10].
Equation (1) indicates that the outage probability of transmission line l is maintained at a relatively low level when it is not overloaded. The outage probability grows until the power flowing on line l reaches the specified threshold. The value of capacity coefficient R ranges from 1.2 to 1.8 [10], [35], [37]. In this paper, R is set as 1.4.

B. EVALUATION MODEL FOR N-K CONTINGENCY
As shown in (1), the outage probability of line l depends on its power flow. Considering the N -k contingency S, the most critical propagation path of outage events in S can be assessed with (2).
During the propagation of cascading failure, the power flowing on the outaged lines is redistributed between the other lines in the grid. In this context, the evaluation of contingency S should consider the power transfer process that causes blackouts with the maximum probability. Thus, a multi-stage evaluation model for N -k contingency is constructed in this section.
Let the first stage of the propagation path denote the precontingency stage. For contingency S, the propagation path can include k stages at most before all the lines are outaged. Suppose line l fails at stage s, namely, line l becomes unavailable since stage s + 1, its failure probability is calculated with the power flow at stage s. A generic objective function is constructed as (3) to assess the multi-stage propagation probability of contingency S.
The binary indicator v s l is used to select the stage at which the outage probability of line l is calculated.
In addition to the objective function (3), the proposed evaluation model includes the power flow constraints and propagation path constraints, which are described as follows.

1) POWER FLOW CONSTRAINTS
The DC power flow is adopted in the proposed evaluation model for the post-contingency analysis. Constraint (4) indicates that the power balance should be satisfied at each stage with the initial power generation and load demand. Constraint (5) sets the power flowing on line l at stage s according to its availability. If u s l is 0, namely, line l is outaged at stage s, F s l is 0. Otherwise, F s l is calculated with the DC power flow.

2) PROPAGATION PATH CONSTRAINTS
For each stage, the outage probability of line l can be assessed with constraints (1) and (5). In addition, the following constraints (6)-(10) need to be considered to describe the propagation path of contingency S.
Constraint (6) sets all the lines in S available at the precontingency stage. Constraint (7) indicates that once line l is out of service at a stage, it cannot be repaired in the following stages. Constraint (8) calculates the status switch of line l at each stage. Constraint (9) indicates that all transmission lines in contingency S should be outaged during the k stages. Constraint (10) indicates that each propagation stage includes one line outage event.
Similar to the assumption in [10], [33], [34], and [35], it should be noted that the N -k contingencies that cause the separation of transmission network are out of the discussion in this paper. The reason is that constraint (1) requires the calculations of DC power flow to assess the outage probability of overloaded transmission lines. However, if the network is separated, the power balance cannot be guaranteed in each sub-grid. The rebalancing of power demand requires generation re-dispatching or even load shedding, which depends on the operation strategy of power systems. Different strategies can lead to different distributions of power flow and thus change the outage probability of transmission lines. In addition, the instability of power system may be caused by the network separation. This may further change the propagation path of contingency from that dominated by the overloading issue [37], [38]. Thus, this paper only considers the N -k contingencies that do not disconnect the power system.

III. SOLUTION METHODOLOGY WITH MULTI-STAGE MILP
The evaluation model proposed in Section II provides an explicit method to assess the severity of N -k contingency as cascading overload failures. However, it is noted that the objective function and constraints of the proposed evaluation model are highly nonlinear, which makes it difficult to solve. Thus, the section proposes several techniques to convert the proposed evaluation model to a multistage MILP for recursive identification of the propagation path of cascading overload failures with the maximum probability level.

A. LINEARIZATION CONVERSION OF THE PROPOSED EVALUATION MODEL
The probability evaluation in the objective function (3) is highly nonlinear. To avoid nonlinear operations, a simplified function that calculates the summation of the increment of power flow in lines before its outage is used in [10] and [33]. Although the line with more increment of power flow is more likely to be outaged, this approximation is not exact enough as equation (1) is a piecewise function. Therefore, a new reformulation of the objective function (3) is proposed in this paper.
Firstly, the logarithmic transformation of (3) is performed to deal with the product of outage probabilities.  (1) is converted into the linear constraints (13)- (20). Here, R s is a scaling coefficient, which is larger than R. In this paper, R s is set as 5.
t s i l, ≤ y The absolute value in (14) can be linearized as constraints (21) and (22).
The product of variables in (5) and (11) can be linearized with the big-M method, which is shown as constraints (23)-(25).
Finally, for the nonlinear logarithmic function in (12), another piecewise linearization function is used. Note that the logarithmic function is convex, its linearization is different from equation (1), and is performed with constraint (26) which requires no integer variable.
where t p is the breakpoint ranging from 0 to 1, and (t p ) is the set of breakpoints.
In conclusion, the linear objective function (27) can be adopted as a substitute for function (3) of outaged lines. This limits the applicability of the proposed model to large contingencies. Thus, a rolling search algorithm is proposed in this section to accelerate the evaluation process.
The idea of the rolling search algorithm is similar to that in the rolling horizon optimization [39], [40]. The rolling horizon is applied to solve multistage or multiperiod optimization problems, for example, the aircraft sequencing problem [39] and the power energy management [40]. The original multistage optimization problem is decomposed into several sub-problems. Each sub-problem covers a prediction horizon T p which is shorter than the total period T . The rolling horizon algorithm starts by solving the first subproblem for the decision at stage 1, and stops until the length of the decision window grows up to T with a step of T p . When the rolling search algorithm is adopted to identify the propagation path of cascading overload failures, the step size is firstly defined as k s . For an N-k contingency S, the number of iterations required for the identification of propagation path is n, for which we have nk s ≤ k < (n + 1) k s .
In iteration 1, the line outages at the first k s −1 stages are identified. And in iteration i (1< i < n), the line outages from stage (i-1)k s to stage ik s -1 are identified. Finally, iteration n identifies the line outages from stage (n-1)k s to stage k. The proposed rolling search algorithm is described as follows: Step 1) Set the step size k s , and initialize the loop flag i as 1.
Step 2) Set the length of propagation path for the N -k contingency as the minimum of ik s and k. When ik s is smaller than k, the evaluation model (M1) is modified as (M2).
Constraints (9) and (30) indicate that the length of propagation path in iteration i is ik s , shorter than k. It sets that the number of line outages at stage ik s is k−ik s +1, and the other stages consist of one line outage. In addition to the line outages at stages 1 to ik s −1, the objective function maximizes the outage probabilities of line outages at stage ik s , for which the outage sequence will be further optimized in the following iterations. This is similar to the decisionmaking with consideration of the further prediction horizons in rolling horizon approaches [39], [40].
Step 3) Solve the evaluation model with the solution as [v * (i) ].
Step 4) Update the loop flag i as i + 1. Step 5) Reconstruct the evaluation model (M2) with the updated step size ik s . In addition, set the status of outaged lines at the first (i-1)k s -1 stages according to the solution v * (i−1) from the last iteration. This is performed by adding the constraint (29) in the evaluation mode (M2). Step 6) The steps 3-5 loop until the maximum length of propagation path reaches k. Fig. 1 gives an illustrative example of the proposed rolling search algorithm with k s initialized as 3.

IV. CASE STUDY
In this section, the proposed evaluation model is tested on the IEEE 118-bus system. The original system data is modified as [10]. Taking the N -10 contingency and N -20 contingency considered in [10] as an example, the proposed evaluation model is implemented on a personal computer with an i9-10850k CPU and 48G RAM. The detailed analysis is as follows.
In the first iteration, the propagation path identified by the proposed model is The computation time for the first and second iterations is 107.95 seconds and 45.17 seconds, respectively. Table 2 lists the power flow results for the identified propagation path.
The outage probability of lines in this N -10 contingency is assessed by substituting the power flow results in Table 2 into equation (1). Table 3 lists the outage probabilities of lines in the identified propagation path. Table 3 also compares the evaluation results of the proposed method with the original outage probabilities and those in [10]. The original   [10] for the N-10 contingency.
outage probabilities are assessed with the pre-contingency power flow. The fifth column of Table 3 shows the ratio of outage probabilities from the proposed method and those in [10].
As shown in Table 3, compared with the propagation path in [10], four lines, including 25-27, 38-37, 64-65, and 22-23, get a larger outage probability in the proposed method. In particular, the outage probability of line 38-37 increases from 0.01 to 1.00. In [10], as the objective function is to identify the propagation path with the largest power increment, line 38-37 is identified to be outaged at the first stage so that the power of 257.63 MW flowing on it is redistributed among other lines. However, the propagation path with the largest power increment does not always imply that with the highest probability. It can be found from Table 2 that the power flowing on line 38-37 increases from 73.61% to 158.65% following the outage of lines 8-5 and 30-17. Thus, the proposed model performs a more explicit evaluation for multiple line outages than [10] by avoiding the approximation of the outage probability of transmission lines with its power flow increment.
By multiplying the probabilities of ten line outages together, the original probability of this N -10 contingency without cascading failures is 2.91 × 10 −17 . After the consideration of the propagation of cascading failures, the probability of this N -10 contingency is assessed as 2.17 × 10 −7 by the proposed method and 3.60 × 10 −12 in [10]. The probability of this N -10 contingency increases in both propagation paths by revealing the interdependency between line outages. In addition, although the computation time for the MILP model in [10] is 117.04 seconds and shorter than that for the proposed method, the propagation path identified by the proposed method reveals a probability of 6.03 × 10 4 times larger than that by the MILP model in [10]. The main reason that the proposed method takes more time is the incremental number of integer variables for the accurate modeling of the probability of cascading overload failure. The results indicate that the improved probability model enables the identification of closer interactions between line outages in the propagation of cascading failures.

B. N-20 CONTINGENCY
The proposed model is also adopted to evaluate the N -20 contingency in [10] Table 4 compares the outage probabilities of this N -20 contingency.
It can be found from Table 4 that eight lines in this N -20 contingency get larger outage probabilities in the propagation path identified by the proposed model than that in [10]. The outage probabilities of lines 22-23, 38-37, 49-50 increase from 0.01 to 1.00 in the proposed model. In [10], these three lines are identified to be outaged at the first stage so that the power flow increment at the first stage can be maximized.  However, this approximation underestimates the interaction between outaged lines and may miss the critical propagation path.
Furthermore, the probability of propagation path of this N -20 contingency is assessed as 8.45 × 10 −15 by the proposed model, while the original probability is 1.09 × 10 −35 and that of the propagation path in [10] is 4.72 × 10 −24 . The computation time for the MILP model in [10] is 226.10 seconds. Although the MILP model in [10] takes less computation time, the probability identified by it still underestimates the propagation probability of this N -20 contingency. This further verifies the effectiveness of the proposed model.

C. COMPARISONS WITH OTHER ALGORITHMS
To further verify the advantages of the proposed MILP model, this section considers another two algorithms that are not  [10] for the N-20 contingency. modeled with MILP for comparisons. The first one is the random chemistry algorithm in [22]. For an N -k contingency S, we prune individual line outages from the remaining outage set until set S contains no outages. In each iteration, the line outage whose removal causes the least increment of outage probability of contingency S is pruned from the current contingency set. The propagation path of contingency S is identified as the reverse sequence of the pruning process. The second is the greedy algorithm, in which we also use k iterations to identify the propagation path. In each iteration, the line outage whose removal causes the maximum increment of outage probability of contingency S is pruned. Finally, the propagation path is the same as the sequence of the pruning process.
The random chemistry algorithm and the greedy algorithm are firstly used to evaluate the N -10 contingency in Section IV.A. The propagation path identified by the random chemistry algorithm is The propagation probabilities of the N -10 contingency in the random chemistry algorithm and the greedy algorithm are 9.58 × 10 −13 and 1.55 × 10 −7 , respectively. Table 5 compares the propagation path from different methods.
Compared with the proposed method and the greedy algorithm, the propagation path in the random chemistry algorithm is rather different. As shown in Table 5, most probabilities of line outages in the random chemistry algorithm are lower than those in the other two methods, except for line . The reason is that line 8-5 is selected by the random chemistry algorithm at the sixth stage of the propagating path, which increases its power from an initial status of 409.59 MW to a highly overloading status of 651.29 MW. It can also be found that the set of line outages at the first five stages of the greedy algorithm is the same as that of the proposed method. Table 5 shows that most line outages in the greedy algorithm and the proposed method have close outage probabilities, except for line 64-65. In the greedy algorithm, line 64-65 is selected at the first stage of the propagation path, which indicates it is not impacted by other line outages. However, it is selected at the fourth stage in the proposed method. As a result, its power flow increases from 444.79 MW to 536.29 MW, and the corresponding outage probability increases from 0.68 to 1.00.
The comparisons are also performed for the N -20 contingency in Section IV.B. The propagation path in the random chemistry algorithm is which has a probability of 4.93 × 10 −20 . And the propagation path in the greedy algorithm is which has a probability of 6.22 × 10 −16 . The difference in evaluation results also verifies the effectiveness of the proposed model for the identification of propagation paths with higher risk.

V. DISCUSSIONS
In this paper, the evaluation of cascading overload failures is improved by the MILP model of the propagation probability. In future studies, this model has two limitations, which deserve further attention.
Firstly, as mentioned in Section II, it does not model the N -k contingency that separates the network. The network separation may cause a power imbalance in sub-grids and make it infeasible for power flow calculations. Future work should focus on the modeling of the power-rebalance constraints to deal with the system separation scenarios. For example, when a power imbalance occurs, generator dispatching or load shedding should be performed.
In addition, as the proposed MILP evaluation model is introduced with the DC power flow constraints, it does not include the bus voltage constraints. However, the voltage violation may activate the under-voltage load shedding or the cascading tripping of renewable generators in some N -k contingencies. Thus, the AC power flow can be considered in the evaluation model and new linearization techniques are required for the construction of MILP model.

VI. CONCLUSION
In this paper, we proposed a multi-stage mixed-integer linear programming to evaluate the propagation probability of N -k contingency as cascading overload failures. The proposed model takes the probability of propagation path of N -k contingency as the objective function. The nonlinear terms of the evaluation model are linearized with different techniques with respect to its properties. Case studies on N -10 and N -20 contingencies of the IEEE 118-bus system verify the effectiveness of the proposed model. The detailed simulation analysis proves that the proposed model can avoid the limitations in the approximation of line outage probability with its increment of power flow. The propagation paths identified by the proposed model can reveal more interdependency and interactions between the multiple line outages in N -k contingencies.
In the future, the proposed model may be further applied in the evaluation framework of high-risk contingencies with high impact and low probability. For example, when a bilevel optimization framework is considered, the probability can be evaluated in the upper-level problem and the impact can be assessed in the lower-level problem. In the upper-level problem, the proposed MILP model can provide a tool to evaluate the probability of potential propagation paths due to the cascading overload failures for contingencies induced by extreme disasters or physical attacks. This can help the power system planner and operator perform a more comprehensive control of reliability risk.