Vulnerability and co-susceptibility determine the size of network cascades

5.2 MB·

by  Yang Yang

Vulnerability and co-susceptibility determine the size of network cascades
Yang Yang,1 Takashi Nishikawa,1,2 and Adilson E. Motter1,2
1Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
2Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA
Inanetwork,alocaldisturbancecanpropagateandeventuallycauseasubstantialpartofthesys-
tem to fail, in cascade events that are easy to conceptualize but extraordinarily difficult to predict.
Here,wedevelopastatisticalframeworkthatcanpredictcascadesizedistributionsbyincorporating
twoingredientsonly: thevulnerabilityofindividualcomponentsandtheco-susceptibilityofgroups
of components (i.e., their tendency to fail together). Using cascades in power grids as a represen-
tative 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
7
studies,theseresultsprovideinsightsandanewmodelingframeworkforunderstandingcascadesin
1
financial systems, food webs, and complex networks in general.
0
2

n Thestabilityofcomplexnetworksislargelydetermined (a) (b)
1.0
a by their ability to operate close to equilibrium—a condi-
J
tionthatcanbecompromisedbyrelativelysmallpertur- 0.8
6 n
bations that can lead to large cascades of failures. Cas- o
2 cades are responsible for a range of network phenomena, 0.6elati
] from power blackouts [1] and air traffic delay propaga- 0.4orr
h C
tion [2] to secondary species extinctions [3, 4] and large
p 0.2
socialriots[5,6]. Evidentinnumerouspreviousmodeling
-
c efforts [3–12] is that dependence between components is 0.0
o
the building block of these self-amplifying processes and
s FIG. 1. Example of co-susceptibility of cascading failures
. canleadtocorrelationsamongeventualfailuresinacas-
s in a power grid. (a) Color-coded network of positive correla-
c cade. tionsbetweenfailuresoftransmissionlines(bluedots)inthe
i A central metric characterizing a cascade is its size. Texas network. In the gray background network, the nodes
s
y While the suitability of a size measure depends on the represent (all) transmission lines, while each link represents
h contextandpurpose,aconvenientmeasureisthenumber a direct physical connection through a common substation.
p of network components (nodes or links) participating in (b) Network in (a) after a correlation-based repositioning of
[ the nodes.
thecascade(e.g.,failedpowerlines,delayedairplanes,ex-

1 tinctspecies). Sincetherearemanyknownandunknown
v factors that can affect the details of cascade dynamics, ures can induce a sequence of failures in other com-
0
the main focus in the literature has been on characteriz- ponents. Here we assume that the initial failures and
9
7 ing the statistics of cascade sizes rather than the size of the propagation of failures can be modeled as stochastic
8 individual events. This leads to a fundamental question: and deterministic processes, respectively (although the
0 what determines the distribution of cascade sizes? framework also applies if the propagation or both are
1. In this Letter, we show that cascading failures (and stochastic). Thus, the cascade size N, defined here as
0 hence their size distributions) are often determined pri- thetotalnumberofcomponentsthatfailaftertheinitial
7 marily by two key properties associated with failures of failures, is a random variable that can be expressed as
1
the system components: the vulnerability, or the failure
: n
v probability of each component, and the co-susceptibility, (cid:88)
N = F , (1)
i orthetendencyofagroupofcomponentstofailtogether. (cid:96)
X
(cid:96)=1
The latter is intimately related to pairwise correlations
r
a betweenfailures,aswewillseebelow. Weprovideacon- where F is a binary random variable representing the
(cid:96)
crete algorithm for identifying groups of co-susceptible failure status of component (cid:96) (i.e., F = 1 if component
(cid:96)
componentsfor anygiven network. We demonstratethis (cid:96) fails during the cascade, and F =0 otherwise). While
(cid:96)
using the representative example of cascades of overload the n components may be connected by physical links,
failuresinpowergrids(Fig.1). Basedonourfindings,we a component may fail as the cascade propagates even if
develop the co-susceptibility model—a statistical model- none of its immediate neighbors have failed [6, 13, 14].
ing framework capable of accurately predicting the dis- Forexample,inthecaseofcascadingfailuresoftransmis-
tribution of cascade sizes, depending solely on the vul- sionlinesinapowergrid,thefailureofonelinecancause
nerability and co-susceptibility of component failures. a reconfiguration of power flows across the network that
We consider a system of n components subject to cas- leads to the overloading and subsequent failure of other
cading failures, in which a set of initial component fail- lines away from the previous failures [14, 15].
2
A concrete example network we analyze throughout dardize F by subtracting the average and dividing it by
(cid:96)
this Letter using the general setup above is the Texas thestandarddeviation. Accordingtorandommatrixthe-
powergrid,forwhichwehave24snapshots,representing ory, the probability density of eigenvalues of the correla-
on- and off-peak power demand in each season of three tion matrix computed from K samples of T independent
consecutiveyears. Eachsnapshotcomprisesthetopology random variables follow the Marchenko-Pastur distribu-
(cid:112)
of the transmission grid, the capacity threshold of each tion [25], ρ(λ) = K (λ −λ)(λ−λ )/(2πλT), where
+ −
line,thepowerdemandofeachloadnode,andthepower λ = (cid:2)1±(cid:112)T/K(cid:3)2. Since those eigenvalues falling be-
±
supply of each generator node (extracted from the data tween λ and λ can be considered contributions from
− +
reportedtoFERC[16]). Foreachsnapshotweuseaphys-
the noise, the sample correlation matrix C(cid:98) can be de-
icalcascademodeltogenerateK =5,000cascadeevents. composedasC(cid:98) =C(cid:98)(ran)+C(cid:98)(sig), whereC(cid:98)(ran) andC(cid:98)(sig)
Inthismodel(whichisavariantofthatinRef.[15]with
are its random and significant parts, respectively, which
the power re-balancing scheme from Ref. [17]), an ini-
can be determined from the eigenvalues and the associ-
tialperturbationtothesystem(underagivencondition)
ated eigenvectors [26]. In the network visualization of
is modeled by the removal of a set of randomly selected Fig. 1(a), we show the correlation coefficients C(cid:98)(sig) be-
lines. Acascadefollowingtheinitialfailuresisthenmod- (cid:96)(cid:96)(cid:48)
tween components (cid:96) and (cid:96)(cid:48) estimated from the cascade-
eled as an iterative process. In each step, power flow
event data for the Texas grid under the 2011 summer
is redistributed according to Kirchhoff’s law and might
on-peak condition. Note that we compute correlation
therefore cause some lines to be overloaded and removed
only between those components that fail more than once
(i.e., to fail) due to overheating. The temperature of the
in the cascade events. As this example illustrates, we
transmission lines is described by a continuous evolution
observe no apparent structure in a typical network vi-
model and the overheating threshold for line removal is
sualization of these correlations. However, as shown in
determined by the capacity of the line [15]. When a fail-
Fig. 1(b), after repositioning the nodes based on correla-
ure causes part of the grid to be disconnected, we re-
tion strength, we can identify clusters of positively and
balance power supply and demand under the constraints
strongly correlated components—those that tend to fail
of limited generator capacity [17]. A cascade stops when
together in a cascade.
no more overloading occurs, and we define the size N
To more precisely capture this tendency of simulta-
of the cascade as the total number of removed lines (ex-
neous failures, we define a notion of co-susceptibility: a
cluding the initial failures). This model [18], accounting
given subset of m components I ≡ {(cid:96) ,...,(cid:96) } is said
forseveralphysicalpropertiesoffailurepropagation,sits 1 m
to be co-susceptible if
relativelyhighinthehierarchyofexistingpower-gridcas-
cademodels[19–22],whichrangesfromthemostdetailed (cid:104)N |N (cid:54)=0(cid:105)−n¯
γ ≡ I I I >γ , (2)
engineering models to simplest graphical or stochastic I m−n¯ th
I
models. The model has also been validated against his-
torical data [23]. where N ≡(cid:80)m F is the number of failures in a cas-
I j=1 (cid:96)j
Ingeneral, mutualdependenceamongthevariablesF cade event among the m components, (cid:104)N |N (cid:54)= 0(cid:105) de-
(cid:96) I I
may be necessary to explain the distribution of the cas- notes the average number of failures among these com-
cade size N. We define the vulnerability p ≡ (cid:104)F (cid:105) of ponents given that at least one of them fails, n¯ ≡
(cid:96) (cid:96) I
component (cid:96) to be the probability that this component (cid:80)m p /(cid:2)1−(cid:81)m (1−p )(cid:3)≥1isthevalue(cid:104)N |N (cid:54)=
j=1 (cid:96)j k=1 (cid:96)k I I
failsinacascadeevent(includingeventswithN =0). If 0(cid:105) would take if F ,...,F were independent. Here we
(cid:96)1 (cid:96)m
therandomvariables F areuncorrelated(andthushave set the threshold in Eq. (2) to be γ = σ /(m−n¯ ),
zerocovariance),thenN(cid:96) wouldfollowPoisson’sbinomial where σ2 ≡ (cid:80)m p (1−p )/(cid:2)1t−h(cid:81)mNI(1−p )(cid:3)I−
σd˜i2st=rib(cid:80)ut(cid:96)ipo(cid:96)n(1[2−4p],(cid:96))w.iHthowaevveerra,gtehµe˜a=ctu(cid:80)al(cid:96)vpa(cid:96)riaanndcevσa2rioafnNce ns¯t2Iat(cid:81)ismkt=ic1aN(l1Ily−inpd(cid:96)ke)pj=eisn1dthe(cid:96)njetvFar(cid:96)1ia,n.(cid:96)j.c.e,oFf(cid:96)mN.IBgykiv=de1nefiNniIti(cid:54)=o(cid:96)nk0,tfhoer
observed in the cascade-event data is significantly larger co-susceptibility measure γ equals zero if F ,...,F
I (cid:96)1 (cid:96)m
thanthecorrespondingvalueσ˜2 undertheno-correlation are independent. It satisfies −(n¯ −1)/(m−n¯ )≤γ ≤
I I I
assumption for all 24 snapshots of the Texas power grid 1, where the (negative) lower bound is achieved if multi-
(with the relative difference, σ¯2 ≡(σ2−σ˜2)/σ˜2, ranging plefailuresneveroccurandtheupperboundisachievedif
fromaround0.18tonearly39). Thus,themutualdepen- allmcomponentsfailwheneveroneofthemfails. Thus,a
dencemustcontributetodeterminingthedistributionof setofco-susceptiblecomponentsarecharacterizedbysig-
N in these examples. nificantly larger number of simultaneous failures among
Part of this dependence is captured by the correlation these components, relative to the expected number for
matrix C, whose elements are the pairwise Pearson cor- statistically independent failures. While γ can be com-
I
relationcoefficientsamongthefailurestatusvariablesF . putedforagivensetofcomponents,identifyingsetsofco-
(cid:96)
When the correlation matrix is estimated from cascade- susceptible components in a given network from Eq. (2)
eventdata,ithasnoiseduetofinitesamplesize,whichwe becomes infeasible quickly as n increases due to combi-
filter out using the following procedure. First, we stan- natorial explosion.
3
(a) Partition (b) Agglomeration Line index Line index
1 95 1801 95 180
1
Q1 step 1: R1 R2 R3 R4 R5 R6 (a) 1 (b) 1.0
Q5 … …… ex 50 00..0068..68on
Q3 Q4 step 4: R1 R2 R3 R5 R6 Line ind95 100 00..0024..24Correlati
step 5: R1 R2 R3 R6 150 00.0
Q6 Q2 step 6: R1 R2 R3 R6 01.8002 1 50 100 150 -−00.2.2
(c) (d)
0.01 µ
FIG. 2. Two-stage algorithm for identifying sets of co-
susceptible components. (a) First stage, in which we par- 0
titionthereducedunweightedcorrelationgraphintomultiple 2010 Spring (e) (f)
2010 Summer
cliques {Q }, and index the cliques in descending order of
k 2010 Fall
size. (b) Second stage, in which we find the co-susceptible 2011 Winter
groups {R } by recursively agglomerating cliques that have 2011 Spring
k 2011 Summer
enough connections. 2011 Fall
2012 Winter ndex
22001122 SSpumrinmger LineI
Here we propose an efficient two-stage algorithm for 2012 Fall
2013 Winter
identifying co-susceptible components [18]. The algo-
0 50 100 150 200 -2 0 2 4
LineIndex
rithm is based on partitioning and agglomerating the Number of transmission lines
vertices of the auxiliary graph G in which vertices rep-
0 FIG. 3. Correlations between transmission line failures in
resent the components that fail more than once in the
the Texas power grid. (a–d) Estimated correlation matrix
cascade-event data, and (unweighted) edges represent C(cid:98)(sig) [(a),(b)] and vulnerabilities p(cid:96) [(c),(d)] under the 2011
the dichotomized correlation between these components. summer on-peak condition. In (a) and (c), the transmission
HereweuseC(cid:98)(cid:96)((cid:96)s(cid:48)ig) >0.4asthecriteriaforhavinganedge linesareindexedsothatp(cid:96) isdecreasingin(cid:96). In(b)and(d),
between vertices (cid:96) and (cid:96)(cid:48) in G . In the first stage [illus- theyareindexedsothatthesetsofco-susceptiblelinesappear
0
as diagonal blocks. (e) Sizes of the groups of co-susceptible
trated in Fig. 2(a)], G is divided into non-overlapping
0
transmission lines, indicated by the lengths of the individual
cliques—subgraphswithinwhichanytwoverticesaredi-
segmentsoftheredportionofthetopbar(on-peaksnapshot)
rectly connected—using the following iterative process.
andoftheblueportionofthebottombar(off-peaksnapshot)
In each step k = 1,2,..., we identify a clique of the for each season. The gray portion of each bar corresponds
largest possible size (i.e., the number of vertices it con- to groups of fewer than three lines among groups R . (f)
k
tains), denote this clique as Q , remove Q from the Relative difference σ¯2 between the variance of the cascade
k k
graph G , and then denote the remaining graph by sizeanditscounterpartundertheno-correlationassumption.
k−1
G . Repeating this step for each k until G is empty,
k k
we obtain a sequence Q ,Q ,...,Q of non-overlapping
1 2 m
cliquesinG ,indexedintheorderofnon-increasingsize.
0
In the second stage, we agglomerate these cliques, as il- however, that individual component vulnerabilities do
lustrated in Fig. 2(b). Initially, we set Rk =Qk for each not necessarily correlate with the co-susceptibility group
k. Then, for each k =2,3,...,m, we either move all the structure [see Fig. 3(d), in comparison with Fig. 3(c)].
vertices in Rk to the largest group among R1,··· ,Rk−1 Wefindthatthesizesofthegroupsofco-susceptiblecom-
for which at least 80% of all the possible edges between ponents vary significantly across the 24 snapshots of the
that group and Rk actually exist, or we keep Rk un- Texas power grid, as shown in Fig. 3(e). The degree of
changed if no group satisfies this criterion. Among the co-susceptibility, as measured by the total number of co-
resultinggroups, wedenotethosegroupswhosesizeisat susceptible components, is generally lower under an off-
leastthreebyR1,R2,...,Rm(cid:48),m(cid:48) ≤m. Akeyadvantage peak condition than the on-peak counterpart [Fig. 3(e)].
of our method over applying community-detection algo- This is consistent with the smaller deviation from the
rithms [26] is that the edge density threshold above can no-correlation assumption observed in Fig. 3(f), where
be optimized for the accuracy of cascade size prediction. this deviation is measured by the relative difference in
We test the effectiveness of our general algorithm on the variance, σ¯2 (defined above). Since high correlation
the Texas power grid. As Figs. 3(a) and 3(b) show, the within a group of components implies a high probability
block-diagonalstructureofC(cid:98)(sig) indicatinghighcorrela- that many of them fail simultaneously, the groups iden-
tion within each group and low correlation between dif- tified by our algorithm tend to have high values of γ .
I
ferent groups becomes evident when the components are Indeed, ourcalculation showsthat Eq.(2)is satisfiedfor
reindexed according to the identified groups. We note, all the 171 co-susceptible groups found in the 24 snap-
4
1.0 (a) 0.06 (b) more quantitatively in the inset; the (binned) probabil-
10−1 ity pN(x) that x≤N ≤x+∆x for the co-susceptibility
10−2 0.04 model is plotted against the corresponding probability
0.9 10−3 pN(x) for the cascade-event data, using a bin size of
Off-peak ∆x = N /20, where N denotes the maximum cas-
10−4 max max
0.02
cadesizeobservedinthecascade-eventdata. Themajor-
10−4 10−3 10−2 10−1
ityofthepointsliewithinthe95%confidenceintervalfor
On-peak
0.8 0.00 pN(x),computedusingtheestimatedpN(x). Tovalidate
0 10 20 30 40 2 4 6 8 theco-susceptibilitymodelacrossall24snapshots,weuse
Cascade size, Load (MW) x104
the Kolmogorov-Smirnov (KS) test [29]. Specifically, for
FIG.4. Validationoftheco-susceptibilitymodel. (a)Cumu- eachsnapshotwetestthehypothesisthatthesamplesof
lative distributions SN(x) and SN(x) of cascade sizes N and N andthecorrespondingsamplesofN arefromthesame
N fromthecascade-eventdataandfromtheco-susceptibility
distribution. Figure 4(b) shows the measure of distance
model, respectively, under the 2011 summer on-peak condi-
between two distributions, sup |S (x)−S (x)|, which
tion. Inset: Binned probabilities from the two distributions x N N
underlies the KS test, as a function of the total amount
plotted against each other. The shaded area indicates the
95% confidence interval. (b) Distance between two distribu- ofelectricalloadinthesystem. Wefindthatthenullhy-
tions as a function of the total load in megawatts (MW) for pothesis cannot be rejected at the 5% significance level
all 24 snapshots. 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
shots. the dashed line), all corresponding to high stress (i.e.,
Given the groups of components generated through high load) conditions. We also see that more stressed
our algorithm, the co-susceptibility model is defined as systems are associated with larger distances between the
the set of binary random variables F(cid:96) (different from distributions, and a higher likelihood of being able to re-
F(cid:96)) following the dichotomized correlated Gaussian dis- jectthenullhypothesis. Webelievethisismainlydueto
tribution [27, 28] whose marginal probabilities (i.e., the higher-order correlations not captured by p and C.
(cid:96)
probabilities that F =1) equal the estimates of p from
(cid:96) (cid:96) The identification of co-susceptibility as a key ingredi-
thecascade-eventdataandwhosecorrelationmatrixC is
ent in determining cascade sizes leads to two new ques-
given by
tions: (1) What gives rise to co-susceptibility? (2) How
(cid:40) to identify the co-susceptible groups? While the first
C = C(cid:98)(cid:96)((cid:96)s(cid:48)ig) if (cid:96),(cid:96)(cid:48) ∈Rk for some k ≤m(cid:48), (3) question opens an avenue for future research, the second
(cid:96)(cid:96)(cid:48)
0 otherwise. question is addressed by the algorithm developed here
(for which we provide a ready-to-use software [18]). The
We are thus approximating the correlation matrix C by co-susceptibilitymodelisgeneralandcanbeusedforcas-
theblockdiagonalmatrixC,wheretheblockscorrespond cades of any type (of failures, information, or any other
to the sets of co-susceptible components. In terms of the spreadable attribute) for which information is available
correlationnetwork, thiscorrespondstousingonlythose on the correlation matrix and the individual “failure”
linkswithinthesamegroupofco-susceptiblecomponents probabilities. Such information can be empirical, as in
for predicting the distribution of cascade sizes. Since the financial data studied in Ref. [30], or generated from
individual groups are assumed to be uncorrelated, this first-principle models, as in the power-grid example used
can be interpreted as model dimensionality reduction, in here. Our approach accounts for correlations (a strength
which the dimension reduces from n to the size of the shared by some other approaches, such as the one based
largest group. We sample F(cid:96) using the code provided in on branching processes [31]), and does so from the fresh,
Ref.[28]. Inthisimplementation,thecomputationaltime network-based perspective of co-susceptibility. Finally,
for sampling scales with the number of variables with an since co-susceptibility is often a nonlocal effect, our re-
exponent of 3, so factors of 2.0 to 15.2 in dimensionality sults suggest that we may need nonlocal strategies for
reduction observed for the Texas power grid correspond reducing the risk of cascading failures, which bears im-
to a reduction of computational time by factors of more plications for future research.
than 8 to more than 3,500.
This work was supported by ARPA-E Award No. de-
We now validate the co-susceptibility model for the
ar0000702.
Texas grid. We estimate the cumulative distribution
(cid:80)
function S (x) of cascade size, N ≡ F , using 3,000
N (cid:96) (cid:96)
samplesgeneratedfromthemodel. AsshowninFig.4(a),
this function matches well with the cumulative distri-
bution function SN(x) of cascades size N computed di- [1] P. Hines, J. Apt, and S. Talukdar, Large blackouts in
rectly from the cascade-event data. This is validated North America: Historical trends and policy implica-
5
tions, Energ. Policy 37, 5249 (2009). Energy Regulatory Commission (FERC) Form 715.
[2] P.Fleurquin,J.J.Ramasco,andV.M.Eguiluz,Systemic [17] P.Hines, E.Cotilla-Sanchez, andS.Blumsack, Topolog-
delay propagation in the US airport network, Sci. Rep. ical models and critical slowing down: Two approaches
3, 1159 (2013). to power system blackout risk analysis, Proc. 44th Int.
[3] J. A. Estes, J. Terborgh, J. S. Brashares, M. E. Power, Conf. Syst. Sci. HICSS’11, Kauai, HI, USA, 1 (2011).
J.Berger,W.J.Bond,andD.A.Wardle,Trophicdown- [18] Source code is available for download at https:
grading of planet Earth, Science 333, 301 (2011). //github.com/yangyangangela/determine_cascade_
[4] S.SahasrabudheandA.E.Motter,Rescuingecosystems sizes
from extinction cascades through compensatory pertur- [19] B.A.Carreras,D.E.Newman,I.Dobson,andN.S.De-
bations, Nature Commun. 2, 170 (2011). gala, Validating OPA with WECC data, Proc. 46th Int.
[5] D. J. Watts, A simple model of global cascades on ran- Conf. Syst. Sci. HICSS’13, Maui, HI, USA, 2197 (2013).
dom networks, Proc. Natl. Acad. Sci. USA 99, 5766 [20] P.Henneaux,P.E.Labeau,J.C.Maun,andL.Haarla,A
(2002). two-level probabilistic risk assessment of cascading out-
[6] C. D. Brummitt, G. Barnett, and R. M. D’Souza, ages, IEEE T. Power Syst. 31, 2393 (2016).
Coupled catastrophes: sudden shifts cascade and hop [21] P. D. Hines, I. Dobson, and P. Rezaei, Cascading power
amonginterdependentsystems,J.Roy.Soc.Interface12, outages propagate locally in an influence graph that is
20150712 (2015). not the actual grid topology, arXiv:1508.01775 (2016).
[7] R.Kinney,P.Crucitti,R.Albert,andV.Latora,Model- [22] B. A. Carreras, D. E. Newman, and I. Dobson, Deter-
ingcascadingfailuresintheNorthAmericanpowergrid, miningthevulnerabilitiesofthepowertransmissionsys-
Euro. Phys. J. B 46, 101 (2005). tem, Proc. 45th Int. Conf. Syst. Sci. HICSS’12, Maui,
[8] S.V.Buldyrev,R.Parshani,G.Paul,H.E.Stanley,and HI, USA, 2044 (2012).
S.Havlin,Catastrophiccascadeoffailuresininterdepen- [23] Y. Yang, T. Nishikawa, and A. E. Motter (to be pub-
dent networks, Nature 464, 1025 (2010). lished).
[9] K.I.Goh,D.S.Lee,B.Kahng,andD.Kim,Sandpileon [24] Y.H.Wang,Onthenumberofsuccessesinindependent
scale-free networks, Phys. Rev. Lett. 91, 148701 (2003). trials, Stat. Sinica. 3, 295 (1993).
[10] A. E. Motter, Cascade control and defense in complex [25] M. L. Mehta, Random Matrices, 3rd edn. (Academic
networks, Phys. Rev. Lett. 93, 098701 (2004). Press, 2004).
[11] I. Dobson, B. A. Carreras, V. E. Lynch, and D.E. New- [26] M. MacMahon and D. Garlaschelli, Community Detec-
man, Complex systems analysis of series of blackouts: tion for Correlation Matrices, Phys. Rev. X 5, 021006
Cascading failure, critical points, and self-organization, (2015).
Chaos 17, 026103 (2007). [27] L. J. Emrich and M. R. Piedmonte, A method for gen-
[12] P.Bak,C.Tang,andK.Wiesenfeld,Self-organizedcriti- erating high-dimensional multivariate binary variates,
cality: Anexplanationofthe1/f noise,Phys.Rev.Lett. Amer. Statist. 45, 302 (1991).
59, 381 (1987). [28] J. H. Macke, P. Berens, A. S. Ecker, A. S. Tolias, and
[13] D. Witthaut and M. Timme, Nonlocal effects and coun- M. Bethge, Generating spike trains with specified corre-
termeasures in cascading failures, Phys. Rev. E 92, lation coefficients, Neural Comp. 21, 397 (2009).
032809 (2015). [29] F.J.Massey,Jr.,TheKolmogorov-Smirnovtestforgood-
[14] I. Dobson, B. A. Carreras, D. E. Newman, and ness of fit, J. Amer. Statist. Assoc. 46, 68 (1951).
J. M. Reynolds-Barredo, Obtaining statistics of cascad- [30] V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Ama-
inglineoutagesspreadinginanelectrictransmissionnet- ral, T. Guhr, and H. E. Stanley, Random matrix ap-
workfromstandardutilitydata,IEEET.PowerSyst.99, proachtocrosscorrelationsinfinancialdata,Phys.Rev.
1 (2016). E 65, 066126 (2002).
[15] M. Anghel, K. A. Werley, and A. E. Motter, Stochastic [31] I.Dobson,Estimatingthepropagationandextentofcas-
model for power grid dynamics, Proc. 40th Int. Conf. cading line outages from utility data with a branching
Syst. Sci. HICSS’07, Big Island, HI, USA, Vol. 1, 113 process, IEEE T. Power Syst. 27, 2146 (2012).
(2007).
[16] The data for the snapshots are obtained from Federal