Vulnerability and co-susceptibility determine the size of network cascades

TL;DR

This paper introduces a statistical framework that predicts the size of network cascades by analyzing component vulnerability and co-susceptibility, with power grids as a key example, offering insights applicable to various complex systems.

Contribution

The paper develops a novel framework that incorporates vulnerability and co-susceptibility to predict cascade sizes, revealing structured groups of co-failure in networks.

Findings

Correlations between component failures form large co-susceptible groups.

Structured co-susceptibility significantly influences cascade sizes.

Framework applicable to diverse complex systems beyond power grids.

Abstract

In a network, a local disturbance can propagate and eventually cause a substantial part of the system to fail, in cascade events that are easy to conceptualize but extraordinarily difficult to predict. Here, we develop a statistical framework that can predict cascade size distributions by incorporating two ingredients only: the vulnerability of individual components and the co-susceptibility of groups of components (i.e., their tendency to fail together). Using cascades in power grids as a representative example, we show that correlations between component failures define structured and often surprisingly large groups of co-susceptible components. Aside from their implications for blackout studies, these results provide insights and a new modeling framework for understanding cascades in financial systems, food webs, and complex networks in general.

Equations

$$N=\sum _{\ell =1}^n{F}_{\ell },$$

$${\gamma }_{\mathcal{I}}\equiv \frac{⟨N_{\mathcal{I}}|N_{\mathcal{I}}\ne 0⟩-{\bar{n}}_{\mathcal{I}}}{m-{\bar{n}}_{\mathcal{I}}}>{\gamma }_{\text{th}},$$

$${\mathcal{C}}_{\ell {\ell }^{\prime }}=\{\begin{array}{ll}{\hat{C}}_{\ell {\ell }^{\prime }}^{(\text{sig})}& \text{if ℓ,ℓ′∈Rk for some k≤m′,}\\ 0& \text{otherwise.}\end{array}$$

Full text

The stability of complex networks is largely determined by their ability to operate close to equilibrium—a condition that can be compromised by relatively small perturbations that can lead to large cascades of failures. Cascades are responsible for a range of network phenomena, from power blackouts hines2009large and air traffic delay propagation fleurquin2013 to secondary species extinctions sahasrabudhe2011rescuing ; eco:11 and large social riots watts2002 ; brummitt2015 . Evident in numerous previous modeling efforts bak_1987 ; kinney_2005 ; buldyrev_2010 ; watts2002 ; goh_2003 ; motter_2004 ; dobson2007complex ; brummitt2015 ; sahasrabudhe2011rescuing ; eco:11 is that dependence between components is the building block of these self-amplifying processes and can lead to correlations among eventual failures in a cascade.

A central metric characterizing a cascade is its size. While the suitability of a size measure depends on the context and purpose, a convenient measure is the number of network components (nodes or links) participating in the cascade (e.g., failed power lines, delayed airplanes, extinct species). Since there are many known and unknown factors that can affect the details of cascade dynamics, the main focus in the literature has been on characterizing the statistics of cascade sizes rather than the size of individual events. This leads to a fundamental question: what determines the distribution of cascade sizes?

In this Letter, we show that cascading failures (and hence their size distributions) are often determined primarily by two key properties associated with failures of the system components: the vulnerability, or the failure probability of each component, and the co-susceptibility, or the tendency of a group of components to fail together. The latter is intimately related to pairwise correlations between failures, as we will see below. We provide a concrete algorithm for identifying groups of co-susceptible components for any given network. We demonstrate this using the representative example of cascades of overload failures in power grids (Fig. 1). Based on our findings, we develop the co-susceptibility model—a statistical modeling framework capable of accurately predicting the distribution of cascade sizes, depending solely on the vulnerability and co-susceptibility of component failures.

We consider a system of $n$ components subject to cascading failures, in which a set of initial component failures can induce a sequence of failures in other components. Here we assume that the initial failures and the propagation of failures can be modeled as stochastic and deterministic processes, respectively (although the framework also applies if the propagation or both are stochastic). Thus, the cascade size $N$, defined here as the total number of components that fail after the initial failures, is a random variable that can be expressed as

[TABLE]

where $F_{\ell}$ is a binary random variable representing the failure status of component $\ell$ (i.e., $F_{\ell}=1$ if component $\ell$ fails during the cascade, and $F_{\ell}=0$ otherwise). While the $n$ components may be connected by physical links, a component may fail as the cascade propagates even if none of its immediate neighbors have failed brummitt2015 ; Witthaut:2015 ; dobson_2016 . For example, in the case of cascading failures of transmission lines in a power grid, the failure of one line can cause a reconfiguration of power flows across the network that leads to the overloading and subsequent failure of other lines away from the previous failures dobson_2016 ; anghel2007stochastic .

A concrete example network we analyze throughout this Letter using the general setup above is the Texas power grid, for which we have 24 snapshots, representing on- and off-peak power demand in each season of three consecutive years. Each snapshot comprises the topology of the transmission grid, the capacity threshold of each line, the power demand of each load node, and the power supply of each generator node (extracted from the data reported to FERC FERC ). For each snapshot we use a physical cascade model to generate $K=5{,}000$ cascade events. In this model (which is a variant of that in Ref. anghel2007stochastic with the power re-balancing scheme from Ref. Hines2011 ), an initial perturbation to the system (under a given condition) is modeled by the removal of a set of randomly selected lines. A cascade following the initial failures is then modeled as an iterative process. In each step, power flow is redistributed according to Kirchhoff’s law and might therefore cause some lines to be overloaded and removed (i.e., to fail) due to overheating. The temperature of the transmission lines is described by a continuous evolution model and the overheating threshold for line removal is determined by the capacity of the line anghel2007stochastic . When a failure causes part of the grid to be disconnected, we re-balance power supply and demand under the constraints of limited generator capacity Hines2011 . A cascade stops when no more overloading occurs, and we define the size $N$ of the cascade as the total number of removed lines (excluding the initial failures). This model source-code , accounting for several physical properties of failure propagation, sits relatively high in the hierarchy of existing power-grid cascade models opaModel ; henneaux2016 ; dobson2012_vulnerability ; hines_2015 , which ranges from the most detailed engineering models to simplest graphical or stochastic models. The model has also been validated against historical data Yang:2016 .

In general, mutual dependence among the variables $F_{\ell}$ may be necessary to explain the distribution of the cascade size $N$. We define the vulnerability $p_{\ell}\equiv\langle F_{\ell}\rangle$ of component $\ell$ to be the probability that this component fails in a cascade event (including events with $N=0$). If the random variables $F_{\ell}$ are uncorrelated (and thus have zero covariance), then $N$ would follow Poisson’s binomial distribution Wang:1993 , with average $\tilde{\mu}=\sum_{\ell}p_{\ell}$ and variance $\tilde{\sigma}^{2}=\sum_{\ell}p_{\ell}(1-p_{\ell})$. However, the actual variance $\sigma^{2}$ of $N$ observed in the cascade-event data is significantly larger than the corresponding value $\tilde{\sigma}^{2}$ under the no-correlation assumption for all $24$ snapshots of the Texas power grid (with the relative difference, $\bar{\sigma}^{2}\equiv(\sigma^{2}-\tilde{\sigma}^{2})/\tilde{\sigma}^{2}$, ranging from around $0.18$ to nearly $39$). Thus, the mutual dependence must contribute to determining the distribution of $N$ in these examples.

Part of this dependence is captured by the correlation matrix $C$, whose elements are the pairwise Pearson correlation coefficients among the failure status variables $F_{\ell}$. When the correlation matrix is estimated from cascade-event data, it has noise due to finite sample size, which we filter out using the following procedure. First, we standardize $F_{\ell}$ by subtracting the average and dividing it by the standard deviation. According to random matrix theory, the probability density of eigenvalues of the correlation matrix computed from $K$ samples of $T$ independent random variables follow the Marchenko-Pastur distribution mehta2004random , $\rho(\lambda)=K\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}/(2\pi\lambda T)$, where \lambda_{\pm}=\bigl{[}1\pm\sqrt{T/K}\,\bigr{]}^{2}. Since those eigenvalues falling between $\lambda_{-}$ and $\lambda_{+}$ can be considered contributions from the noise, the sample correlation matrix $\widehat{C}$ can be decomposed as $\widehat{C}=\widehat{C}^{(\text{ran})}+\widehat{C}^{(\text{sig})}$, where $\widehat{C}^{(\text{ran})}$ and $\widehat{C}^{(\text{sig})}$ are its random and significant parts, respectively, which can be determined from the eigenvalues and the associated eigenvectors MacMahon:2015 . In the network visualization of Fig. 1(a), we show the correlation coefficients $\widehat{C}_{\ell\ell^{\prime}}^{\text{(sig)}}$ between components $\ell$ and $\ell^{\prime}$ estimated from the cascade-event data for the Texas grid under the 2011 summer on-peak condition. Note that we compute correlation only between those components that fail more than once in the cascade events. As this example illustrates, we observe no apparent structure in a typical network visualization of these correlations. However, as shown in Fig. 1(b), after repositioning the nodes based on correlation strength, we can identify clusters of positively and strongly correlated components—those that tend to fail together in a cascade.

To more precisely capture this tendency of simultaneous failures, we define a notion of co-susceptibility: a given subset of $m$ components $\mathcal{I}\equiv\{\ell_{1},\ldots,\ell_{m}\}$ is said to be co-susceptible if

[TABLE]

where $N_{\mathcal{I}}\equiv\sum_{j=1}^{m}F_{\ell_{j}}$ is the number of failures in a cascade event among the $m$ components, $\langle N_{\mathcal{I}}|N_{\mathcal{I}}\neq 0\rangle$ denotes the average number of failures among these components given that at least one of them fails, \bar{n}_{\mathcal{I}}\equiv\sum_{j=1}^{m}p_{\ell_{j}}/\bigl{[}1-\prod_{k=1}^{m}(1-p_{\ell_{k}})\bigr{]}\geq 1 is the value $\langle N_{\mathcal{I}}|N_{\mathcal{I}}\neq 0\rangle$ would take if $F_{\ell_{1}},\ldots,F_{\ell_{m}}$ were independent. Here we set the threshold in Eq. (2) to be $\gamma_{\text{th}}=\sigma_{N_{\mathcal{I}}}/(m-\bar{n}_{\mathcal{I}})$, where \sigma_{N_{\mathcal{I}}}^{2}\equiv\sum_{j=1}^{m}p_{\ell_{j}}(1-p_{\ell_{j}})/\bigl{[}1-\prod_{k=1}^{m}(1-p_{\ell_{k}})\bigr{]}-\bar{n}_{\mathcal{I}}^{2}\prod_{k=1}^{m}(1-p_{\ell_{k}}) is the variance of $N_{\mathcal{I}}$ given $N_{\mathcal{I}}\neq 0$ for statistically independent $F_{\ell_{1}},\ldots,F_{\ell_{m}}$. By definition, the co-susceptibility measure $\gamma_{\mathcal{I}}$ equals zero if $F_{\ell_{1}},\ldots,F_{\ell_{m}}$ are independent. It satisfies $-(\bar{n}_{\mathcal{I}}-1)/(m-\bar{n}_{\mathcal{I}})\leq\gamma_{\mathcal{I}}\leq 1$, where the (negative) lower bound is achieved if multiple failures never occur and the upper bound is achieved if all $m$ components fail whenever one of them fails. Thus, a set of co-susceptible components are characterized by significantly larger number of simultaneous failures among these components, relative to the expected number for statistically independent failures. While $\gamma_{\mathcal{I}}$ can be computed for a given set of components, identifying sets of co-susceptible components in a given network from Eq. (2) becomes infeasible quickly as $n$ increases due to combinatorial explosion.

Here we propose an efficient two-stage algorithm for identifying co-susceptible components source-code . The algorithm is based on partitioning and agglomerating the vertices of the auxiliary graph $G_{0}$ in which vertices represent the components that fail more than once in the cascade-event data, and (unweighted) edges represent the dichotomized correlation between these components. Here we use $\widehat{C}^{(\text{sig})}_{\ell\ell^{\prime}}>0.4$ as the criteria for having an edge between vertices $\ell$ and $\ell^{\prime}$ in $G_{0}$. In the first stage [illustrated in Fig. 2(a)], $G_{0}$ is divided into non-overlapping cliques—subgraphs within which any two vertices are directly connected—using the following iterative process. In each step $k=1,2,\ldots$, we identify a clique of the largest possible size (i.e., the number of vertices it contains), denote this clique as $Q_{k}$, remove $Q_{k}$ from the graph $G_{k-1}$, and then denote the remaining graph by $G_{k}$. Repeating this step for each $k$ until $G_{k}$ is empty, we obtain a sequence $Q_{1},Q_{2},\ldots,Q_{m}$ of non-overlapping cliques in $G_{0}$, indexed in the order of non-increasing size. In the second stage, we agglomerate these cliques, as illustrated in Fig. 2(b). Initially, we set $R_{k}=Q_{k}$ for each $k$. Then, for each $k=2,3,\ldots,m$, we either move all the vertices in $R_{k}$ to the largest group among $R_{1},\cdots,R_{k-1}$ for which at least $80\%$ of all the possible edges between that group and $R_{k}$ actually exist, or we keep $R_{k}$ unchanged if no group satisfies this criterion. Among the resulting groups, we denote those groups whose size is at least three by $R_{1},R_{2},\ldots,R_{m^{\prime}}$, $m^{\prime}\leq m$. A key advantage of our method over applying community-detection algorithms MacMahon:2015 is that the edge density threshold above can be optimized for the accuracy of cascade size prediction.

We test the effectiveness of our general algorithm on the Texas power grid. As Figs. 3(a) and 3(b) show, the block-diagonal structure of $\widehat{C}^{\text{(sig)}}$ indicating high correlation within each group and low correlation between different groups becomes evident when the components are reindexed according to the identified groups. We note, however, that individual component vulnerabilities do not necessarily correlate with the co-susceptibility group structure [see Fig. 3(d), in comparison with Fig. 3(c)]. We find that the sizes of the groups of co-susceptible components vary significantly across the $24$ snapshots of the Texas power grid, as shown in Fig. 3(e). The degree of co-susceptibility, as measured by the total number of co-susceptible components, is generally lower under an off-peak condition than the on-peak counterpart [Fig. 3(e)]. This is consistent with the smaller deviation from the no-correlation assumption observed in Fig. 3(f), where this deviation is measured by the relative difference in the variance, $\bar{\sigma}^{2}$ (defined above). Since high correlation within a group of components implies a high probability that many of them fail simultaneously, the groups identified by our algorithm tend to have high values of $\gamma_{\mathcal{I}}$. Indeed, our calculation shows that Eq. (2) is satisfied for all the $171$ co-susceptible groups found in the $24$ snapshots.

Given the groups of components generated through our algorithm, the co-susceptibility model is defined as the set of binary random variables ${\cal F}_{\ell}$ (different from $F_{\ell}$) following the dichotomized correlated Gaussian distribution emrich1991method ; macke2009generating whose marginal probabilities (i.e., the probabilities that ${\cal F}_{\ell}=1$) equal the estimates of $p_{\ell}$ from the cascade-event data and whose correlation matrix ${\cal C}$ is given by

[TABLE]

We are thus approximating the correlation matrix $C$ by the block diagonal matrix ${\cal C}$, where the blocks correspond to the sets of co-susceptible components. In terms of the correlation network, this corresponds to using only those links within the same group of co-susceptible components for predicting the distribution of cascade sizes. Since individual groups are assumed to be uncorrelated, this can be interpreted as model dimensionality reduction, in which the dimension reduces from $n$ to the size of the largest group. We sample ${\cal F}_{\ell}$ using the code provided in Ref. macke2009generating . In this implementation, the computational time for sampling scales with the number of variables with an exponent of $3$, so factors of $2.0$ to $15.2$ in dimensionality reduction observed for the Texas power grid correspond to a reduction of computational time by factors of more than $8$ to more than $3{,}500$.

We now validate the co-susceptibility model for the Texas grid. We estimate the cumulative distribution function $S_{\cal N}(x)$ of cascade size, ${\cal N}\equiv\sum_{\ell}{\cal F}_{\ell}$, using $3{,}000$ samples generated from the model. As shown in Fig. 4(a), this function matches well with the cumulative distribution function $S_{N}(x)$ of cascades size $N$ computed directly from the cascade-event data. This is validated more quantitatively in the inset; the (binned) probability $p_{\cal N}(x)$ that $x\leq{\cal N}\leq x+\Delta x$ for the co-susceptibility model is plotted against the corresponding probability $p_{N}(x)$ for the cascade-event data, using a bin size of $\Delta x=N_{\max}/20$, where $N_{\max}$ denotes the maximum cascade size observed in the cascade-event data. The majority of the points lie within the 95% confidence interval for $p_{\cal N}(x)$, computed using the estimated $p_{N}(x)$. To validate the co-susceptibility model across all 24 snapshots, we use the Kolmogorov-Smirnov (KS) test massey1951kolmogorov . Specifically, for each snapshot we test the hypothesis that the samples of ${\cal N}$ and the corresponding samples of $N$ are from the same distribution. Figure 4(b) shows the measure of distance between two distributions, $\sup_{x}|S_{N}(x)-S_{\cal N}(x)|$, which underlies the KS test, as a function of the total amount of electrical load in the system. We find that the null hypothesis cannot be rejected at the 5% significance level for most of the cases we consider [$21/24=87.5$%, blue dots in Fig. 4(b)]; it can be rejected in only three cases (red triangles, above the threshold distance indicated by the dashed line), all corresponding to high stress (i.e., high load) conditions. We also see that more stressed systems are associated with larger distances between the distributions, and a higher likelihood of being able to reject the null hypothesis. We believe this is mainly due to higher-order correlations not captured by $p_{\ell}$ and $C$.

The identification of co-susceptibility as a key ingredient in determining cascade sizes leads to two new questions: (1) What gives rise to co-susceptibility? (2) How to identify the co-susceptible groups? While the first question opens an avenue for future research, the second question is addressed by the algorithm developed here (for which we provide a ready-to-use software source-code ). The co-susceptibility model is general and can be used for cascades of any type (of failures, information, or any other spreadable attribute) for which information is available on the correlation matrix and the individual “failure” probabilities. Such information can be empirical, as in the financial data studied in Ref. Plerou:2002 , or generated from first-principle models, as in the power-grid example used here. Our approach accounts for correlations (a strength shared by some other approaches, such as the one based on branching processes dobson2012 ), and does so from the fresh, network-based perspective of co-susceptibility. Finally, since co-susceptibility is often a nonlocal effect, our results suggest that we may need nonlocal strategies for reducing the risk of cascading failures, which bears implications for future research.