Supplemental Material for:
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)
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 [1–4].
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
.
(1)
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.
(1)
, the line power flow
ˆ
f
can be written as a linear
transformation of the power injections
p
, i.e.
ˆ
f
=
BCL
+
p
.
(2)
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
2
the
m
×
m
diagonal matrix
W
=
diag
(
C
−
1
1
,...,C
−
1
m
).
In view of Eq.
(2)
, the normalized power flows
f
can be
expressed in terms of the power injections
p
as
f
=
Vp
,
where
V
=
WBCL
+
∈
R
m
×
n
.
Stochastic and deterministic injections
We now briefly outline how the model presented above
can be extended to a setting where only a subset of nodes
houses stochastic power injections (modeling wind and
solar parks), while the other nodes house deterministic in-
jections (corresponding to conventional controllable power
plants).
First, we introduce the following notation: if
z
is a
n
-dimensional multivariate Gaussian random vector with
mean
λ
and covariance matrix
Λ
, it will be denoted by
z
∼N
n
(
λ
,
Λ
).
Define the following:
n
s
number of stochastic buses
,
n
d
number of deterministic buses
,
I
s
⊆{
1
,...,n
}
indices of stochastic buses
,
I
d
⊆{
1
,...,n
}
indices of deterministic buses
,
p
s
= (
p
i
)
i
∈I
s
∈
R
n
s
stochastic power injection
,
p
d
= (
p
i
)
i
∈I
d
∈
R
n
d
deterministic power injection
,
V
s
∈
R
m
×
n
s
matrix consisting of the
columns of
V
indexed by
I
s
,
V
d
∈
R
m
×
n
d
matrix consisting of the
columns of
V
indexed by
I
d
,
f
s
=
V
s
p
s
∈
R
m
stochastic component of
f,
f
d
=
V
d
p
d
∈
R
m
deterministic component of
f.
If a bus hosts both stochastic and deterministic genera-
tors, it is considered a stochastic bus. Stochastic power
injections are modelled by mean of a
n
s
-dimensional mul-
tivariate Gaussian random vector with mean
μ
s
∈
R
n
s
and covariance matrix
Σ
p
∈
R
n
s
×
n
s
, which we denote by
p
s
∼N
n
s
(
μ
s
,ε
Σ
p
)
,
With the previous notation, the normalized power flows
can be decomposed as
f
=
f
s
+
f
d
=
V
s
p
s
+
f
d
, where
f
s
∼N
m
(
ν
s
,ε
Σ
f
)
,
ν
s
=
V
s
μ
s
,
Σ
f
=
V
s
Σ
p
V
>
s
.
(3)
The nominal power flows values are thus equal to
ν
=
ν
s
+
f
d
. The decay rate for an overload in line
`
, analogously
to formula (6) in the Main Body of the paper, is given by
I
`
=
inf
p
s
∈
R
n
s
:
|
ˆ
e
>
`
(
V
s
p
s
+
f
d
)
|≥
1
1
2
(
p
s
−
μ
s
)
>
Σ
−
1
p
(
p
s
−
μ
s
)
.
Provided that
ν
`
6
= 0, the solution is unique and reads
p
(
`
)
s
=
(sign(
ν
`
)
−
ν
`
)
σ
2
`
Σ
p
V
>
s
ˆ
e
`
+
μ
s
∈
R
n
s
,
(4)
where
σ
2
`
= (
Σ
f
)
`,`
. The corresponding most likely real-
ization for power flows reads
f
(
`
)
=
V
s
p
(
`
)
s
+
f
d
=
(sign(
ν
`
)
−
ν
`
)
σ
2
`
V
s
Σ
p
V
>
s
ˆ
e
`
+
ν
s
+
f
d
∈
R
m
.
(5)
In the next section we prove these claims for the particular
case of
n
s
=
n
.
LARGE DEVIATIONS PRINCIPLES FOR
FAILURE EVENTS
Gaussian case
In this section we provide proofs for Eqs. (3)-(7) in
the Main Body. For the sake of clarity we present here
only the proofs for the case
n
=
n
s
, and we remark that
Eqs.
(4)
-
(5)
in the Supplemental Material can be proved
along similar lines. In the following, we write
p
ε
and
f
ε
to stress the dependence of the power injections and of
the line power flows on the noise parameter
ε
.
Proposition 1.
Assume that
max
j
=1
,...,m
|
ν
j
|
<
1
. Then,
for every
`
= 1
,...,m
, the sequence of line power flows
(
f
ε
)
ε>
0
satisfies the large deviations principle
lim
ε
→
0
ε
log
P
(
|
(
f
ε
)
`
|≥
1) =
−
(1
−|
ν
`
|
)
2
2
σ
2
`
=
−
I
`
.
(6)
The most likely power injection configuration
p
(
`
)
∈
R
n
given the event
|
(
f
ε
)
`
|≥
1
is the solution of the variational
problem
p
(
`
)
=
arg inf
p
∈
R
n
:
|
ˆ
e
>
`
Vp
|≥
1
1
2
(
p
−
μ
)
>
Σ
−
1
p
(
p
−
μ
)
,
(7)
which, when
ν
`
6
= 0
, can be explicitly computed as
p
(
`
)
=
μ
+
(sign(
ν
`
)
−
ν
`
)
σ
2
`
Σ
p
V
>
ˆ
e
`
.
The next proposition shows that the conditional distri-
bution of
p
ε
, given
|
(
f
ε
)
`
|≥
1, gets concentrated around
p
(
`
)
exponentially fast as
ε
→
0, motivating the interpreta-
tion of
p
(
`
)
as the most likely power injection configuration
given the failure of line
`
.
Proposition 2.
Assume that
max
k
=1
,...,m
|
ν
k
|
<
1
, and
that
ν
`
6
= 0
. Then, for al l nodes
i
= 1
,...,n
, and for all
δ >
0
,
lim
ε
→
0
ε
log
P
((
p
ε
)
i
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
∣
∣
|
(
f
ε
)
`
|≥
1)
<
0
.
3
The line power flows corresponding to the power injec-
tion configuration
p
(
`
)
can be calculated as
f
(
`
)
=
Vp
(
`
)
=
ν
+
(sign(
ν
`
)
−
ν
`
)
σ
2
`
V
Σ
p
V
>
ˆ
e
`
∈
R
m
.
We observe that the vectors
p
(
`
)
and
f
(
`
)
are equal to the
conditional expectation of the power injections
p
ε
and
power flows
f
ε
, respectively, conditional on the failure
event
f
`
= sign(
ν
`
), namely
p
(
`
)
=
E
[
p
ε
|
(
f
ε
)
`
= sign(
ν
`
)]
,
(8)
f
(
`
)
=
E
[
f
ε
|
(
f
ε
)
`
= sign(
ν
`
)]
.
In particular, for every
k
= 1
,...,m
,
f
(
`
)
k
=
ν
k
+ (sign(
ν
`
)
−
ν
`
)
Cov(
f
`
,f
k
)
Var(
f
`
)
.
Note that the case
ν
`
= 0 has been excluded only for
compactness. Indeed, in this special case the variational
problem
(7)
has two solutions,
p
(
`,
+)
and
p
(
`,
−
)
. This
can be easily explained by observing that if the power
flow on line
`
has mean
ν
`
= 0, then it is equally likely
for the overload event
{|
f
`
| ≥
1
}
to occur as
{
f
`
≥
1
}
or as
{
f
`
≤ −
1
}
and the most likely power injection
configurations that trigger them can be different.
The previous proposition immediately yields the large
deviations principle also for the first line failure event
‖
f
ε
‖
∞
≥
1, which reads
lim
ε
→
0
ε
log
P
(
||
f
ε
||
∞
≥
1) =
−
min
`
=1
,...,m
(1
−|
ν
`
|
)
2
2
σ
2
`
.
Indeed, the decay rate for the event that at least one
line fails is equal to the minimum of the decay rates for
the failure of each line. The most likely power injections
configuration that leads to the event
‖
f
ε
‖
∞
≥
1 is
p
(
`
∗
)
with
`
∗
= arg min
`
=1
,...,m
(1
−|
ν
`
|
)
2
2
σ
2
`
.
Proof of Proposition 1.
Let (
Z
(
i
)
)
i
∈
N
be a sequence
of i.i.d.
m
-dimensional multivariate normal vectors
Z
(
i
)
∼ N
m
(
ν
,
Σ
f
), and let
S
k
=
1
k
∑
k
i
=1
Z
(
i
)
be the se-
quence of the partial sums. By setting
ε
=
1
k
, it immedi-
ately follows that that
f
ε
d
=
S
k
, where
d
=
denotes equality
in distribution. Denote
g
(
p
) =
1
2
(
p
−
μ
)
>
Σ
−
1
p
(
p
−
μ
).
Following [5, Section 3.D], we get
lim
ε
→
0
ε
log
P
((
f
ε
)
`
≥
1) = lim
k
→∞
1
k
log
P
((
S
k
)
`
≥
1) =
=
−
inf
p
∈
R
n
:
ˆ
e
>
`
Vp
≥
1
g
(
p
) =
−
(1
−
ν
`
)
2
2
σ
2
`
,
(9)
lim
ε
→
0
ε
log
P
((
f
ε
)
`
≤−
1) = lim
k
→∞
1
k
log
P
((
S
k
)
`
≤−
1) =
=
−
inf
p
∈
R
n
:
ˆ
e
>
`
Vp
≤−
1
g
(
p
) =
−
(
−
1
−
ν
`
)
2
2
σ
2
`
.
(10)
The optimizers of problems
(9)
and
(10)
are easily com-
puted respectively as as
p
(
`,
+)
=
μ
+
(1
−
ν
`
)
σ
2
`
Σ
p
V
>
ˆ
e
`
,
p
(
`,
−
)
=
μ
+
(
−
1
−
ν
`
)
σ
2
`
Σ
p
V
>
ˆ
e
`
.
Note that trivially
inf
p
∈
R
n
:
|
ˆ
e
>
`
Vp
|≥
1
g
(
p
) =
= min
{
inf
p
∈
R
n
:
ˆ
e
>
`
Vp
≥
1
g
(
p
)
,
inf
p
∈
R
n
:
ˆ
e
>
`
Vp
≤−
1
g
(
p
)
}
,
and thus identities (6) and (7) immediately follow.
Proof of Proposition 2.
We have
log
P
((
p
ε
)
i
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
∣
∣
|
(
f
ε
)
`
|≥
1)
= log
P
((
p
ε
)
i
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
,
|
(
f
ε
)
`
|≥
1)
−
log
P
(
|
(
f
ε
)
`
|≥
1)
.
Denote
g
(
p
) =
1
2
(
p
−
μ
)
>
Σ
−
1
p
(
p
−
μ
). From large devia-
tions theory, it holds that that
lim
ε
→
0
ε
log
P
(
|
(
f
ε
)
`
|≥
1) =
−
inf
p
∈
R
n
:
|
ˆ
e
>
`
V
p
|≥
1
g
(
p
)
(11)
lim
ε
→
0
ε
log
P
((
p
ε
)
i
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
,
|
(
f
ε
)
`
|≥
1) =
=
−
inf
p
∈
R
n
:
|
ˆ
e
>
`
V
p
|≥
1
,
|
p
i
−
p
(
`
)
i
|≥
δ
g
(
p
)
.
(12)
Define the corresponding decay rates as
I
`
=
inf
p
∈
R
n
:
|
ˆ
e
>
`
Vp
|≥
1
g
(
p
)
, J
`
=
inf
p
∈
R
n
:
|
ˆ
e
>
`
Vp
|≥
1
,
|
p
i
−
p
(
`
)
i
|≥
δ
g
(
p
)
.
Then we can rewrite
lim
ε
→
0
ε
log
P
((
p
ε
)
k
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
∣
∣
|
(
f
ε
)
`
|≥
1) =
−
J
`
+
I
`
,
and, therefore, the claim is equivalent to proving that
J
`
> I
`
. Notice that the feasible set of the minimization
problem
(12)
is strictly contained in that of the problem
(11), implying that
J
`
≥
I
`
.
Recall that
p
(
`
)
is the unique optimal solution of
(11)
,
and let
ˆ
p
(
`
)
be an optimal solution of
(12)
. Clearly
ˆ
p
(
`
)
is feasible also for problem
(11)
. If it was the case that
J
`
=
I
`
, then
ˆ
p
(
`
)
would be an optimal solution for
(11)
,
and thus by uniqueness (
g
(
p
) is strictly convex)
ˆ
p
(
`
)
=
p
(
`
)
. But this leads to a contradiction, since
ˆ
p
(
`
)
is by
construction such that
|
ˆ
p
i
−
ˆ
p
(
`
)
i
|≥
δ
. Hence
J
`
> I
`
and
we conclude that
lim
ε
→
0
ε
log
P
((
p
ε
)
i
/
∈
(
p
(
`
)
i
−
δ,p
(
`
)
i
+
δ
)
∣
∣
|
(
f
ε
)
`
|≥
1)
<
0
.
4
Extension to non-Gaussian case
In this section we briefly describe how to extend the
analyis to the non-Gaussian scenario. Consider a model
for the power injection vector given by
p
ε
=
μ
+
√
ε
X
,
where
μ
∈
R
n
and
X
= (
X
1
,...,X
n
) is a random vector
with mean 0 and log-moment generating function
log
M
(
s
) = log
E
[
e
〈
s
,
X
〉
]
.
The power flows vector is thus given by
f
ε
=
V p
ε
. Define
the Fenchel-Legendre (also known as the convex conju-
gate) transform of log
M
(
s
), i.e.
Λ
∗
(
x
) = sup
s
∈
R
n
(
〈
s
,
x
〉−
log
M
(
s
))
.
Then, for every
`
= 1
,...,m
, the sequence (
f
ε
)
ε>
0
satisfies
the large deviations principle (see [6])
lim
ε
→
0
ε
log
P
(
|
(
f
ε
)
`
|≥
1) =
−
inf
x
∈
R
n
:
|
ˆ
e
>
`
V
(
μ
+
x
)
|≥
1
Λ
∗
(
x
)
,
and the most likely power injection configuration
p
(
`
)
∈
R
n
given the event
|
(
f
ε
)
`
|≥
1 is
p
(
`
)
=
μ
+
arg inf
x
∈
R
n
:
|
ˆ
e
>
`
V
(
μ
+
x
)
|≥
1
Λ
∗
(
x
)
.
The rest of the analysis can then be carried out along
similar lines as we did for the Gaussian case.
POWER FLOW REDISTRIBUTION
For every line
`
define
J
(
`
) to be the collection of lines
that fail jointly with
`
as
J
(
`
) =
{
k
:
|
f
(
`
)
k
|≥
1
}
.
Let
j
(
`
) =
|J
(
`
)
|
be its cardinality and note that
j
(
`
)
≥
1
as trivially
`
always belongs to
J
(
`
). Denote by
̃
G
(
`
)
the
graph obtained from
G
by removing all the lines in
J
(
`
).
Let us focus first on the case of the isolated failure
of line
`
, that is when
J
(
`
) =
{
`
}
. In this case
̃
G
(
`
)
=
G
(
V,E
\{
`
}
) is the graph obtained from
G
after removing
the line
`
= (
i,j
). Provided that the power injections
remain unchanged, the power flows redistribute among the
remaining lines. Using the concept of
effective resistance
matrix
R
∈
R
n
×
n
and under the DC approximation, in [
7
–
9
] it is proven that alternative paths for the power to flow
from node
i
to
j
exist (i.e.,
̃
G
(
`
)
is still connected) if and
only if
β
i,j
R
i,j
6
= 1. In other words,
β
i,j
R
i,j
= 1 can only
occur in the scenario where line
`
= (
i,j
) is a
bridge
, i.e.,
its removal results in the disconnection of the original
graph
G
in two components. If
̃
G
(
`
)
is still a connected
graph, the power flows after redistribution
̄
f
(
`
)
∈
R
m
−
1
are related with the original line flows
f
∈
R
m
in the
network
G
by the relation
̄
f
(
`
)
k
=
f
k
+
f
(
`
)
`
φ
(
`
)
k
,
for every
k
6
=
`,
where
f
(
`
)
`
=
±
1 depending on the way the power flow on
line
`
exceeded the threshold 1. If
`
= (
i,j
) and
k
= (
a,b
)
the coefficient
φ
(
`
)
k
∈
R
can be computed as
φ
(
`
)
k
=
φ
(
i,j
)
,
(
a,b
)
=
β
k
·
C
`
C
k
·
R
a,j
−
R
a,i
+
R
b,i
−
R
b,j
2(1
−
β
`
R
i,j
)
,
(13)
The ratio
C
`
/C
k
appears in the latter formula since we
work with normalized line power flows and we correspond-
ingly defined
φ
(
`
)
=
{
φ
(
`
)
k
}
k
6
=
`
to be the normalized ver-
sion of the classical line outage distribution factors (LODF,
[
10
]). Moreover, we define the
most likely
power flows
configuration
̃
f
(
`
)
∈
R
m
−
1
after redistribution as
̃
f
(
`
)
k
=
f
(
`
)
k
+
f
(
`
)
`
φ
(
`
)
k
,
for every
k
6
=
`.
(14)
Ring topology
We now focus on a particular topology, namely the
ring on
n
nodes, which we use as an illustrative example
to show the non-locality of cascades in the Main Body.
In this topology, nodes are placed on a ring and each
node is connected to its previous and subsequent neigh-
bor. Denote the set of nodes as
N
=
{
1
,...,m
}
and
the set of lines
L
=
{
l
1
,...,l
n
}
, where
l
1
= (
n,
1)
,l
2
=
(1
,
2)
,...,l
n
= (
n
−
1
,n
). It is easy to prove that, in a ring
network with homogeneous line thresholds and unitary
susceptances,
φ
`,k
=
−
1 for every
`
6
=
k
.
Lemma 1.
Consider a ring network with homogeneous
line thresholds (
C
`
=
C
for every line
`
) and homogeneous
unitary susceptances (
β
`
= 1
for every line
`
). Then
i)
The effective resistance between a pair of nodes
i,j
is given by
R
i,j
=
|
j
−
i
|
(
n
−|
j
−
i
|
)
n
.
(15)
ii)
For every pair of lines
l
k
= (
k
−
1
,k
)
,l
`
= (
`
−
1
,`
)
,
with
k
6
=
`
, the LODF is constant and equal to
φ
(
k
−
1
,k
)
,
(
`
−
1
,`
)
=
−
1
.
Proof.
i) See identity (4) in [
11
]. ii) First, observe that
the effective resistance between two adjacent nodes
i
and
j
=
i
+ 1 in a circuit graph is equal to
R
i,j
=
n
−
1
n
, thanks
to Eq.
(15)
. After a straightforward calculation, and using
that
β
`
= 1
,C
`
=
C
for all lines
`
, Eq. (13) becomes
φ
(
k
−
1
,k
)
,
(
`
−
1
,`
)
=
−
2
2
n
(
1
−
n
−
1
n
)
=
−
1
.
5
General topology
Going back to the case of a general network topology
and any type of failures, isolated or joint, the power flows
after redistribution
̄
f
(
`
)
∈
R
m
−
j
(
`
)
are related with the
power injections
p
∈
R
n
by the relation
̄
f
(
`
)
=
̃
V
(
`
)
p
,
where the (
m
−
j
(
`
))
×
n
matrix
̃
V
(
`
)
can be constructed
analogously to
V
, but considering the altered graph
̃
G
(
`
)
instead of
G
. We define the
most likely
power flow con-
figuration
̃
f
(
`
)
after redistribution as
̃
f
(
`
)
=
̃
V
(
`
)
p
(
`
)
,
which generalizes Eq.
(14)
to any kind of failure, isolated
of joint. The next proposition shows that it is enough
to look at the vector
̃
f
(
`
)
to determine whether a line
that survived at the first cascade stage (i.e., that did not
fail jointly with
`
) will fail with high probability or not
after the power redistribution (i.e., at the second cascade
stage).
Proposition 3.
Assume that
max
k
=1
,...,m
|
ν
k
|
<
1
, and
that
ν
`
6
= 0
. Then, for al l lines
k
∈
E
\J
(
`
)
, and for al l
δ >
0
,
lim
ε
→
0
ε
log
P
((
̄
f
(
`
)
ε
)
k
/
∈
(
̃
f
(
`
)
k
−
δ,
̃
f
(
`
)
k
+
δ
)
∣
∣
|
(
f
ε
)
`
|≥
1)
<
0
.
In particular, if
|
̃
f
(
`
)
k
|≥
1
, then
P
(
|
(
̄
f
(
`
)
ε
)
k
|≥
1
∣
∣
|
(
f
ε
)
`
|≥
1)
→
1
as
ε
→
0
,
exponential ly fast in
1
/ε
.
Proof.
Let
A
ε,k
denote the event
A
ε,k
=
{
(
̄
f
(
`
)
ε
)
k
/
∈
(
̃
f
(
`
)
k
−
δ,
̃
f
(
`
)
k
+
δ
)
}
, and define
Q
=
lim
ε
→
0
ε
log
P
(
A
ε,k
∣
∣
|
(
f
ε
)
`
|≥
1). The proof that
Q <
0 is
analogous to the proof of Prop. 2. For the second part,
it follows from
lim
ε
→
0
ε
log
P
(
A
ε,k
∣
∣
|
(
f
ε
)
`
|≥
1) =
Q
that
for every
η >
0 there exists a
̄
ε
such that, for every
ε <
̄
ε
,
Q
−
η
≤
ε
log
P
(
A
ε,k
∣
∣
|
(
f
ε
)
`
|≥
1)
≤
Q
+
η,
and thus
exp
(
Q
−
η
ε
)
≤
P
(
A
ε,k
∣
∣
|
(
f
ε
)
`
|≥
1)
≤
exp
(
Q
+
η
ε
)
.
Let
A
c
ε,k
denote the complementary event of
A
ε,k
. Since
|
̃
f
(
`
)
k
|≥
1, for
δ
sufficiently small we have
A
c
ε,k
=
{|
(
̄
f
(
`
)
ε
)
k
−
̃
f
(
`
)
k
|
< δ
}⊆{|
(
̄
f
(
`
)
ε
)
k
|≥
1
}
,
yielding
P
(
|
(
̄
f
(
`
)
ε
)
k
|≥
1
∣
∣
|
(
f
ε
)
`
|≥
1)
≥
P
(
A
c
ε,k
∣
∣
|
(
f
ε
)
`
|≥
1)
≥
1
−
exp
(
Q
+
η
ε
)
.
Since
Q <
0, the result follows.
APPLICATION: SCIGRID GERMAN NETWORK
We now demonstrate our methodology in the case of a
real-world power grid and a realistic system state.
Dataset Description
We perform our experiments using PyPSA, a free soft-
ware toolbox for power system analysis [
12
]. We use the
dataset described in [
13
,
14
], which provides a model of
the German electricity system based on SciGRID and
OpenStreetMap [15, 16].
The dataset includes load/generation time series and
geographical locations of the nodes, differentiating be-
tween renewable and conventional generation. It also
provides data for transmission lines limits, transform-
ers, generation capacity and marginal costs, allowing us
to couple our theoretical analysis with realistic Optimal
Power Flow (OPF, [
17
]) computations. The time-series
provide hourly data for the entire year 2011. For more
technical information on the dataset, we refer to [
13
,
14
].
The SciGRID German network consists of 585 buses,
1423 generators including conventional power plants and
wind and solar parks, 38 pump storage units, 852 lines
and 96 transformers. For the analysis carried out in
this paper, storage units are not included and we ex-
clude transformer failures. The renewable generators
are divided in three classes, solar, wind onshore and
wind offshore. Each bus can house multiple generators,
both renewable and conventional, but it is limited to
at most one renewable generator for each class. Let
N
w.off
,
N
w.on
,
N
sol
denote the set of buses housing, re-
spectively, wind offshore, wind onshore and solar genera-
tors, with
|N
w.off
|
= 5
,
|
N
w.on
|
= 488
,
|
N
sol
|
= 489, and
N
w.off
⊆ N
w.on
⊆ N
sol
. The remaining 96 buses house
441 conventional generators.
Let
n
s
= 489 denote the total number of buses housing
renewable generators. If a bus houses both renewable and
conventional generators, it will be considered a
stochas-
tic
bus for our decomposition formulation. We model
stochastic net power injections by means of a multivariate
Gaussian random vector
p
s
∼N
n
s
(
μ
s
,ε
Σ
p
).
The distinction between the noise parameter
ε
and the
covariance matrix
Σ
p
is relevant only for the theoretical
analysis (where we take the limit
ε
→
0 while the matrix
Σ
p
is fixed), since as far as the numerical case study is
concerned, all the results are obtained by using the prod-
uct
ε
Σ
p
, which is directly estimated from the SciGRID
data. In the following, we will thus take
ε
= 1 and refer
to the covariance matrix of
p
s
simply as
Σ
p
.
6
Data-based model for
μ
s
In order to get a realistic nominal line flows value
ν
,
we perform a linear OPF relative to the day 01
/
01
/
2011,
for different hours of the day. A linear OPF consists of
minimizing the total cost of generation, subject to energy
balance, generation and transmission lines constraints, un-
der the assumptions of the DC approximations. In order
to model a heavily-loaded but not overloaded system, in
the OPF we scale the true line limits
C
`
by a contingency
factor of
λ
= 0
.
7. This is a common practice in power
engineering that allows room for reactive power flows and
stability reserve.
More precisely, let
g
(
i
) be the generation at bus
i
as
outputted by the OPF for a given hour, and let us write it
as
g
(
i
) =
g
r
(
i
) +
g
d
(
i
), with
g
r
(
i
) the power produced by
renewable generators attached to the bus, and
g
d
(
i
) the
power supplied by conventional generators. If the demand
at bus
i
is given by
d
(
i
), then the average stochastic power
injection vector
μ
s
∈
R
n
s
is modeled as
(
μ
s
)
i
=
g
r
(
i
) +
g
d
(
i
)
−
d
(
i
)
, i
= 1
,...,n
s
,
while the deterministic power injection reads
p
i
=
g
d
(
i
)
−
d
(
i
) for
i
∈I
d
.
Data-based model for Σ
p
In order to model the fluctuations of renewable genera-
tion around the nominal values, and thus estimate
Σ
p
, we
use realistic hourly values of wind and solar energy pro-
duction to fit a stochastic model. We then use the steady-
state covariance of the model residuals as an estimate
for
Σ
p
. Following [
18
], we choose to use AutoRegressive-
Moving-Average (ARMA) models, which we describe in
details below.
Note that we do not aim to find the best possible
stochastic model for renewable generation, which is be-
yond the scope of this paper, but instead to provide an
estimate for the covariance matrix
Σ
p
in order to validate
our theoretical results, which are asymptotically valid in a
small-noise regime. We speculate that more sophisticated
models, and/or data on smaller time-scales, may lead to
smaller values for the correlations in
Σ
p
, thus getting
closer to the small noise limit.
We now describe the estimation procedure for
Σ
p
(as
mentioned before we normalize
ε
= 1 in our empirical
study). The SciGRID dataset contains time series
y
w.off
∈
R
M
×
5
,
y
w.on
∈
R
M
×
488
,
y
sol
∈
R
M
×
489
,
for the available power output of wind offshore, wind
onshore and solar generators, for each hour of the year
2011, accounting for a total of
M
= 8760 measurements
for each generator [
13
]. For each time series,
y
(
·
)
(
t,j
)
denote the available power output at time
t
for the
j
-th
generator of a given type, in MW units.
Wind power model
As a pre-processing step, we merge together the two
time series
y
w.off
,
y
w.on
by summing up the onshore and
offshore wind power at the buses
N
w.off
⊆ N
w.on
. This
yields the time series of wind power production
y
w
(
t,j
) =
y
w.on
(
t,j
) +
1
j
∈N
w.off
y
w.off
(
t,j
)
,
where
1
{}
is the indicator function of the event in the
bracket, taking value 1 if the event is satisfied, and 0,
otherwise.
We select one portion of the data
{
1
,...,T
} ⊆
{
1
,...,M
}
, corresponding to the month of January, to be
used to fit the model. For each windpark
j
, following [
18
]
we consider an ARMA(1,24) model of the form
x
(
t,j
) =
a
1
,j
x
(
t
−
1
,j
) +
e
(
t,j
)
+
m
1
,j
e
(
t
−
1
,j
) +
...
+
m
24
,j
e
(
t
−
24
,j
)
,
where
x
(
t
−
1
,j
) is the auto-regressive term, and
e
(
t
−
k,j
),
k
= 1
,...,
24, are the white-noise error terms. For each
windpark
j
, we fit the above model to the wind power data
{
y
w
(
t,j
)
}
t
=1:
T
in R using the function
arima
, and con-
sider the time series of the residuals
e
w
(1
,j
)
,...,e
w
(
T,j
).
The empirical variance of the residuals is used as proxy
for the variance of the output of windpark
j
, namely
(
Σ
w
)
jj
=
̂
Var(
e
w
(1
,j
)
,...,e
w
(
T,j
))
,
where
̂
Var
denotes the empirical variance. In a similar
way, the empirical covariance of the residuals is used to
model the covariance between the output of windparks
i
and
j
,
i
6
=
j
, namely
(
Σ
w
)
ij
=
̂
Cov
(
{
e
w
(
t,i
)
}
t
=1:
T
,
{
e
w
(
t,j
)
}
t
=1:
T
)
,
where
̂
Cov denotes the empirical covariance.
Solar power model
State-of-the-art models for solar irradiance often com-
bine statistical techniques with cloud motion analysis and
numerical weather prediction (NWP) models, see [
19
] for
a review. Since the available data in our case study are
limited to historical records for power production of solar
generators, and do not include any weather data, we used
the purely statistical model ARMA(
p
,
q
), which has been
used succesfully in [20].
Regarding the orders
p,q
of the ARMA model, af-
ter some exploratory analysis we decided to use an
ARMA(24,24) model with all parameters fixed to 0, ex-
cept for the ones corresponding to the seven hours before,
and the one corresponding to twenty-four hours before.
The rationale behind this choice is that by using the value