OPTIMIZATION STRATEGIES FOR THE VULNERABILITY
ANALYSIS OF THE ELECTRIC POWER GRID
∗
ALI PINAR
†
, JUAN MEZA
‡
, VAIBHAV DONDE
§
, AND BERNARD LESIEUTRE
¶
Abstract. Identifying small groups of lines, whose removal would cause a severe blackout, is
critical for the secure operation of the electric power grid. We show how power grid vulnerability
analysis can be studied as a mixed integer nonlinear programming (minlp) problem. Our analysis
reveals a special structure in the formulation that can be exploited to avoid nonlinearity and ap-
proximate the original problem as a pure combinatorial problem. The key new observation behind
our analysis is the correspondence between the Jacobian matrix (a representation of the feasibility
boundary of the equations that describe the flow of power in the network) and the Laplacian matrix
in spectral grap h theo ry (a represe ntation of the graph of the power grid). The reduced combinatorial
problem is known as the network inhibition problem, for which we present a mixed integer linear
programming formulation. Our experiments on benchmark power grids show that the reduced com-
binatorial model provides an accurate approximation, to enable vulnerability analyses of real-sized
problems with more than 10,000 power lines.
Key words. mixed integer nonlinear programmi ng, network inhibition, network flow, mixed
integer linear programming, electric power flow, network vulnerability, graph theory
AMS subject classifications. 90C11, 90C27, 90C90, 90C30
1. Introduction. Robust operation of a power grid requires anticipation of com-
ponent outages that could lead to dramatic blackouts. The current practice is to check
for single contingencies to ensure the system stays intact after a single line outage.
However, a small number of line outages (e.g., 3–5) can cause catastrophic blackouts,
as evidenced by the Northeast Blackout in August 2003. In this article, we consider
the power network vulnerability analysis problem, which aims to find small groups of
lines, whose loss can cause a severe blackout. Specifically, we pose the following two
related optimization problems: 1) compute the minimum number of line failures that
will cause a damage of at least a specified severity and 2) compute a combination of
a specified numb er of lines, whose loss will cause the maximum damage.
We consider the problem in a static sense by examining the relation between the
operating point, which describes the current generation and consumption at each node
in the network, and the feasibility boundary of the power flow equations. The severity
of the events we identify could be different when dynamics and cascading events are
considered. Our main focus here, therefore, is to identify simple events that can trigger
a severe blackout, not to analyze its consequences, which requires solving differential
algebraic equations with discrete variables. Cascading events start with a significant
disturbance that forces system elements to operate beyond their capabilities. For
this reason, we look for minimal changes in the network topology that push the
∗
This work was supported by the Director, O ffic e of Science, Division of Mathematical, In-
formation, and Computational Sciences of U.S. Department of Energy unde r contract DE-AC03-
76SF00098.
†
Correspon ding author. High Performance Computing Research Department, Lawrence Berkeley
National Laboratory (apinar@lbl.gov).
‡
High Performance Computing Research Department, Lawrence Berkeley National Laboratory
(JCMeza@lbl.gov).
§
Corporate research center, ABB Inc., Raleigh NC (vaibhav.d.donde@us.abb.com). Work per-
formed at Lawrence Berkeley National Laboratory.
¶
Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory
and Electrical and Computer Engineering Department, University of Wisconsin, Madison.
(BCLesieutre@lbl.gov)
1
Work also completed under Department of Energy Contract No. DE-AC02-05CH11231.
2 Pinar, Meza, Donde, and Lesieutre
current operating point significantly outside of the feasibility region of the power
flow equations. This problem statement leads to a bi-level optimization problem,
since we are looking for minimal changes in network topology that maximize the
distance between the current operating point and the new feasibility region. Moreover,
the problem combines nonlinearity due to the power flow equations, with discrete
variables, due to changes in the network topology.
In this article, we propose a mixed integer nonlinear programming (minlp) formu-
lation for the power network vulnerability analysis problem. To measure the severity
of the disturbance to the system, we use a load shedding mechanism, which opti-
mally decreases the generation and consumption in the system to restore feasibility.
We avoid solving nested optimization problems by replacing the inner optimization
problem that compute the distance between the current operating point and the new
feasibility region, with its Karush-Kuhn-Tucker conditions. Next, we analyze the
structure of a feasible solution to our minlp formulation to reveal a special structure
that can b e exploited to reduce the problem to a pure combinatorial problem. We
show that at a feasible solution to our minlp formulation, the power network will be
divided into two groups: one with excess generation and one with excess load, and
the optimal load shedding strategy requires that in the load-rich region, we decrease
only the consumption and keep the generation as is. Similarly in the generation-rich
region, we decrease only the generation and keep the consumption as is. Moreover,
we prove that at least one line that connects these two regions works at its maxi-
mum capacity to transfer power from the generation-rich side to the other. This clear
combinatorial structure of a feasible solution me ans that an optimal solution seeks a
decomposition with maximum load/generation mismatch and minimum transmission
capability between the the two regions. This observation leads to our major result:
the original minlp problem can be reduced, after some realistic simplifications, to
a pure combinatorial problem, namely the network inhibition problem. With this
reduction, we directly seek the values of discrete variables in the formulation with-
out solving the nonlinear equations, simplifying the problem complexity both in a
theoretical and practical sense.
Identification of multiple contingencies has recently drawn much interest both
from the optimization and power systems communities. Salmeron, Wood, and Baldick
[25] employed a linearized power flow model and used a bilevel optimization framework
along with mixed-integer programming to analyze the security of the electric grid. The
critical elements of the grid were identified by maximizing the long-term disruption in
the power system operation. The bilevel optimization framework has also been used
by Arroyo and Galiana [18]. In all of these formulations the optimization framework
appears promising for such types of problems where the critical system elements
that make the system vulnerable to failures must be identified. Donde et al. [12],
prop os ed a method that connected the feasibility boundary of power flow equations
with spectral graph theory, when voltages are fixed at their nominal values, and only
active power flow constraints are considered. Later, Donde et al. [11] extended
their approach to include reactive p ower and proposed a mixed integer nonlinear
programming formulation to identify the most significant blackout that can be caused
by a specified number of lines or to identify the minimum number of lines to cause
a blackout of specified severance. More recently, Lesieutre et al. [19, 20] approached
this problem from a graph theoretical perspective, by looking for subgraphs in a
given graph that are loosely connected to the rest of the graph and have a significant
load/generation mismatch. Grijalva and Sauer [15, 16] related topological cuts in the
Vu lnera bility Analysis of the Power Grid 3
power network with the static collapse based on branch complex flows. He et. al. [17]
used a voltage stability margin index to identify weak locations in a power network.
Bienstock and Mattia used the direct current power flow model and mixed integer
linear programming to find the most cost-effective way to increase edge capacities
to avoid cascading outages for a given set of failure scenarios [3]. Oliviera et al.
have used similar models and techniques to study how to add power lines to improve
system resilience [21]. In addition to these largely static analyses, system dynamics for
cascading events has also drawn a lot of interest. In [4,6,8] Dobson et al. used a long-
term model of the grid to study how failure of a component affects other components
in the system, to reveal failure s tatistics consistent with those observed in the power
grid. The same authors have also studied probabilistic models with the aim to better
understand cascade propagation [5,9, 10].
The remainder of this article is organized as follows. Section 2 reviews matrix
representations of graphs and the basics of sp ec tral graph theory that are relevant
to this article. In Section 3, we present a minlp formulation for the power network
vulnerability analysis problem. The structure of a feasible solution to this problem
and how this structure can be exploited to reduce the minlp formulation to a pure
combinatorial problem are discussed in Section 4. We describe the network inhibition
problem and its integer programming formulation in Section 5. Section 6 presents our
experimental results, and we conclude with Section 7.
2. Graphs and matrices. Matrix representations of graphs have long been
used to apply algebraic techniques to analyze graphs. Here we review the node-arc
incidence matrix and the Laplacian matrix, as two of the commonly used representa-
tions for graphs. The node-arc incidence matrix of a graph is used in flow problems,
and we will use this representation to present power flow equations. The Laplacian
matrix for graphs on the other hand, underlies spectral graph theory, which can be
used to analyze the connectedness of graphs. Let G = (V, E) be a graph with n
vertices and m edges. We use (v
i
, v
j
) to denote an edge that goes from vertex v
i
to
vertex v
j
. The node-arc incidence matrix, A, of this graph is an m × n matrix, where
the j-th column of A represents the j-th vertex, v
j
, and the i-th row represents the
the i-th edge, e
i
, in G. Each row has only two nonzeros at the columns that represent
the e nd vertices of the respective edge. T he entry is -1 or 1, depending on whether
the respective edge is directed from or to the corresponding vertex, respectively. For-
mally, we use a
ij
to denote the matrix entry at the i-th row and the j-th column of
A, which is defined as follows.
a
ij
=
−1 if e
i
= (v
j
, u) ∈ E
1 if e
i
= (u, v
j
) ∈ E
0 otherwise
The node-arc incidence matrix A of the graph in Fig. 2.1 is as follows.
A =
−1 1
−1 1
1 −1
−1 1
−1 1
−1 1
−1 1
−1 1
4 Pinar, Meza, Donde, and Lesieutre
Fig. 2.1. A sample directed graph
The Laplacian of a graph G = (V, E) is an n × n matrix, where each row and
column represents a vertex in the graph. The diagonal entry is equal to the degree of
the associated vertex. An off-diagonal entry is -1, if the associated vertices of the row
and column are connected in the graph, and 0 otherwise. Formally, let d
i
denote the
degree of vertex v
i
, and let l
ij
denote the entry of the Laplacian matrix at the i-th
row and the j-th column, which we define as follows.
l
ij
=
d
i
if i = j
−1 if (v
i
, v
j
) ∈ E or (v
j
, v
i
) ∈ E
0 otherwise
The Laplacian of the graph in Fig. 2.1 is
L =
2 −1 −1
−1 4 −1 −1 −1
−1 −1 3 −1
−1 2 −1
−1 −1 3 −1
−1 −1 2
.
We note that L can also be defined as
L = A
T
A, (2.1)
where A is the node-arc incidence matrix of the graph. This property holds regardless
of the directions of edges in G. It is possible to add edge weights to the definition of
Laplacian of a graph. In this case, the diagonal entry becomes the sum of weights of
edges adjacent to the respective vertex, as oppose d to the degree of this vertex, and
the negative of the edge weight replaces “-1” as the off-diagonal entries. In this case,
Eq. (2.1) can be rephrased as
L
w
= A
T
D
w
A, (2.2)
where D
w
is a diagonal matrix so that the i-th diagonal is the weight of edge e
i
, and
L
w
is the weighted Laplacian. Observe that a zero diagonal entry on D
w
corresponds
to removing a line from the graph.
Vu lnera bility Analysis of the Power Grid 5
The Laplacian of a graph is the basic element of spectral graph theory. Let
λ
0
≤ λ
1
≤ . . . ≤ λ
n−1
be the eigenvalues of L. The Laplacian matrix is symmetric and
semi-definite, and thus all eigenvalues are real and nonnegative. It is easy to see that
λ
0
= 0, since all rows and columns of L add up to zero, and thus the vector, e, whose
entries are all the same and nonzero, is a singular vector for L. The smallest nontrivial
eigenvalue λ
1
is more interesting due to its applications. Fiedler called λ
1
the algebraic
connectivity of G [13], as it provides a metric for the connectedness of a graph. If
the graph inherently involves two loosely coupled sub-graphs, then λ
1
will be small.
Fiedler also proved that λ
1
will decrease as we remove edges from the original graph,
and it will be zero when the graph is decoupled into two disconnected components. A
fundamental result in spectral graph theory generalizes this observation so that the
multiplicity of the eigenvalue 0 gives the number of connected components in G.
Lemma 2.1. Let L be the Laplacian of graph G, and let λ
0
≤ λ
1
≤ . . . ≤ λ
n−1
be its eigenvalues. If λ
i
= 0 and λ
i+1
6= 0, then G has exactly i + 1 connected
components.
The multiplicity of eigenvalue 0 determines the number of connected components
in a graph, associated eigenvectors identify these connected components. For an
eigenvalue λ
i
= 0, the corresponding eigenvector v
i
has the sam e value for all vertices
in a component, and a different value for each one of the i + 1 components. This
result underlies our analysis of the structure of an optimal solution in Section 4.2.3.
3. Problem formulation. Our focus in this work is to identify simple events
that can trigger a c asc ading event, not to analyze consequences of cascading. Cas-
cading events start with a significant disturbance to the system, and continue with
failures of other system components, as these components are pushed beyond their
capabilities, while the system is trying to avert a blackout. It will be the initial sig-
nificant disturbance that we seek in this work, and thus we focus on static p ower flow
analysis. Below, we first describe our power flow model and then describe how we
measure the significance of an event. Finally, we cast the power grid vulnerability
problem as a mixed integer nonlinear programming (minlp) problem. In [11], a simi-
lar formulation is presented for a full power flow model with active and reactive power
equations, and a slightly different load-shedding model.
3.1. Power system model. We consider a loss les s power system network with
m buses (nodes) and n lines (edges). We assume the voltages at the buses are fixed,
and thus the dependence of real power injections at buses on the phase angle vari-
ables θ can be fully described by active power constraints, making the reactive power
constraints unnecessary. The power flowing through the lines can be expressed as
P
line
= B sin(Aθ),
where P
line
is a vector of power flows over the lines, B is a diagonal matrix whose
diagonal entries correspond to line admittances, A is a node-arc incidence matrix that
represents the power network, and sin(Aθ) denotes a vector whose i-th component
is sin((Aθ)
i
). A vector of powe r injections P is then obtained by adding the power
flowing out of the buses into the network.
A
T
B sin(Aθ) − P = 0, (3.1)
with Aθ taking values between −π/2 and π/2, as required for steady state stability.
Here, we will work with a given topology of the power grid and inve stigate the
endurance of the grid to changes in topology. To extend the power flow equations
