Risk of Cascading Blackouts Given Correlated Component Outages

TL;DR

This paper investigates how correlations in component failures, modeled with copulas, influence the risk of cascading blackouts in power systems, revealing that higher correlation increases blackout risk and the significance of higher order contingencies.

Contribution

It introduces a copula-based method to quantify the impact of correlated component outages on blackout risk, extending analysis to large-scale power system models.

Findings

Correlation increases blackout risk.

Higher order contingencies become more significant with correlation.

New bounding approaches for large systems.

Abstract

Cascading blackouts typically occur when nearly simultaneous outages occur in k out of N components in a power system, triggering subsequent failures that propagate through the network and cause significant load shedding. While large cascades are rare, their impact can be catastrophic, so quantifying their risk is important for grid planning and operation. A common assumption in previous approaches to quantifying such risk is that the $k$ initiating component outages are statistically independent events. However, when triggered by a common exogenous cause, initiating outages may actually be correlated. Here, copula analysis is used to quantify the impact of correlation of initiating outages on the risk of cascading failure. The method is demonstrated on two test cases; a 2383-bus model of the Polish grid under varying load conditions and a synthetic 10,000-bus model based on the geography of the Western US. The large size of the Western US test case required development of new approaches for bounding an estimate of the total number of N-3 blackout-causing contingencies. The results suggest that both risk of cascading failure, and the relative contribution of higher order contingencies, increase as a function of spatial correlation in component failures.

Tables

Table 1. TABLE I: Risk Attributable to N − 2 𝑁 2 N-2 and N − 3 𝑁 3 N-3 Malignancies in the Polish Test Case for Varying Levels of Spatial Correlation.

	$L$ (km)
${\rho }_o$	0	100	200	300
0.00	0.0394	-	-	-
0.05	-	0.0586	0.0675	0.0720
0.10	-	0.0876	0.1142	0.1290
0.15	-	0.1314	0.1918	0.2293

Table 2. TABLE II: Estimated Lower Bounds ( L ​ B 𝐿 𝐵 LB ) and Upper Bounds ( U ​ B 𝑈 𝐵 UB ) on Risk Attributable to N − 2 𝑁 2 N-2 and N − 3 𝑁 3 N-3 Malignancies in the Western US Test Case for Varying Levels of Correlation.

		$L$ (km)
${\rho }_o$		0	100	200	300
0.00	$LB$ $UB$	0.0654 0.0665	-	-	-
0.05	$LB$ $UB$	-	0.0846 0.0864	0.0950 0.0974	0.1019 0.1048
0.10	$LB$ $UB$	-	0.1148 0.1181	0.1444 0.1502	0.1654 0.1735
0.15	$LB$ $UB$	-	0.1631 0.1701	0.2293 0.2445	0.2801 0.3036

Equations

$$R_{\omega }=p_{\omega }{s}_{\omega }$$

$$R_k=\sum _{\omega \in {\Omega }_k}{R}_{\omega }$$

$$\widehat{R_k}=\frac{|{\Omega }_k|}{|{\Omega }_k^{sampled}|}\sum _{\omega \in {\Omega }_k^{sampled}}{R}_{\omega }$$

$$\hat{R}=\sum _{k\in \{\mathrm{2..}{k}_{max}\}}\widehat{R_k}$$

$$F_X(\mathbf{T})=\mathrm{Pr}\left(\bigcap _{i=1}^k{X}_i\le t_i\right)$$

$$F_{X_i}(x_i)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x_i-{\mu }_i}{{\sigma }_i\sqrt{2}}\right)\right]$$

$${\sigma }_i=\frac{-1}{{\mathrm{erf}}^{-1}(2p_i-1)\sqrt{2}}$$

$${\rho }_{ij}={\rho }_o{e}^{-d_{ij}/L}$$

$$\mathrm{cov}(i,j)={\rho }_{ij}{\sigma }_i{\sigma }_j$$

$$f(\mathbf{x})=\frac{1}{\sqrt{(2\pi {)}^k|\mathbf{C}|}}\exp \left\{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\top }{\mathbf{C}}^{-1}(\mathbf{x}-\mathbf{\mu })\right\}$$

$$F_X(0)=\int {\text{displaylimits}}_{-\infty }^0\int {\text{displaylimits}}_{-\infty }^0\cdots \int {\text{displaylimits}}_{-\infty }^{0}f(x_1,x_2,\dots ,x_k)\mathrm{d}{x}_1\mathrm{d}{x}_2\dots \mathrm{d}{x}_k$$

$$\mathrm{Dist}(U,V)=\frac{{\sum }_{i=1}^{2}d(u_i,V)+{\sum }_{i=1}^{2}d(v_i,U)}{4}$$

$$\displaystyle d(v_{i},U)=\left\{\begin{array}[]{ll}||v_{i}-u_{1}||&t\leq 0\\ ||v_{i}-(u_{1}+tm)||&0\leq t\leq 1\\ ||v_{i}-u_{2}||&t\geq 1\end{array}\right.$$

Full text

1 Introduction

Cascading power failures are typically initiated when a small number of $k$ components in a power system of $N$ components disconnect nearly simultaneously, and the subsequent rerouting of power flow triggers additional component outages. This process continues until the system reaches a state of equilibrium. While most cascades do not propagate extensively throughout the network, the rare cases when they do can cause massive blackouts affecting millions of people. Due to their vast size and substantial social and economic costs, the risk they pose to power systems is significant [1, 2, 3].

Networks with heterogeneous load profiles, such as power systems, are particularly prone to cascades; without the right precautions, even a single node may trigger a cascade [4]. To mitigate the risk posed by cascading failure, power systems are required to operate such that no single component outage will cause a cascade (so-called $N$-1 security). While grid planners and operators are now also obligated to consider the risk of cascading failure due to multiple contingencies $(k>1)$ [5], it is not yet clear how to estimate this risk. For brevity, minimal $N-k$ contingencies that result in a cascading blackout are referred to as “malignancies”, while contingencies that do not cause a blackout are referred to as “benign” [6]. By “minimal”, we mean that no smaller subset of outages results in a blackout. Analysis of high order malignancies is challenging due to the nonlinear ways in which cascades propagate, the vast number of $N-k$ malignancies, and the combinatorial search space of possible contingencies.

In addition to helping to quantify risk of cascading failure, studying $N-k$ malignancies may potentially inform mitigating actions. For example, prior research into a simple model of cascading overloads in communication networks [7] suggests that the intentional removal of key components directly after initiating sets of outages may reduce the size of subsequent cascades. In a power system model of the Polish grid, optimally dispatching generation assuming a 50% reduction in line limits on the 3 branches that contribute the most to the risk of cascades from $N-2$ malignancies dramatically reduces the overall risk of cascading failure with only a modest increase in operational costs [8].

Many previous approaches to cascading failure risk analysis (including our own) assumed initiating component outages to be independent events [1, 8, 9, 10, 11]. However, $N-k$ malignancies triggered by the same exogenous event, or “common cause”, represent a significant source of risk to power systems [12], and can result in spatial correlation in initiating outages. For example, extreme weather events can result in spatially correlated damage [13], protection system failures can sometimes cause multiple outages within a small geographic region [14], and terrorist attacks may be spatially localized, such as in the 2013 sniper attack on the Metcalf substation near San Jose, CA, where the perpetrators shot 17 transformers at the same substation [15]. Non-spatial attributes, such as component age, may also induce correlations in component failures [16, 17].

There is a dramatically increasing computational burden to assessing risk for $N-k$ malignancies as $k$ increases, so it is important to understand the degree to which higher-order ($k>2$) malignancies contribute to risk, and thus how important it is to consider them in risk estimation. Even though there are many more $N-3$ than $N-2$ malignancies for any given system, when there is no correlation in initiating outages the probability of $N-3$ malignancies occurring is so much lower than that of $N-2$ malignancies that the impact of $N-3$ malignancies on risk is negligible [9]. However, as correlation in component outages increases, the impact of higher order malignancies on risk will increase. To what degree should risk analysis take into account the conditional probability of component failure, given a common cause?

There has been some prior work on ways to incorporate correlation into risk analysis. In [18], correlation was incorporated by assuming 100% correlation of outages within a fixed radius. In [13], spatial correlation was achieved by determining outage rates of lines adjacent to initial failures probabilistically, according to a Poisson process. In [19], a random field with spatial autocorrelation was used in a cascade model to assess risk from common-cause events. Others [20] have simulated the impact of hidden relay failures on cascading failure risk by allowing proximate lines to trip probabilistically.

Another approach to incorporating correlation into risk estimation is via copula analysis. Popularized in the field of finance [21], copulas have been used in a wide variety of disciplines to model the co-dependence of multiple variables [22, 23]. Within the realm of power systems, copulas are a popular tool for uncertainty analysis. They have been used in the modelling of stochastic generation, such as wind [24, 25, 26]. The impacts of variable infeeds on security assessment have also been considered using copulas [27]. Li [28] suggests copulas as a useful way to incorporate correlation between random variables in power systems risk analysis.

A flexible and generalizable approach to risk estimation given correlated component outages was presented in [29] and used to estimate risk due to $N-2$ malignancies in a 2383-bus model of the Polish grid. This paper extends that work in several significant ways including: (i) incorporating the effects of $N-3$ malignancies, (ii) studying how the risk due to $N-3$, relative to $N-2$, malignancies changes as a function of correlation in outage probabilities, and (iii) applying the method to a much larger and more geographically realistic 10,000-bus test case, which necessitated (iv) development of new methods for estimating the total number of $N-3$ malignancies.

This paper is organized as follows: methods for risk estimation using samples of $N-k$ malignancies, and the computationally efficient “Random Chemistry” (RC) sampling method used in this work, are reviewed in Sections 2.1 and 2.2, respectively. In Section 2.3 a method using copula analysis to incorporate initiating outage correlations into risk estimation is presented, and in Section 2.4 an approach to quantifying distance between transmission lines, when considering spatial correlation, is described. The two test cases used to demonstrate the method are described in Section 2.5; new methods for bounding the total number of $N-3$ malignancies in large systems are described in Section 2.6. Results and discussion are presented in Sections 3 and 4, respectively.

2 Methods

2.1 Estimating Risk of Cascading Failure

This study uses the method for estimating risk of cascading failure from sampled $N-k$ malignancies presented in [11, 9], briefly reviewed below.

The risk due to a set of branches (transmission lines or transformers) $\omega$ can be calculated as  [30]:

[TABLE]

where $p_{\omega}$ is the joint probability of the branches in $\omega$ failing and $s_{\omega}$ is the size of the resultant blackout. Note that $p_{\omega}$ is itself a function of $p_{i}$, the independent outage probability for each branch $i\in\omega$, as well as any effect of correlation among branch outage probabilities (as further defined in Section 2.3).

Blackout size $s_{\omega}$ is quantified as the total power (MW) unserved due to load shedding. In this work, a cascading blackout is considered to have occurred when 5% or more of the total load is shed in DCSIMSEP, a simulator of cascading outages in power systems [31]. The risk posed to the system by all $N-k$ malignancies comprising branches $\omega$, for a given $k$, is then:

[TABLE]

where $\Omega_{k}$ is the complete set of all $N-k$ malignancies for the specified $k$. For realistically-sized power systems it is not computationally tractable to find the entire set $\Omega_{k}$ for $k>2$. However, if $\Omega_{k}^{sampled}\subset\Omega_{k}$ is a large and representative subset of size $|\Omega_{k}^{sampled}|$, comprising all unique $N-k$ malignancies found by many iterations of some sampling strategy, and if the size of the complete set of $N-k$ malignancies $|\Omega_{k}|$ can be estimated, then risk $\hat{R_{k}}$ associated with $N-k$ malignancies, for a given $k$, can be estimated as follows:

[TABLE]

Estimating $|\Omega_{k}|$ for $k>2$ on large systems is itself a very challenging problem, as discussed in Section 2.6.

Considering only malignancies with $k\leq k_{max}$ and assuming that non-minimal supersets of malignancies do not substantially change the amount of load shed, as justified in  [9], the risk of cascading failure can be approximated as:

[TABLE]

In this work, $k_{max}=3$.

2.2 Random Chemistry Sampling

For each $k$, there are ${N\choose k}$ possible $N-k$ contingencies, only a small proportion of which are malignancies. While exhaustive search may be feasible (albeit time consuming) for $k=2$, it is computationally intractable for $k>2$ in large power systems. Thus, many iterations of the Random Chemistry (RC) sampling method were used to efficiently identify large sets of $N-k$ malignancies in each test case. RC is a stochastic set size reduction algorithm for identifying a small minimal set of initiating events that trigger some outcome of interest [32], and was first applied for identifying $N-k$ malignancies in power systems in [6]. For the reader’s convenience, the RC algorithm is briefly reviewed below.

The RC algorithm uses a subset reduction scheme $\{a_{1},a_{2},\dots,a_{final}\}$. Subsets of size $a_{1}$ are randomly sampled from a universal set of $N$ system components until one such set is found that causes a blackout; if $a_{1}$ is relatively large, this typically requires few tries. A set of size $a_{i+1}$ is then randomly sub-sampled from the preceding set of size $a_{i}$ until a set is found that causes a blackout (or some maximum number of sub-samples is tried, in which case the algorithm aborts), and so on for each subsequent set size in the scheme.

If $a_{i+1}=a_{i}/c$, for some constant $c$, then the algorithm requires only $O(N\log N)$ time to identify a subset of size $a_{final}$. As in [6], we use $c=2$ from $a_{1}$ down to subsets of size 20 and then use $c=1.5$ down to $a_{final}$. A bottom-up brute force search of all subsets of a given size $k$ is subsequently applied (conducted in randomized order, starting from $k=2$), exiting when the first minimal malignancy of size $k=\{2,3,\dots,a_{final}\}$ is identified.

Repeated sampling with independent RC trials is performed to compile large subsets of $N-k$ malignancies ($\Omega^{RC}_{k}$). Risk due to each $k\leq k_{max}$ is then calculated using $\Omega^{sampled}_{k}=\Omega^{RC}_{k}$ in (3) for estimating system risk with (4).

A comparison of risk estimation using RC sampling vs. Monte Carlo (MC) sampling on a model of the Polish power system at peak winter load [33] showed that the RC approach was at least two orders of magnitude faster than MC on this system, and did not introduce measurable bias into the estimate [11, 9]. However, these previous studies assumed branch outages were uncorrelated. Under that assumption, $N-3$ malignancies contribute relatively little to the risk, despite the fact that $|\Omega_{3}^{RC}|\gg|\Omega_{2}^{RC}|$, since their probability of occurrence is so much smaller than that of the $N-2$ malignancies [9].

In this study, the universal set is assumed to comprise the set of $N$ branches in each test case, $a_{final}=5$, and up to 20 sub-samples at each set size were allowed before aborting an RC trial, as in [11, 9]. The specific set reduction scheme used for each test case is given in Section 2.5. Note that $k_{max}<a_{final}$ because, with the number of RC trials performed, $|\Omega^{RC}_{k}|$ was insufficient for estimating $|\Omega_{k}|$ for $k>3$.

2.3 Copula Analysis for Correlation

Copula functions couple multivariate distributions to the marginal distributions of individual variables [34]. Given a set of $k$ random variables, $\mathbf{X}=[X_{1},X_{2},\dots,X_{k}]$, where $\Pr(X_{i}\leq t_{i})$ is the marginal probability that branch $i$ fails, for some threshold $t_{i}$, then

[TABLE]

where $\mathbf{T}=\left[t_{1},t_{2},\cdots,t_{k}\right]$, represents the joint probability that $k$ branches fail together. Without loss of generality, it is assumed that $t_{i}=0$ for all $i$.

There are numerous classes of copula functions in popular use. For this demonstration of the method, a Gaussian copula was assumed, but alternative distributions may be assumed where appropriate. Here, it is assumed that the inverse stress on a transmission line $i$ is a univariate Gaussian random variable $X_{i}=\mathcal{N}(\mu_{i},\sigma_{i})$, with mean $\mu_{i}$ and standard deviation $\sigma_{i}$, with the cumulative distribution function:

[TABLE]

Given the independent probability $p_{i}$ of branch $i$ going out, $\mu_{i}$ and $\sigma_{i}$ are chosen such that when $X_{i}<0$, the branch goes out. In other words, $\mu_{i}$ and $\sigma_{i}$ are chosen such that $F_{X_{i}}(0)=p_{i}$ for each branch $i$. Without loss of generality, it is assumed that $\mu_{i}=1,\forall i$, and each $\sigma_{i}$ is then solved for as follows:

[TABLE]

A multivariate normal distribution $\mathbf{X}=\mathcal{N}(\boldsymbol{\mu},\mathbf{C})$ with mean $\boldsymbol{\mu}=[\mu_{1},\mu_{2},\dots,\mu_{k}]$ and covariance matrix $\mathbf{C}$ is then used as the copula function to couple these univariate marginal distributions (Fig. 1).

In this study, it is assumed that the correlation between outages in branches $i$ and $j$ decays exponentially with the distance between them $d_{ij}$, according to:

[TABLE]

where $\rho_{o}$ represents the maximum possible correlation coefficient (at distance zero) and $L$ represents the characteristic length, which controls the decay rate of the correlation; $L$ can be interpreted as the distance at which $\rho_{ij}$ reaches $\rho_{o}/e$ (i.e., $\approx 36.8\%$ of $\rho_{o}$) (Fig. 2).

Eq. (8) can be adjusted to represent a wide range of exogenous common cause events individually, or in combination, by adjusting the parameters $\rho_{o}$ and $L$ to align with data for a particular set of threats. The exponential decay form captures the spatially decaying nature of earthquakes [35], and can approximately capture the impact of other threats that are likely to be geographically correlated, such as tornados [36] and hurricanes [37].

The resulting correlation coefficient $\rho_{ij}$ calculated by (8), and the standard deviations $\sigma_{i}$ and $\sigma_{j}$ calculated by (7), are used to calculate the pairwise covariance between branches $i$ and $j$ as:

[TABLE]

Using (9) to find each element of the covariance matrix $\mathbf{C}$, the probability density function of the multivariate normal distribution (10) is used to form the copula.

[TABLE]

Integrating (10) over the region in the joint distribution that represents outages of all system components gives:

[TABLE]

where $F_{X}(\mathbf{0})$ represents the joint outage probability of $k$ components $\Pr(\mathbf{X}\leq\mathbf{0})$. The multiple-integral in (11) represents the generalized solution for arbitrary $k$ and is equivalent to the cumulative distribution function of the multivariate normal distribution, which can be solved efficiently in MATLAB using methods described in [38]. In this work, $k\in\{2,3\}$.

2.4 Defining Inter-Branch “Distance”

The definition of “distance” will vary based on the type of common cause threatening the system. Assuming spatial correlation in branch outages here, without loss of generality a modified version of the inter-branch distance metric defined in [29] is used. Branches are assumed to be straight lines between the buses that form their endpoints. Given branch $U$ with endpoints ($u_{1},u_{2}$) and branch $V$ with endpoints ($v_{1},v_{2}$), let the distance from $U$ to $V$ be defined as

[TABLE]

where $d(v_{i},U)$ is the minimum Euclidean distance from the point $v_{i}$ to the line segment $U=(u_{1},u_{2})$, calculated as:

[TABLE]

where $m=u_{2}-u_{1}$ and t=\frac{(v_{i}-u_{1})\,\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}\,m}{||m||^{2}}, as illustrated in Figure 3.

This formulation defines a semi-metric since the triangle inequality does not hold in some cases. However, all other formal requirements of a metric are met. Specifically:

$Dist(U,V)\geq 0$ 2. 2.

$Dist(U,V)=Dist(V,U)$ 3. 3.

$Dist(U,V)=0\iff U=V$

The third identity implies branches $U$ and $V$ share the same endpoints, thus are parallel. This $Dist(U,V)$ measure is consistent with what would be intuitively expected when considering spatially correlated damage. For example, in Fig. 4, ${Dist}(A,B)>{Dist}(C,D)$ and ${Dist}(E,F)>{Dist}(G,H)$.

Using the $Dist(U,V)$ measure, it is apparent that branch pairs that form malignancies are much closer together than those of benign contingency pairs in both test cases (Fig. 5). This property will exacerbate the effects of spatial correlation on risk of cascading failure.

2.5 Case Studies

This risk estimation approach is demonstrated on two publicly available test cases, modeling the Polish and Western United States (US) transmission systems.

The Polish test case, examined in previous work on risk estimation [8, 9, 29], contains 2383 buses and 2896 branches at a peak winter load and is distributed with the MATPOWER simulation package [33]. The true spatial locations of branches and buses are not publicly available for this test case, so hypothetical locations were inferred based on a graph layout of the grid topology, assuming branches are straight lines between buses (Fig. 6). This layout was then scaled to $670\times 670$ km, the approximate width/height of Poland, to simulate geographic distances. Some of the transmission lines were overloaded in the Polish test case provided by [33], so the adjusted base case described in [8] was used. Unless otherwise stated, references to the “Polish test case” refer to this adjusted base case. As in [9], different load levels were modeled in the Polish test case from 80% to 115% of the adjusted base case by multiplying all line loads by a scalar factor and then re-running the security constrained optimal power flow, to ensure the pre-contingency system at each load level is $N-1$ secure.

The Western US test case is a synthetic network based on the footprint of the western Unites States and comes via the Electric Grid Test Case Repository [39]. This test case is much larger than the Polish test case, with 10,000 buses and 12,706 branches, and has a more realistic geographic layout (Fig. 7). As with the Polish test case, some transmission lines were overloaded for the Western US test case provided by [39], and so adjustments were made as described in [8]. Since the case did not include short and long-term emergency flow limits (“RateB” and “RateC”), they were synthesized to be 110% and 150% of normal (“RateA”) limits, respectively.

Independent branch outage rates were not available for either the Polish or Western US test cases. For the results presented here, all independent outage rates were set equal to the mean outage rate of 0.9158 hours per year provided by the RTS-96 test case [40]. These independent outage rates were deliberately assumed identical for all branches in order to more clearly elucidate the impact of spatial correlations in outage rates, as assessed using (4), for all combinations of $L\in\{0,100,200,300\}$ km and $\rho_{o}\in\{0.00,0.05,0.10,0.15\}$.

For the results shown here, the RC algorithm used the subset reduction scheme $\{80,40,20,14,10,7,5\}$ for the Polish model, as in [6, 11, 8, 9, 29]. For the larger Western US model, the initial RC subset size $a_{1}$ was raised to 320 to increase the probability that the initial subset causes a blackout; thus, the Western US test case used the subset reduction scheme $\{320,160,80,40,20,14,10,7,5\}$. We did 1,000,000 RC trials for the Polish test case and 704,400 RC trials for the Western US model; fewer RC trials were used for the larger test case because computation time was much higher than for the Polish model (averaging 9.5 seconds per RC trial for the Western US model vs. 2.35 seconds for the Polish model, on an Intel Core i5-3470 CPU @ 3.2GHz with 8 GB of RAM).

2.6 Estimating $|\Omega_{k}|$

As described in Section 2.1, this approach to risk estimation requires an estimate of the total number of $N-k$ malignancies $|\Omega_{k}|$, for $k\leq k_{max}$. There was no need to estimate $|\Omega_{2}|$, since RC sampling identified the complete set of $N-2$ malignancies $\Omega_{2}$ in both test cases, as evidenced by the flattening in the accumulation curves (Fig. 8 (top)), and later verified through brute force search for the Polish test case. The set of unique $N-2$ malignancies was complete after only 5,090 non-unique $N-2$ malignancies had been found by RC sampling in the Polish test case (of 4,191,960 possible $N-2$ contingencies) and after only 9,364 non-unique $N-2$ malignancies had been found by RC sampling in the Western US test case (of 80,714,865 possible $N-2$ contingencies).

However, obtaining the entire set of N-3 malignancies is not computationally tractable in either test case, due to the sheer size of these sets. It was initially argued (incorrectly) in [6] that, if one has already identified $i$ of the $N-k$ malignancies using independent RC trials, then the probability that the next identified $N-k$ malignancy has not previously been found is $(|\Omega_{k}|-i)/|\Omega_{k}|$, so one could infer $|\Omega_{k}|$ from the observed frequency with which unique malignancies were found (assuming sufficient curvature in the accumulation curve). However, the assumption that independent RC trials uniformly sample from the $\Omega_{k}$ sets has since proven false. In subsequent studies it was discovered that the accumulation curves were not exponential (as they would be if the sampling were uniform), but could be better fit with a 4-parameter exponential Weibull curve to estimate $|\Omega_{k}|$ [9]. While this non-linear curve-fitting approach works for estimating $|\Omega_{3}|$ in the Polish test case, the Western US test case is so much larger that there is insufficient curvature in the $N-3$ accumulation curve (Fig. 8 (bottom)) to reliably fit a curve.

It has previously been noted that the frequency of occurrence of individual branches in $N-2$ malignancies is heavy-tailed [6, 41]. A similarly heavy-tailed distribution is apparent in the distribution of occurrences of specific branch pairs in $N-3$ malignancies, in both the Polish and Western US (Fig. 9).

Further examination reveals that the set reduction scheme used in RC does, indeed, introduce a sampling bias when sampling from such heavy-tailed distributions. Specifically, the branch pairs that appear in disproportionately large numbers of $N-3$ malignancies are systematically under-sampled by the RC algorithm. To illustrate this, the 20 most frequently occurring branch pairs found in $|\Omega^{RC}_{3}|$ for the Western US test case were selected, and a brute force search of all possible $N-3$ contingencies that included each of these top 20 branch pairs (requiring $O(N)$ computation time for each branch pair) was performed. In Fig. 10, the proportion of $N-3$ malignancies found by RC that contain each of these branch pairs is plotted as a function of the observed number of occurrences of the branch pairs in $\Omega^{RC}_{3}$. A clear negative trend is present, with $N-3$ malignancies containing the most frequently occurring branch pairs severely under-sampled relative to $N-3$ malignancies containing less frequently occurring branch pairs. While a thorough explanation of the causes of this sampling bias are beyond the scope of this paper, here the bias is exploited to estimate both lower and upper bounds on $|\Omega_{3}|$. (It is worth noting that, as $|\Omega^{RC}_{3}|$ approaches $|\Omega_{3}|$, the sampling bias of branch pairs found in $N-3$ malignancies decreases. However, for large networks such as the Western US test case, it is not computationally feasible to sample more than a small fraction of the $N-3$ malignancies, so the sampling bias remains high.)

Given that sampling probabilities are unequal, this problem is analogous to the common conservation biology task of estimating population sizes via mark-and-recapture surveys in closed populations with heterogeneous sampling probabilities. There are numerous techniques that have been developed for this kind of problem [42]. Here, Chao’s method [43] is used, because it is known to be particularly robust to heterogeneous sampling probabilities. In the power system context, the “population” under consideration is $\Omega_{k}$, the set of all $N-k$ malignancies. Chao’s estimate is calculated as $|\Omega_{3}|^{Chao}=|\Omega_{3}^{RC}|+n_{1}^{2}/(2n_{2})$, where $n_{1}$ is the number of $N-3$ malignancies found exactly once by RC sampling and $n_{2}$ is the number of $N-3$ malignancies found exactly twice by RC sampling. Chao’s method produces a lower-bound on the population size within a fixed confidence interval [43], so it is assumed that $|\Omega_{3}|^{Chao}\leq|\Omega_{3}|$.

An upper-bound on $|\Omega_{3}|$ can be estimated by taking advantage of the two observations demonstrated above: (i) certain branch pairs appear disproportionately often in $N-3$ malignancies (Fig. 9), and (ii) the most frequently occurring branch pairs are under-sampled, relative to less frequent branch pairs (Fig. 10). Based on these observations, the Random Chemistry Proportional (RCP) method is proposed as a way to estimate an upper bound on $|\Omega_{3}|$, as follows: (i) apply RC sampling for a sufficient number of trials such that the identity of the most frequently occurring branch-pair ($Pair_{max}$) in the $N-3$ malignancies of the growing set $\Omega^{RC}_{3}$ becomes stable (for the Western US test case, $Pair_{max}$, indicated by the star in Fig. 10, became stable after about 7000 RC trials), (ii) perform a brute force search of all possible $N-3$ contingencies that include $Pair_{max}$ (this requires only $O(N)$ simulations) to determine the true number of $N-3$ malignancies that include $Pair_{max}$, (iii) compute what proportion of all $N-3$ malignancies that include $Pair_{max}$ were found by RC sampling, and finally (iv) we estimate the total number of $N-3$ malignancies (referred to as $|\Omega_{3}|^{RCP}$) by assuming that all other less-frequently occurring branch-pairs have found this same proportion of the total number of $N-3$ malignancies in which they occur. Assuming that $Pair_{max}$ is under-sampled, this method provides an overestimate, and hence an upper-bound, on $|\Omega_{3}|$; i.e., it is expected that $|\Omega_{3}|^{RCP}>|\Omega_{3}|$.

As the number of $N-3$ malignancies found by RC sampling increases, $|\Omega_{3}|^{Chao}$ and $|\Omega_{3}|^{RCP}$ can be seen to be converging (Fig. 11), thus increasing the confidence in these as lower and upper bounds on the true value of $|\Omega_{3}|$. Risk estimates are calculated for the rightmost values of $|\Omega_{3}|^{Chao}$ and $|\Omega_{3}|^{RCP}$ shown in Fig. 11, to obtain approximate bounds on risk due to $N-3$ malignancies for the Western US test case.

Similar approaches could conceivably be applied for estimating $|\Omega_{k}|$ for $k>3$, however in this study $|\Omega^{RC}_{4}|$ and $|\Omega^{RC}_{5}|$ were insufficient to support this.

3 Results

3.1 Set sizes of $\Omega_{2}$ and $\Omega_{3}$

Brute force search was used to verify that RC sampling found all $N-2$ malignancies in the Polish base test case, with $|\Omega_{2}|=540$. It is assumed that RC sampling also found all $N-2$ malignancies in the Western US test case with $|\Omega_{2}|=564$, since the accumulation curve became flat (Fig. 8(top)). Using the non-linear curve-fitting method of [9], $|\Omega_{3}|$ is estimated to be $\approx 6.4\times 10^{4}$ in the Polish test case at the base load. Using the Chao lower-bounding method [43] and the RCP upper-bounding method (described in Section 2.6), it is estimated that $2.0\times 10^{5}\leq|\Omega_{3}|<2.9\times 10^{5}$ in the Western US test case.

3.2 Impact of $N-2$ Correlation and Load Level on Risk

As shown in prior work [8, 9], the load levels on the Polish grid can greatly affect the vulnerability of the network to cascading power failure due to $N-2$ malignancies. As noted in [9], risk varies non-linearly and non-monotonically with load, in part due to variations in the proximity of generation to demand that result from optimal power flow dispatch at different load levels. Risk actually tends to drop at very high load levels because the security constrained optimal power flow results in more local generation, thus reducing the flow on critical long-distance transmission lines that can participate in many $N-k$ malignancies when overloaded. For a direct comparison to the results presented in [9], the impact of spatial correlation in $N-2$ malignancies on risk was assessed as a function of load in the Polish test case.

Changes in the system risk as a function of load at $L=300$ km (the longest characteristic correlation length tested) for 3 values of $\rho_{o}$ are illustrated in Fig. 12. Risk increased faster than linearly as a function of linearly increasing $\rho_{o}$, at each of the given load levels. The largest percentage increase in system risk occurred at load level 114%, while the greatest absolute increases in risk occur between load levels of 97%-112%, where there are the most $N-2$ malignancies. In general, while introducing correlation in initiating outages magnifies the risk of cascading blackouts, it does not fundamentally alter the overall shape of the risk curve as a function of load at $L=300$ km.

When $\rho_{o}=0.15$ (the largest $\rho_{o}$ tested) and $L$ was varied from 0 to 300 km, results were superficially similar to those in Fig. 12, in that higher correlation increases risk without changing the overall shape of the risk curve as a function of load. As was the case when $L$ was fixed, the largest percentage increase in system risk was found to occur at load level 114%, and greatest absolute increases in risk occured between load levels of 97%-112%. However, in this case risk increases slower than linearly with linear increases in $L$, with the largest increases occurring for intermediate values of $L$ (Fig. 13). This occurs because increasing $L$ beyond a certain point has diminishing impact on correlation, as $L$ approaches the radius of the network.

The super-linear increases in risk as a function of $\rho_{o}$ and sub-linear increases in risk as a function of $L$ at the base load are clearly illustrated in Fig. 14.

3.3 Risk from $N-2$ and $N-3$ Malignancies

Risk of cascading blackouts posed by $N-2$ and $N-3$ malignancies was computed for both the Polish base test case and the Western US test case, over all values of $L$ and $\rho_{o}$ tested, using the set size estimates given in Sec. 3.1.

For the Polish test case (Table I), the increase in estimated risk due to spatial correlation ranges from 149% for the most modest level of correlation tested ($L=100$ km, $\rho_{o}=0.05$) to 582% in the most extreme case tested ($L=300$ km, $\rho_{o}=0.15$), relative to the uncorrelated case.

For the Western US test case, (Table II), the increase in lower (upper) bounds on risk estimates varied from 129% (130%) for the most modest level of correlation tested ($L=100$ km, $\rho_{o}=0.05$) to 428% (456)% in the most extreme case tested ($L=300$ km, $\rho_{o}=0.15$), relative to the uncorrelated case.

For both test cases, the general effect of $L$ and $\rho_{o}$ on risk that is described in Section 3.2 also holds in these results. That is, risk tends to grow faster than linearly with respect to $\rho_{o}$ and slower than linearly with respect to $L$. The larger proportionate increases in the Polish test case, relative to the Western US test case, occur because the average distance between branches in malignancies are shorter than in the Western US test case (Fig. 5), thus magnifying the impacts of spatial correlation.

3.4 Relative Risk of $\boldsymbol{N-2}$ vs. $\boldsymbol{N-3}$ Malignancies

It is expected that $N-3$ malignancies will contribute more to risk when there is spatial correlation in initiating outages, but it is not clear to what degree. There are several factors that could potentially disproportionately affect the impact of $N-3$ malignancies on risk when there is spatial correlation, relative to that of $N-2$ malignancies, including: (i) size of blackouts caused by $N-3$ vs. $N-2$ malignancies; (ii) the independent probability of branch outages in $N-3$ vs. $N-2$ malignancies; (iii) the distance between branches in $N-3$ vs. $N-2$ malignancies. These factors are each discussed in more detail below.

3.4.1 Blackout Sizes

If the sizes of blackouts resulting from $N-3$ malignancies were larger than $N-2$ malignancies, this could disproportionately increase the relative contribution of $N-3$ malignancies to risk when spatial correlation is present. However, we have observed that the sizes of cascading blackouts tend to follow similarly shaped distributions, independent of the number of component outages in the triggering event, due to similar patterns of network separation. This is illustrated by the distributions of blackout sizes (as estimated by DCSIMSEP) caused by all $N-2$ malignancies, and the subsets of identified $N-3$ malignancies found by RC sampling, for both the Polish and Western US test cases (Fig. 15). In both test cases, the median blackout size for the identified $N-3$ malignancies was actually lower than those caused by the $N-2$ malignancies. Specifically, for the Polish test case, the median blackout size caused by $N-2$ malignancies was 7,624 MW vs. 3,372 MW for those caused by identified $N-3$ malignancies. In the Western US test case, the median blackout size from $N-2$ malignancies was 10,473 MW whereas from the $N-3$ malignancies it was 10,382 MW. These patterns and trends continue in the identified sets of $N-4$ and $N-5$ malignancies (Fig. 15).

3.4.2 Independent Branch Outage Rates

In this study, independent outage rates were assumed to be homogeneous for all branches. However, in a real system the distribution of independent outage rates will be heterogeneous. If branches that are typically involved in $N-3$ malignancies are independently more likely to fail than those involved in $N-2$ malignancies, this could inflate the relative risk of $N-3$ malignancies when spatial correlation is present. While there is no obvious rationale for why this might be true, the observation that branches that occur frequently in $N-2$ malignancies also appear frequently in $N-3$ malignancies is indirect evidence against this. For example, in the Polish network, 8 of the 10 most frequently occurring branches in $N-2$ and $N-3$ malignancies are shared, accounting for 44% and 24% of all $N-2$ and $N-3$ malignancies found, respectively. Likewise, for the Western US test case, 9 of the 10 most frequently occurring branches in $N-2$ malignancies are also in the top 10 most frequently occurring $N-3$ malignancies, accounting for 49% and 29% of all $N-2$ and $N-3$ malignancies, respectively.

3.4.3 Distance Between Branches

The distance between branches in $N-3$ vs. $N-2$ malignancies will obviously impact the degree to which spatial correlation will increase their relative contributions to risk. In both the Polish and Western US test cases, median distances between all pairs of branches occurring in identified $N-k$ malignancies increases with $k\in\{2,3,4,5\}$ (Fig. 16). Specifically, in the Polish test case the median distances between pairs of branches were 76.1 km and 121.8 km, in $N-2$ and $N-3$ malignancies, respectively; in the Western US test case the medians were 169.4 km and 494.6 km in the $N-2$ and $N-3$ malignancies, respectively. This helps to mitigate the increase in relative risk from $N-3$ vs. $N-2$ malignancies that occurs as a result of spatial correlation.

3.4.4 Comparing Relative Risk with Correlation

For the Polish test case, $<1\%$ of risk can be attributed to $N-3$ malignancies when there is no correlation whereas under the highest level of correlation considered ($L=300$ km, $\rho_{o}=0.15$), the share of risk associated with $N-3$ malignancies rises to around 9% (Fig. 17). Similarly, for the Western US test case, $N-3$ malignancies account for 3%-5% of risk when there is no correlation, but between 16%-24% under the maximal correlation ($L=300$ km, $\rho_{o}=0.15$) considered (Fig. 18).

4 Discussion

Previous research into cascading failure risk demonstrated that $N-3$ malignancies constitute a relatively low proportion of risk compared to $N-2$ malignancies, assuming initiating branch outage independence [9]. This suggests that, if initiating outages are generally caused by independent events, limiting risk analysis to the more computationally tractable $N-2$ malignancies may be sufficient to capture the majority of risk. However, in reality, common causes such as relay failures, weather disturbances, earthquakes, fire, or spatially localized terrorist attacks may trigger multiple near-simultaneous outages in geographic proximity that could potentially result in cascading blackouts. For these cases an assumption of independence will under-estimate cascading blackout risk.

This paper presents a method that uses copula analysis as a flexible, customizable approach for incorporating correlation into risk calculations, building on preliminary work presented in [29]. The impact of spatial correlation in $N-2$ and $N-3$ initiating outages on risk of cascading blackouts is assessed in the Polish test case as well as the much larger, and geographically more realistic, Western US test case.

The Western US test case has over four times as many branches as the Polish test case, and the increased computational cost of performing cascade simulations combined with the substantially larger search space of $N-3$ contingencies rendered previously developed methods [6, 9] for estimating the set size of $N-3$ malignancies ineffective. Thus, extending the approach to include risk from $N-3$ malignancies in the Western US test case required new methods to estimate lower- and upper-bounds on the total number of $N-k$ malignancies, for $k=3$.

The results indicate that when spatial correlation is present in initiating outages, the relative contribution of $N-3$ malignancies to risk of cascading blackouts increases, although the increase is partially mitigated by the fact that pairwise distances between branches in $N-3$ malignancies are greater than in $N-2$ malignancies.

It is expected that the impact of even higher-order malignancies will similarly increase with increasing spatial correlation, even though median pairwise distances between branches in malignancies continue to increase with $k$. In principle, the approaches to estimating lower and upper bounds on $|\Omega_{3}|$ presented here for the Western US test case could also be applied for estimating $|\Omega_{k}|$ for $k>3$, given a sufficient number of RC trials. We are currently exploring these and other methods on large synthetic networks to establish their limitations as a function of network size and functional heterogeneity. While ignoring higher-order $N-k$ malignancies when component outages are correlated will likely underestimate the magnitude of the risk of cascading failures, estimating $|\Omega_{k}|$ for $k$ higher than 3 (or possibly 4) on large networks may not be computationally tractable, due to the sheer number of these high-order malignancies. Preliminary work indicates that, for the Western US test case, $|\Omega_{4}|$ may be at least two orders of magnitude higher than $|\Omega_{3}|$.

The lack of accurate data regarding independent transmission line outage rates and the impacts of common cause events on those rates are an important limitation to applying these methods in practice. Some such data is available to industry through systems such as the NERC TADS database, but these data are typically unavailable for research purposes. Increasingly, efforts are being made to better predict how common-cause events such as weather-related events will impact the grid [44]. Such knowledge will inform specific applications of the generalized framework introduced herein.

The methods presented in this paper should work with any simulator (AC, DC, or even something more sophisticated like a full dynamics cascading failure simulator). However, more complicated simulation models require much larger input datasets and tuning these models to get accurate results is a longer process (see [45, 46] for illustrations of the challenges associated with dynamic modeling of cascading failure). For example, in order to get an accurate model from an AC or dynamic power system simulator one would need to accurately model all of the dynamic reactive power elements, such as synchronous condensers and switched capacitor banks, in order to get accurate results. The impact of these controls is that a system with large amounts of reactive support will act more like a DC model with uniform voltage, than an AC model without reactive support. Since the focus of this paper is primarily on the computational method for risk analysis, rather than precise power systems details, we have used a simulator based on the DC power flow. Based on our experience with more complicated simulation models, we do not expect that a more complicated simulator would produce qualitative differences in the results, although quantitative differences would result since there are more mechanisms of cascading in an AC power flow model.

Future work will study the impact of parametric choices such as distance metrics, correlation functions, and the underlying branch outage probability distributions. The method will also be applied to study the risk of cascading failure in interdependent networks, which are ubiquitous in human-engineered infrastructures [47]. For example, coupling between communication and power networks can substantially impact their robustness to cascading failures [48, 49]; the method presented herein could help to quantify the impact of this coupling on risk in the presence of correlated component outages. As new methods for measuring the risk of cascading failure in systems with correlated initiating event probabilities emerge, there will be a need for comparisons to better understand the relative computational efficiency and accuracy of these approaches.

Acknowledgments

This work was supported in part by the National Science Foundation, award numbers ECCS-1254549, CNS-1735513 and DGE-1144388. The computational resources provided by the University of Vermont Advanced Computing Core (VACC) are gratefully acknowledged.