Vulnerability and co-susceptibility determine the size of network cascades
Yang Yang Takashi Nishikawa Adilson E. Motter





Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA


Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA


Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA


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.


Vulnerability and co-susceptibility determine the size of network cascades
    Adilson E. Motter
    Accepted XXX. Received YYY; in original form ZZZ
==========================================================================



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 <cit.> and air traffic delay propagation <cit.> to secondary species extinctions <cit.> and large social riots <cit.>.
Evident in numerous previous modeling efforts <cit.> 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. <ref>).
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 

    N = ∑_ℓ=1^n F_ℓ,

where F_ℓ is a binary random variable representing the failure status of component ℓ (i.e., F_ℓ = 1 if component ℓ fails during the cascade, and F_ℓ = 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 <cit.>.
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 <cit.>.

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 <cit.>).
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. <cit.> with the power re-balancing scheme from Ref. <cit.>),
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 <cit.>.
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 <cit.>.
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 <cit.>, accounting for several physical properties of failure propagation, sits relatively high in the hierarchy of existing power-grid cascade models <cit.>, which ranges from the most detailed engineering models to simplest graphical or stochastic models.
The model has also been validated against historical data <cit.>.

In general, mutual dependence among the variables F_ℓ may be necessary to explain the distribution of the cascade size N.
We define the vulnerability p_ℓ≡⟨ F_ℓ⟩ of component ℓ to be the probability that this component fails in a cascade event (including events with N=0).
If the random variables F_ℓ are uncorrelated (and thus have zero covariance), then N would follow   Poisson's binomial distribution <cit.>, with average μ̃=∑_ℓ p_ℓ and variance σ̃^2=∑_ℓ p_ℓ(1-p_ℓ).
However, the actual variance σ^2 of N observed in the cascade-event data is significantly larger than the corresponding value σ̃^2 under the no-correlation assumption for all 24 snapshots of the Texas power grid (with the relative difference, σ̅^2≡(σ^2 - σ̃^2)/σ̃^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_ℓ.
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_ℓ 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 <cit.>, ρ(λ) = K√((λ_+ - λ)(λ - λ_-))/(2 πλ T), where λ_±  = [1 ±√(T/K) ]^2. 
Since those eigenvalues falling between λ_- and λ_+ can be considered contributions from the noise, the sample correlation matrix C can be decomposed as C  = C^(ran)+C^(sig), where C^(ran) and C^(sig) are its random and significant parts, respectively, which can be determined from the eigenvalues and the associated eigenvectors <cit.>.
In the network visualization of Fig. <ref>(a), we show the correlation coefficients C_ℓℓ'^(sig) between components ℓ and ℓ' 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. <ref>(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 ℐ≡{ℓ_1,…,ℓ_m} is said to be co-susceptible if 

    γ_ℐ≡⟨ N_ℐ| N_ℐ≠0 ⟩ - n̅_ℐ/m - n̅_ℐ > γ_th,

where N_ℐ≡∑_j=1^m F_ℓ_j is the number of failures in a cascade event among the m components, ⟨ N_ℐ| N_ℐ≠0 ⟩ denotes the average number of failures among these components given that at least one of them fails, n̅_ℐ≡∑_j=1^m p_ℓ_j/[1-∏_k=1^m(1-p_ℓ_k)] ≥ 1 is the value ⟨ N_ℐ| N_ℐ≠0 ⟩ would take if F_ℓ_1,…,F_ℓ_m were independent.
Here we set the threshold in Eq. (<ref>) to be γ_th=σ_N_ℐ/(m - n̅_ℐ), where σ_N_ℐ^2 ≡∑_j=1^m p_ℓ_j(1-p_ℓ_j)/[1-∏_k=1^m(1-p_ℓ_k)] - n̅_ℐ^2 ∏_k=1^m(1-p_ℓ_k) is the variance of N_ℐ given N_ℐ≠0 for statistically independent F_ℓ_1,…,F_ℓ_m.
By definition, the co-susceptibility measure γ_ℐ equals zero if F_ℓ_1,…,F_ℓ_m are independent.
It satisfies -(n̅_ℐ-1)/(m-n̅_ℐ) ≤γ_ℐ≤ 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 γ_ℐ can be computed for a given set of components, identifying sets of co-susceptible components in a given network from Eq. (<ref>) becomes infeasible quickly as n increases due to combinatorial explosion.




Here we propose an efficient two-stage algorithm for identifying co-susceptible components <cit.>.
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 C^(sig)_ℓℓ'>0.4 as the criteria for having an edge between vertices ℓ and ℓ' in G_0.
In the first stage [illustrated in Fig. <ref>(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,…, 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, …, 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. <ref>(b).
Initially, we set R_k = Q_k for each k. 
Then, for each k=2,3,…,m, we either move all the vertices in R_k to the largest group among  R_1,⋯, 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,…,R_m', m'≤ m.
A key advantage of our method over applying community-detection algorithms <cit.> 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. <ref>(a) and <ref>(b) show, the block-diagonal structure of C^(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. <ref>(d), in comparison with Fig. <ref>(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. <ref>(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. <ref>(e)].
This is consistent with the smaller deviation from the no-correlation assumption observed in Fig. <ref>(f), where this deviation is measured by the relative difference in the variance, σ̅^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 γ_ℐ.
Indeed, our calculation shows that Eq. (<ref>) 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 F_ℓ (different from F_ℓ) following the dichotomized correlated Gaussian distribution <cit.> whose marginal probabilities (i.e., the probabilities that F_ℓ=1) equal the estimates of p_ℓ from the cascade-event data and whose correlation matrix C is given by

    C_ℓℓ' = C^(sig)_ℓℓ'   if ℓ,ℓ'∈ R_k for some k≤ m',
    
    0    otherwise.

We are thus approximating the correlation matrix C by the block diagonal matrix 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 F_ℓ using the code provided in Ref. <cit.>.
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_ N(x) of cascade size, N≡∑_ℓ F_ℓ, using 3,000 samples generated from the model. 
As shown in Fig. <ref>(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_ N(x) that x ≤ N≤ x+Δ 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 Δ 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_ 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 <cit.>.
Specifically, for each snapshot we test the hypothesis that the samples of N and the corresponding samples of N are from the same distribution. 
Figure <ref>(b) shows the measure of distance between two distributions, sup_x |S_N(x)-S_ 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. <ref>(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_ℓ 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 <cit.>). 
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. <cit.>, 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 <cit.>), 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.


This work was supported by ARPA-E Award No. DE-AR0000702.





hines2009large 
P. Hines, J. Apt, and S. Talukdar,
Large blackouts in North America: Historical trends and policy implications,
Energ. Policy 37, 5249 (2009). 

fleurquin2013
P. Fleurquin, J. J. Ramasco, and V. M. Eguiluz,
Systemic delay propagation in the US airport network,
Sci. Rep. 3, 1159 (2013).

eco:11
J. A. Estes, J. Terborgh, J. S. Brashares, M. E. Power, J. Berger, W. J. Bond, and D. A. Wardle,
Trophic downgrading of planet Earth,
Science 333, 301 (2011). 

sahasrabudhe2011rescuing
S. Sahasrabudhe and A. E. Motter,
Rescuing ecosystems from extinction cascades through compensatory perturbations,
Nature Commun. 2, 170 (2011).

watts2002
D. J. Watts,
A simple model of global cascades on random networks,
Proc. Natl. Acad. Sci. USA 99, 5766 (2002). 

brummitt2015
C. D. Brummitt, G. Barnett, and R. M. D'Souza,
Coupled catastrophes: sudden shifts cascade and hop among interdependent systems,
J. Roy. Soc. Interface 12, 20150712 (2015).

kinney_2005
R. Kinney, P. Crucitti, R. Albert, and V. Latora,  
Modeling cascading failures in the North American power grid,
Euro. Phys. J. B 46, 101 (2005). 

buldyrev_2010
S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin,
Catastrophic cascade of failures in interdependent networks,
Nature  464, 1025 (2010). 

goh_2003
K. I. Goh, D. S. Lee, B. Kahng, and D. Kim, 
Sandpile on scale-free networks,
Phys. Rev. Lett. 91, 148701 (2003).

motter_2004
A. E. Motter,  
Cascade control and defense in complex networks,
Phys. Rev. Lett. 93, 098701 (2004).

dobson2007complex
I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, 
Complex systems analysis of series of blackouts: Cascading failure, critical points, and self-organization,
Chaos 17, 026103 (2007).

bak_1987
P. Bak, C. Tang, and K. Wiesenfeld,
Self-organized criticality: An explanation of the 1/f noise,
Phys. Rev. Lett. 59, 381 (1987).

Witthaut:2015
D. Witthaut and M. Timme,
Nonlocal effects and countermeasures in cascading failures,
Phys. Rev. E 92, 032809 (2015).

dobson_2016
I. Dobson, B. A. Carreras, D. E. Newman, and J. M. Reynolds-Barredo,  
Obtaining statistics of cascading line outages spreading in an electric transmission network from standard utility data,
IEEE T. Power Syst. 99, 1 (2016). 

anghel2007stochastic
M. Anghel, K. A. Werley, and A. E. Motter,
Stochastic model for power grid dynamics,
Proc. 40th Int. Conf. Syst. Sci. HICSS'07, 
Big Island, HI, USA,
Vol. 1, 113 (2007).

FERC
The data for the snapshots are obtained from Federal Energy Regulatory Commission (FERC) Form 715.

Hines2011
P. Hines, E. Cotilla-Sanchez, and S. Blumsack,
Topological models and critical slowing down: Two approaches to power system blackout risk analysis,
Proc. 44th Int. Conf. Syst. Sci. HICSS'11, 
Kauai, HI, USA, 1 (2011).

source-code
Source code is available for download at <https://github.com/yangyangangela/determine_cascade_sizes>

opaModel
B. A. Carreras, D. E. Newman, I. Dobson, and N. S. Degala,
Validating OPA with WECC data,
Proc. 46th Int. Conf. Syst. Sci. HICSS'13, 
Maui, HI, USA, 2197 (2013).

henneaux2016
P. Henneaux, P. E. Labeau, J. C. Maun, and L. Haarla,
A two-level probabilistic risk assessment of cascading outages,
IEEE T. Power Syst. 31, 2393 (2016).

hines_2015
P. D. Hines, I. Dobson, and P. Rezaei,  
Cascading power outages propagate locally in an influence graph that is not the actual grid topology,
arXiv:1508.01775 (2016). 

dobson2012_vulnerability
B. A. Carreras, D. E. Newman, and I. Dobson,
Determining the vulnerabilities of the power transmission system,
Proc. 45th Int. Conf. Syst. Sci. HICSS'12,
Maui, HI, USA, 2044 (2012). 

Yang:2016
Y. Yang, T. Nishikawa, and A. E. Motter (to be published).

Wang:1993
Y. H. Wang, 
On the number of successes in independent trials,
Stat. Sinica. 3, 295 (1993). 

mehta2004random
M. L. Mehta,
Random Matrices, 3rd edn.
(Academic Press, 2004).

MacMahon:2015
M. MacMahon and D. Garlaschelli,
Community Detection for Correlation Matrices, 
Phys. Rev. X 5, 021006 (2015).

emrich1991method
L. J. Emrich and M. R. Piedmonte,
A method for generating high-dimensional multivariate binary variates,
Amer. Statist. 45, 302 (1991). 

macke2009generating
J. H. Macke, P. Berens, A. S. Ecker, A. S. Tolias, and M. Bethge,
Generating spike trains with specified correlation coefficients,
Neural Comp. 21, 397 (2009). 

massey1951kolmogorov
F. J. Massey, Jr.,
The Kolmogorov-Smirnov test for goodness of fit,
J. Amer. Statist. Assoc. 46, 68 (1951). 

Plerou:2002
V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Amaral, T. Guhr, and H. E. Stanley,
Random matrix approach to cross correlations in financial data,
Phys. Rev. E 65, 066126 (2002).

dobson2012
I. Dobson, 
Estimating the propagation and extent of cascading line outages from utility data with a branching process,
IEEE T. Power Syst. 27, 2146 (2012). 



