Emergent failures and cascades in power grids: a statistical physics perspective
Tommaso Nesti,
1
Alessandro Zocca,
2
and Bert Zwart
1
1
CWI, Amsterdam 1098 XG, Netherlands
2
California Institute of Technology, Pasadena, California 91125, USA
(Dated: April 23, 2018)
We model power grids transporting electricity generated by intermittent renewable sources as
complex networks, where line failures can emerge indirectly by noisy power input at the nodes. By
combining concepts from statistical physics and the physics of power flows, and taking weather
correlations into account, we rank line failures according to their likelihood and establish the most
likely way such failures occur and propagate. Our insights are mathematically rigorous in a small-noise
limit and are validated with data from the German transmission grid.
PACS numbers: 89.75.Hc,89.20.-a,88.80.-q
Understanding cascading failures in complex networks
is of great importance and has received a lot of atten-
tion in recent years [
1
–
17
]. Despite proposing different
mechanisms for their evolution, a common feature is that
cascades are triggered by some
external
event. This initial
attack is chosen either (i) deliberately, to target the most
vulnerable or crucial network component or (ii) uniformly
at random, to understand the average network reliability.
This distinction led to the insight that complex networks
are resilient to random attacks, but vulnerable to targeted
attacks [
7
,
18
,
19
]. However, both lead to the
direct
failure
of the attacked network component.
In this Letter, we focus on networks in which edge fail-
ures occur in a fundamentally different manner. Specifi-
cally, we consider networks where fluctuations of the node
inputs can trigger edge failures. The realization (which
we call configuration) of the noise at the nodes is not only
the cause of edge failures, but can also impact the way
they propagate in the network.
We present our results in the context of power grids that
transport electricity generated by solar and wind parks.
In power grids, line failures can arise when the network
is driven from a stable state to a critically loaded state
by external factors; intermittent power generation at the
nodes causes random fluctuations in the line power flows,
possibly triggering outages and cascading failures. Thus,
line failures can emerge
indirectly
due to the interplay
between noisy correlated (due to weather) power input at
the nodes, the network structure, and power flow physics.
This interplay is challenging to analyze, yet this problem
is urgent as the penetration of renewable energy sources
is increasing [20, 21].
We analyze this interplay using statistical physics and
large deviations theory. We consider a parsimonious static
stochastic model similar to [
22
], introduce a scaling param-
eter
ε
describing the magnitude of the noise and consider
the regime
ε
→
0. In the limit, we can identify the most
vulnerable lines and explicitly determine the most likely
configuration of power inputs leading to failures and sub-
sequent propagating failures. These results are validated
using real data for the German transmission network.
Previous works applying large-deviations techniques to
problems in complex networks dynamics, such as epidemic
extinction and biophysical networks, include [23, 24].
We model a transmission network by a connected graph
G
with
n
nodes representing the
buses
and
m
directed
edges modeling
transmission lines
. The nominal val-
ues of net power injections at the nodes are given by
μ
=
{
μ
i
}
i
=1
,...,n
. We model the stochastic fluctuation
of the power injections around
μ
, due to variability in
renewable generation, by means of the random vector
p
=
{
p
i
}
i
=1
,...,n
, which is assumed to follow a multivari-
ate Gaussian distribution with density
φ
(
x
) =
exp(
−
1
2
(
x
−
μ
)
T
(
ε
Σ
p
)
−
1
(
x
−
μ
))
(2
π
)
n
2
det(
ε
Σ
p
)
1
2
,
(1)
with
ε
Σ
p
∈
R
n
×
n
being the covariance matrix of
p
. In
our theoretical analysis, we assume that
Σ
p
is known and
let
ε
→
0.
The Gaussian assumption is debatable, both for solar
and wind. While consistent with atmospheric physics [
25
]
and recent wind park statistics [
26
,
27
], different models
are preferred for different timescales [
28
–
31
]. An extension
of our framework to the dynamic model in [
31
] looks
promising (using Freidlin-Wentzell theory as in [
32
]). For
a static non-Gaussian extension, see [33].
Assuming the vector
μ
has zero sum and using the DC
approximation [
20
], the line power flows
f
=
{
f
i
}
i
=1
,...,m
are given by
f
=
Vp
,
(2)
where
V
is an
m
×
n
matrix encoding the grid topol-
ogy and parameters (i.e., line susceptances). The DC
approximation is commonly used in transmission system
analysis [
34
–
37
]. More realistic nonlinear models based
on AC power flows [
38
] may be analyzed leveraging the
contraction principle [39].
The total net power injected in the network
∑
n
i
=1
p
i
is
non-zero as
p
is random. Automated affine response and
redispatch mechanisms take care of this issue in power
grids. Mathematically, this corresponds to a “distributed
arXiv:1709.10166v3 [physics.soc-ph] 20 Apr 2018
2
(a)
Nominal line flows
|
ν
`
|
at 11am.
(b)
True overload probabilities
log
10
P
(
|
f
`
|≥
1) at 11am.
(c)
Top 5% of most likely lines to fail
(red) at 11am, according to (3), and
nominal injections from renewable
sources.
FIG. 1
slack” in our model: the total power injection mismatch
is distributed uniformly among all nodes (the matrix
V
accounts for this; see [33]).
In view of Eqs.
(1)
-
(2)
, the line power flows
f
also fol-
low a multivariate Gaussian distribution with mean
ν
and covariance matrix
ε
Σ
f
. The vector
ν
=
V
μ
∈
R
m
describes the nominal line flows, while the covariance
matrix
ε
Σ
f
=
ε
V
Σ
p
V
T
describes the correlations be-
tween line flows fluctuations, taking into account both
the correlations of the power injections (encoded by
Σ
p
)
and correlations created by the network topology due to
power flow physics (Kirchhoff’s laws) via
V
.
A line
overloads
if the absolute amount of power flowing
in it exceeds a given
line threshold
. We assume that such
overloads immediately lead to the outage of the corre-
sponding line, to which we will henceforth refer simply
as
line failure
. The rationale behind this assumption is
that there are security relays on high voltage transmis-
sion lines performing an emergency shutdown as soon
as the current exceeds a dangerous level. Without such
mechanisms, lines may overheat, sag and eventually trip.
We can express the line flows in units of the line thresh-
old by incorporating the latter in the definition of
V
[
33
],
so that
f
is the vector of
normalized
line power flows and
the failure of line
`
corresponds to
|
f
`
| ≥
1. We let the
power grid operate on average safely by assuming that
max
`
=1
,...,m
|
ν
`
|
<
1, so that only large fluctuations of
line flows lead to failures.
We are most interested in scenarios where power grids
are highly stressed, meaning that the nominal power
injections
{
μ
i
}
i
=1
,...,n
are such that the corresponding
nominal line power flows
{
ν
`
}
`
=1
,...,m
are close to their
thresholds. Such a stress could be caused by very high
wind generation [40].
An illustrative scenario is reported in Fig. 1a, which
depicts a snapshot of nominal line flows on the SciGRID
German network [
41
]. SciGRID is a detailed model of the
actual German transmission network with
n
= 585 buses
and
m
= 852 lines that we use as main illustration. The
dataset includes load/generation time series, line limits,
grid topology and generation costs. In our case study,
we obtain
μ
by solving an Optimal Power Flow problem
(OPF [
42
]) based on realistic data for wind and solar
generation, and we estimate
ε
Σ
p
using ARMA models;
for details see the supplement [
33
], which also describes a
setting covering conventional controllable power plants.
We now turn to the analysis of emergent failures and
their propagation using large deviations theory [
43
]. We
begin by deriving the exponential decay of probabilities
of single line failure events
|
f
`
|≥
1 for
`
= 1
,...,m
. As
line power flows are Gaussian, we obtain, see Example
3.1 in [43], that
I
`
=
−
lim
ε
→
0
ε
log
P
ε
(
|
f
`
|≥
1) =
(1
−|
ν
`
|
)
2
2
σ
2
`
,
(3)
where
σ
2
`
= (
Σ
f
)
``
. We call
I
`
the
decay rate
of the failure
probability of line
`
. Thus, for small
ε
, we approximate
the probability of the emergent failure of line
`
as
P
(
|
f
`
|≥
1)
≈
exp(
−
I
`
/ε
) = exp
(
−
(1
−|
ν
`
|
)
2
2
εσ
2
`
)
,
(4)
and that of the first emergent failure as
P
(max
`
|
f
`
|≥
1)
≈
exp(
−
min
`
I
`
/ε
)
.
(5)
These approximations for failure probabilities may not
be sharp in general, even when
ε
is small, since all terms
that are decaying subexponentially in 1
/ε
are ignored.
Nevertheless, Eq.
(4)
is quite useful for ranking purposes,
allowing to explicitly identify the lines that are most
likely to fail. To verify this empirically, we note that
the expression in Eq.
(4)
only depends on the product
εσ
2
`
=
ε
(
V
Σ
p
V
T
)
``
, and thus, ultimately, only on the
product
ε
Σ
p
, which in our case study we estimate directly
from the SciGRID data, see [33].
3
Fig. S2 shows the heatmap for the exact line failure
probabilities
P
(
|
f
`
| ≥
1), for the same day and hour as
in Fig. 1a: it is clear that a larger
|
ν
`
|
does not neces-
sarily imply a higher chance of failure. Fig. 1c depicts
the 5% most likely lines to fail, ranked according to
I
`
.
The ranking based on the large deviations approximation
successfully recovers the most likely lines to fail, and, in
fact, yields the same ordering as the one based on exact
probabilities [
33
], thus providing an accurate indicator of
system vulnerabilities.
Fig. 1c also illustrates the nominal renewable generation
mix: the buses housing stochastic power injections have
different colors (blue/light blue for wind offshore/onshore,
yellow for solar) and sizes proportional to the absolute
values of the corresponding nominal injections. Many
vulnerable lines are located where the most renewable
energy production occurs. However, the interplay between
network topology, power flows physics and correlation in
power injections caused by weather fluctuations, results
in a spread-out arrangement of vulnerable lines, which is
hard to infer by looking at nominal values only.
We proceed with an analysis of how emergent failures
occur, using again large deviations theory. In particular,
we provide an explicit estimate of the most likely power
injection that caused a specific emergent failure. To this
end, we fix a line
`
and consider the conditional distribu-
tion of
p
, given
|
f
`
| ≥
1. The mean of this distribution
greatly simplifies as
ε
→
0 to
p
(
`
)
=
arg inf
p
∈
R
n
:
|
ˆ
e
T
`
Vp
|≥
1
1
2
(
p
−
μ
)
T
Σ
−
1
p
(
p
−
μ
)
.
(6)
If
ν
`
6
= 0, the solution is unique and reads
p
(
`
)
=
μ
+
(sign(
ν
`
)
−
ν
`
)
σ
2
`
Σ
p
V
T
ˆ
e
`
,
(7)
where
sign
(
a
) = 1 if
a
≥
0 and
−
1 otherwise, and
ˆ
e
`
∈
R
m
is the
`
-th unit vector. As
ε
→
0, the conditional variance
of
p
given
|
f
`
| ≥
1 decreases to 0 exponentially fast in
1
/ε
, yielding that the conditional distribution of
p
given
|
f
`
|≥
1 gets sharply concentrated around
p
(
`
)
[33].
We interpret
p
(
`
)
as the
most likely
power injection pro-
file, conditional on the failure of line
`
. The corresponding
line power flow profile
f
(
`
)
=
Vp
(
`
)
is
f
(
`
)
k
=
ν
k
+
(sign(
ν
`
)
−
ν
`
)
σ
2
`
Cov(
f
`
,f
k
)
,
∀
k
6
=
`.
(8)
As such, our framework provides more explicit informa-
tion than the approach in [
44
], which approximates the
most likely way events happen using the mode, without
leveraging large deviations. In our validation experiments,
we found that the error between
p
(
`
)
and
E
[
p
||
f
`
|≥
1]
is typically less than 1% of the nominal values [
33
]. A
numerical illustration is given in Fig. 2b.
A key finding is that an emergent line failure does not
occur due to large fluctuations only in neighboring nodes,
but as a cumulative effect of small unusual fluctuations in
the entire network “summed up” by power flow physics,
and correlations in renewable energy. Such an emergent
failure requires every line flow to be driven to an unusual
state
f
(
`
)
k
, which deviates from the nominal value
ν
k
by
an amount proportional to the covariance
Cov
(
f
`
,f
k
), in
view of Eq. (7).
We continue by investigating the propagation of failures,
combining our results describing the most likely power
injections configuration leading to the first failure, and
the power flow redistribution in the network afterwards.
To this end, we first differentiate between different types
of line failures, by assessing whether the most likely way
for failure of line
`
to occur is as (i) an
isolated failure
, if
|
f
(
`
)
k
|
<
1
for all line
k
6
=
`
, or (ii) a
joint failure
, if there
exists some other line
k
6
=
`
such that
|
f
(
`
)
k
|≥
1.
Any type of line failure(s) cause(s) a global redistri-
bution of the line power flows according to Kirchhoff’s
laws, which could trigger further outages and cascades. In
our setting, the power injections configuration
p
(
`
)
redis-
tributes across an altered network
̃
G
(
`
)
(a subgraph of the
original graph
G
) in which line
`
(and possible other lines,
in case of a joint failure) has been removed, increasing
stress on the remaining lines. The way this redistribution
happens on
̃
G
(
`
)
is governed by power flow physics and
we assume that it occurs instantaneously. Extending this
to dynamic models [
45
,
46
] is a natural future topic, as
transient oscillatory effects may severe the impact of line
failures.
The power flow redistribution amounts to compute a
new matrix
̃
V
linking the power injections and the new
power flows, which can be constructed analogously to
V
[
33
]. The most likely power flow configuration on
̃
G
(
`
)
after redistribution is
̃
f
(
`
)
=
̃
Vp
(
`
)
.
In the special case of an isolated failure (say of line
`
) it is enough to calculate the vector
φ
(
`
)
∈
R
m
−
1
of
(normalized) redistribution coefficients, known as
line
outage distribution factors
(LODF) [
47
]. The quantity
φ
(
`
)
j
takes values in [
−
1
,
1], and
|
φ
(
`
)
j
|
represents the percentage
of power flowing in line
`
that is redirected to line
j
after the failure of the former. The most likely power
flow configuration on
̃
G
(
`
)
after redistribution then equals
̃
f
(
`
)
=
{
f
(
`
)
k
}
k
6
=
`
+
f
(
`
)
`
φ
(
`
)
,
where
f
(
`
)
`
=
±
1 depending
on the way the power flow is most likely to exceed the
threshold 1. The power flow configuration
̃
f
(
`
)
can be
efficiently used to determine which lines subsequently fail,
by checking for which
k
we have
|
̃
f
(
`
)
k
|≥
1, see [33].
There is much evidence that failures propagate non-
locally in power grids [
48
–
52
]. To analyze this in our
framework we first consider a ring network with
μ
= 0
and
Σ
p
=
I
. In this network there are two paths along
which power can flow between any two nodes, using the
convention that a positive flow corresponds to a counter-
clockwise direction. If line
`
fails, the power originally
flowing on line
`
must now flow on the remaining path
4
(a)
After the emergent failure of line 27
(red) six additional lines (orange) fail,
4pm.
(b)
Most likely power injection
p
(
`
)
causing the isolated failure of line 720 (red),
and subsequent failures (orange). The bus sizes reflect how much
p
(
`
)
deviates
from
μ
at 11am (red positive deviations, blue negative). Left, with correlation in
noise); Right, without correlation in noise (setting to 0 all the off-diagonals of
Σ
p
).
FIG. 2
in the opposite direction. To make this rigorous we show
in [
33
] that
φ
(
`
)
k
=
−
1
for every
k
6
=
`
. As power flows
must sum to zero by Kirchhoff’s law, neighboring lines
tend to have positively correlated power flows, while flows
on distant lines exhibit negative correlations. Hence, the
power injections that make the power flows in line
`
exceed
the line threshold (say by becoming larger than 1) also
make the power flows in the antipodal half of the network
negative. These will go beyond the line threshold
−
1 after
the power flow redistributes, cf. Fig. 3.
1
1
/
7
-13
/
35
-19
/
35
-13
/
35
1
/
7
failed
1
/
7-1
-13
/
35-1
-19
/
35-1
-13
/
35-1
1
/
7-1
FIG. 3:
Left: most likely power injections
p
(
`
)
leading to the
failure of line
`
(orange), visualized using the color and size of the
nodes (red positive deviations, blue negative), together with power
flows
f
(
`
)
k
. Right: situation after the power flow redistribution with
three subsequent failures and the values
̃
f
(
`
)
k
=
f
(
`
)
k
−
1,
k
6
=
`
.
In the SciGRID example, Fig. 2a shows how the emer-
gent isolated failure of line
`
= 27 causes the failure of
six more lines
k
1
,...,k
6
, two of which are far way from
the original failure. For validation purposes, we found
numerically that
P
(
line
k
j
fails
∀
j
= 1
,...,
6
| |
f
27
| ≥
1)
≥
0
.
9987. Conversely, the failure of line 27 under the
nominal power injection profile leads to only two subse-
quent failures. The nontypical input caused other lines to
be more loaded than expected, and these lines get more
vulnerable as the cascades progresses, resulting in more
subsequent failures.
To validate this insight, we have looked at the first two
stages of emergent cascading failures for several IEEE
test networks, and compare them with those of classical
cascading failures, obtained using nominal power injection
values rather than the most likely ones and deterministic
removal of the initial failing line; see [
33
] for a precise de-
scription of the experiment. As before, emergent cascades
tend to lead to a higher number of subsequent failures in
each stage.
A non-diagonal noise matrix
Σ
p
exacerbates these ef-
fects. Experiments (see Fig. 2b) with our SciGRID case
study suggest that, if there is a correlation in noise, for ex-
ample due to fluctuations in weather patterns, the number
of subsequent failures can become higher. Furthermore,
it is easier for a failure to be triggered by many small
disturbances across the network, compared to the case
where these correlations are not taken into account. In
the latter case, we see a more local effect with relatively
larger disturbances.
In conclusion, we illustrated the potential of concepts
from statistical physics and large deviations theory to an-
alyze emergent failures and their propagation in complex
networks. Exogenous noise disturbances at the nodes,
potentially amplified by correlations, push a complex
network into a critical state in which edge failure may
emerge. Large deviations theory provides a tool to rank
such failures according to their likelihood and predicts
how such failures most likely occur and propagate. When
an emergent edge failure occurs, its impact on the network
can be more significant than a purely exogenous failure,
possibly resulting in cascades that propagate quicker than
in classical vulnerability analysis.
The accuracy of the small noise limit has been validated
in our case study, making the case for applying large
deviations techniques to more realistic models. In [
33
]
we propose a promising economic application of our
approach, showing how our framework can shed light
on the trade-off between network reliability and societal
costs.
5
Acknowledgements.
We thank the referees for many
useful comments, in particular for suggesting SciGRID.
NWO Vici 639.033.413 and NWO Rubicon 680.50.1529
grants provided financial support. AZ acknowledges the
support of Resnick Sustainability Institute at Caltech.
[1]
R. Albert, I. Albert, and G. Nakarado, Physical Review
E
69
, 025103 (2004).
[2]
R. Albert and A.-L. Barab ́asi, Reviews of Modern Physics
74
, 47 (2002).
[3]
R. Albert, H. Jeong, and A.-L. Barab ́asi, Nature
406
,
378 (2000).
[4]
P. Crucitti, V. Latora, and M. Marchiori, Physical Review
E
69
, 045104 (2004).
[5]
P. Crucitti, V. Latora, M. Marchiori, and A. Rapisarda,
Physica A: Statistical Mechanics and its Applications
320
, 622 (2003).
[6]
B. Mirzasoleiman, M. Babaei, M. Jalili, and M. Safari,
Physical Review E
84
, 046114 (2011).
[7]
A. Motter and Y.-C. Lai, Physical Review E
66
, 065102
(2002).
[8] A. Motter, Physical Review Letters
93
, 1 (2004).
[9]
D. Heide, M. Sch ̈afer, and M. Greiner, Physical Review
E
77
, 056103 (2008).
[10]
R. Kinney, P. Crucitti, R. Albert, and V. Latora, Euro-
pean Physical Journal B
46
, 101 (2005).
[11]
B. Sch ̈afer, C. Beck, K. Aihara, D. Witthaut,
and
M. Timme, Nature Energy
3
, 119 (2018).
[12]
S. Sun, Z. Liu, Z. Chen, and Z. Yuan, Physica A: Statis-
tical Mechanics and its Applications
373
, 851 (2007).
[13]
Y. Yang, T. Nishikawa, and A. E. Motter, Science
358
,
eaan3184 (2017).
[14]
D. Watts, Proceedings of the National Academy of Sci-
ences
99
, 5766 (2002).
[15]
D. Witthaut and M. Timme, Physical Review E
92
,
032809 (2015).
[16]
D. Witthaut, M. Rohden, X. Zhang, S. Hallerberg, and
M. Timme, Physical Review Letters
116
, 138701 (2016).
[17]
D. Witthaut and M. Timme, The European Physical
Journal B
86
, 377 (2013).
[18]
R. Cohen, K. Erez, D. Ben-Avraham, and S. Havlin,
Physical Review Letters
85
, 4626 (2000).
[19]
R. Cohen, K. Erez, D. Ben-Avraham, and S. Havlin,
Physical Review Letters
86
, 3682 (2001).
[20]
D. Bienstock,
Electrical Transmission System Cascades
and Vulnerability
(SIAM, Philadelphia, 2015) Chap. 4.
[21]
I. Dobson, B. Carreras, V. Lynch, and D. Newman,
Chaos: An Interdisciplinary Journal of Nonlinear Science
17
, 026103 (2007).
[22]
Z. Wang, A. Scaglione, and R. Thomas, in
2012 45th
Hawaii International Conference on System Sciences
(IEEE, 2012) pp. 2115–2124.
[23]
D. K. Wells, W. L. Kath, and A. E. Motter, Phys. Rev.
X
5
, 031036 (2015).
[24]
J. Hindes and I. B. Schwartz, Phys. Rev. Lett.
117
, 028302
(2016).
[25]
D. Bienstock, M. Chertkov, and S. Harnett, SIAM Review
56
, 461 (2014).
[26]
S. Kolumban, S. Kapodistria, and N. Nooraee, ArXiv
e-prints (2017), arXiv:1707.06497.
[27]
J. Berg, A. Natarajan, J. Mann, and E. Patton, Wind
Energy
19
, 1975 (2016).
[28]
A. S. Brouwer, M. van den Broek, A. Seebregts, and
A. Faaij, Renewable and Sustainable Energy Reviews
33
,
443 (2014).
[29]
D. Schlachtberger, S. Becker, S. Schramm,
and
M. Greiner, Energy Conversion and Management
125
,
336 (2016).
[30]
Y. Peings and G. Magnusdottir, Environmental Research
Letters
9
, 034018 (2014).
[31]
P. Milan, M. W ̈achter, and J. Peinke, Phys. Rev. Lett.
110
, 138701 (2013).
[32]
T. Nesti, J. Nair, and B. Zwart, ArXiv e-prints (2016),
arXiv:1606.02986.
[33]
“See Supplemental Material at ... for all the details about
the power flow model, the large deviations principle for
failures and the numerics for the German network, which
includes refs. [53-70],”.
[34]
K. Purchala, L. Meeus, D. Van Dommelen, and R. Bel-
mans, in
IEEE Power Engineering Society General Meet-
ing
(IEEE, 2005) pp. 2457–2462.
[35]
B. Stott, J. Jardim, and O. Alsac, IEEE Transactions on
Power Systems
24
, 1290 (2009).
[36]
L. Powell,
Power system load flow analysis
(McGraw Hill,
2004).
[37]
A. Wood, B. Wollenberg, and G. Sheble,
Power genera-
tion, operation, and control
, 3rd ed. (John Wiley & Sons,
2014).
[38]
D. Mehta, D. Molzahn, and K. Turitsyn, in
2016 Amer-
ican Control Conference (ACC)
(IEEE, 2016) pp. 1753–
1765.
[39]
A. Dembo and O. Zeitouni,
Large Deviations Techniques
and Applications
, Stochastic Modelling and Applied Prob-
ability, Vol. 38 (Springer, Berlin, Heidelberg, 2010).
[40]
T. Pesch, H.-J. Allelein, and J.-F. Hake, The European
Physical Journal Special Topics
223
, 2561 (2014).
[41]
C. Matke, W. Medjroubi, and D. Kleinhans, “SciGRID -
An Open Source Reference Model for the European Trans-
mission Network,”
http://http://scigrid.de
(2015).
[42]
M. Huneault and F. D. Galiana, IEEE Transactions on
Power Systems
6
, 762 (1991).
[43] H. Touchette, Physics Reports
478
, 1 (2009).
[44]
M. Chertkov, F. Pan, and M. Stepanov, IEEE Transac-
tions on Smart Grid
2
, 162 (2011).
[45]
I. Simonsen, L. Buzna, K. Peters, S. Bornholdt, and
D. Helbing, Phys. Rev. Lett.
100
, 218701 (2008).
[46]
B. Sch ̈afer, D. Witthaut, M. Timme, and V. Latora,
ArXiv e-prints (2017), arXiv:1707.08018.
[47]
J. Guo, Y. Fu, Z. Li, and M. Shahidehpour, IEEE Trans-
actions on Power Systems
24
, 1633 (2009).
[48]
D. Jung and S. Kettemann, Phys. Rev. E
94
, 012307
(2016).
[49] S. Kettemann, Phys. Rev. E
94
, 062311 (2016).
[50]
D. Labavi ́c, R. Suciu, H. Meyer-Ortmanns, and S. Ket-
temann, The European Physical Journal Special Topics
223
, 2517 (2014).
[51]
H. Ronellenfitsch, D. Manik, J. Horsch, T. Brown, and
D. Witthaut, IEEE Transactions on Power Systems
32
,
4060 (2017).
[52]
D. Manik, M. Rohden, H. Ronellenfitsch, X. Zhang,
S. Hallerberg, D. Witthaut, and M. Timme, Phys. Rev.
E
95
, 012319 (2017).
[53]
H. Cetinay, F. Kuipers, and P. Van Mieghem, IEEE
6
Systems Journal (2016).
[54]
M. Schaub, J. Lehmann, S. Yaliraki, and M. Barahona,
Network Science
2
, 66 (2014).
[55]
S. Soltan, D. Mazauric, and G. Zussman, IEEE Transac-
tions on Control of Network Systems
4
, 288 (2017).
[56]
R. Bapat, Ramanujan Math. Soc. Lect. Notes Ser
7
, 63
(2008).
[57]
T. Brown, J. H ̈orsch, and D. Schlachtberger, ArXiv e-
prints (2017), arXiv:1707.09913.
[58]
T.
Brown,
https://pypsa.org/examples/
scigrid-lopf-then-pf.html
(2017).
[59]
T. Brown,
https://pypsa.org/examples/add_load_
gen_trafos_to_scigrid.html
(2017).
[60]
OpenStreetMap
contributors,
“Planet
dump
re-
trieved
from
https://planet.osm.org
,”
https:
//www.openstreetmap.org
(2017).
[61]
M. Milligan, M. Schwartz, and Y. Wan, National Renew-
able Energy Laboratory, Golden, CO (2003).
[62]
J. Antonanzas, N. Osorio, R. Escobar, R. Urraca, F. M.
de Pison, and F. Antonanzas-Torres, Solar Energy
136
,
78 (2016).
[63]
R. Huang, T. Huang, R. Gadh, and N. Li, in
2012 IEEE
Third International Conference on Smart Grid Communi-
cations (SmartGridComm)
(2012) pp. 528–533.
[64]
J. Zhang, B.-M. Hodge, and A. Florita, in
ASME 2013
7th International Conference on Energy Sustainability
(2013) p. V001T16A003.
[65]
B. M. Hodge and M. Milligan, in
2011 IEEE Power and
Energy Society General Meeting
(2011).
[66]
F. Schweppe, M. Caramanis, R. Tabors, and R. Bohn,
Spot Pricing of Electricity
(Springer US, 1988).
[67]
P. Hines, I. Dobson, E. Cotilla-Sanchez, and M. Eppstein,
in
2013 46th Hawaii International Conference on System
Sciences
(IEEE, 2013) pp. 2141–2150.
[68]
P. Hines, I. Dobson, and P. Rezaei, IEEE Transactions
on Power Systems
32
, 1 (2016).
[69]
J. Qi, K. Sun, and S. Mei, IEEE Transactions on Power
Systems
30
, 804 (2015).
[70]
Y. Yang, T. Nishikawa, and A. Motter, Physical Review
Letters
118
, 048301 (2017).
1
Supplemental Material for:
Emergent failures and cascades in power grids: a statistical physics perspective
POWER GRID MODEL AND DC
APPROXIMATION
We model the power grid network as a connected
weighted graph
G
with
n
nodes, modeling buses, and
m
edges, representing the transmission lines. We make
use of the
DC approximation
, which is commonly used in
high-voltage transmission system analysis [34–37].
Choosing an arbitrary but fixed orientation of the trans-
mission lines, the network structure is described by the
edge-vertex incidence matrix
C
∈
R
m
×
n
defined as
C
`,i
=
1
if
`
= (
i,j
)
,
−
1 if
`
= (
j,i
)
,
0
otherwise
.
Denote by
β
`
=
β
i,j
=
β
j,i
>
0 the weight of edge
`
= (
i,j
), corresponding to the
susceptance
of that trans-
mission line. By convention, we set
β
i,j
=
β
j,i
= 0 if there
is no transmission line between
i
and
j
. Denote by
B
the
m
×
m
diagonal matrix defined as
B
=
diag
(
β
1
,...,β
m
).
The network topology and weights are simultaneously
encoded in the
weighted Laplacian matrix
of the graph
G
, defined as
L
=
C
>
BC
or entry-wise as
L
i,j
=
{
−
β
i,j
if
i
6
=
j,
∑
k
6
=
j
β
i,k
if
i
=
j.
All the rows of
L
sum up to zero and thus the matrix
L
is
singular. The eigenvalue zero has multiplicity one (thanks
to the assumption that the graph
G
is connected) and
the corresponding eigenvector is
1
. Denote by
v
2
,...,
v
n
the remaining eigenvectors of
L
, which are orthogonal to
1
and thus have all zero sum.
According to the
DC approximation
, the relation be-
tween any zero-sum vector of power injections
p
∈
R
n
and the phase angles
θ
∈
R
n
they induce in the network
nodes can be written in matrix form as
p
=
L
θ
,
Defining
L
+
∈
R
n
×
n
as the
Moore-Penrose
pseudo-inverse
of
L
, we can rewrite this as
θ
=
L
+
p
.
(S1)
This latter identity is particularly useful in our context,
since it holds for any vector of power injections
p
∈
R
n
,
even if it has no zero sum. Indeed, decomposing the
vector
p
using the basis of eigenvectors
1
,
v
2
,...,
v
n
of
L
+
one notices that the only component of
p
with non-
zero sum belongs to the null space of
L
+
(generated by
the eigenvector
1
).
This mathematical fact corresponds to the assumption
that the power grid has automatic redispatch/balancing
mechanisms, in which the total power injection mismatch
is distributed uniformly among all the nodes, thus ensur-
ing that the total net power injection is always zero.
Denote by
J
∈
R
n
×
n
the matrix with all entries equal
to one. Exploiting the eigenspace structure of
L
,
L
+
can
be calculated as
L
+
=
(
L
+
1
n
J
)
−
1
−
1
n
J
,
In the literature, instead of
L
+
it is commonly used
another matrix
̄
L
, calculated using the inverse of the
(
n
−
1)
×
(
n
−
1) sub-matrix obtained from
L
by means
of deleting the first row and first column. In our method
we are implicitly choosing an average value of zero as
a reference for the nodes voltage phase angles, while in
the classical one the first node is used as reference by
setting is phase angle equal to zero. We remark that
these two procedure are equivalent if one is interested in
the line power flows, as these latter depend only on the
phase angle differences. However, the matrix
̄
L
does not
account for the distributed slack, which needs to added
by post-multiplying by the matrix
S
=
I
−
1
n
J
∈
R
n
×
n
.
The real line power flows
ˆ
f
are related with the phase
angles
θ
via the linear relation
ˆ
f
=
BC
θ
. In view of
Eq.
(S1)
, the line power flow
ˆ
f
can be written as a linear
transformation of the power injections
p
, i.e.
ˆ
f
=
BCL
+
p
.
(S2)
It is convenient to look at the
normalized line power flow
vector
f
∈
R
m
, defined component-wise as
f
`
=
ˆ
f
`
/C
`
for
every
`
= 1
,...,m
, where
C
`
is the
line threshold
of line
`
,
which is assumed to be given. Line thresholds are in place
because a protracted current overload would heat up the
line, causing sag, loss of tensile strength and eventually
mechanical failure. If this happens, the failure may cause
a global redistribution of the line power flows which could
trigger cascading failures and blackouts.
The relation between line power flows and normalized
power flows can be rewritten as
f
=
W
ˆ
f
, where
W
is
the
m
×
m
diagonal matrix
W
=
diag
(
C
−
1
1
,...,C
−
1
m
). In
view of Eq.
(S2)
, the normalized power flows
f
can be
expressed in terms of the power injections
p
as
f
=
Vp
,
where
V
=
WBCL
+
∈
R
m
×
n
.