Probabilistic N-k Failure-Identification for Power Systems
Kaarthik Sundar, Carleton Coffrin, Harsha Nagarajan, Russell Bent∗
Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM

arXiv:1704.05391v8 [] 19 Mar 2018

Abstract
This paper considers a probabilistic generalization of the N-k failure-identification problem in power transmission networks, where the
probability of failure of each component in the network is known a priori and the goal of the problem is to find a set of k components that
maximizes disruption to the system loads weighted by the probability of simultaneous failure of the k components. The resulting problem
is formulated as a bilevel mixed-integer nonlinear program. Convex relaxations, linear approximations, and heuristics are developed to
obtain feasible solutions that are close to the optimum. A general cutting-plane algorithm is proposed to solve the convex relaxation and
linear approximations of the N-k problem. Extensive numerical results corroborate the effectiveness of the proposed algorithms on small-,
medium-, and large-scale test instances; the test instances include the IEEE 14-bus system, the IEEE single-area and three-area RTS96
systems, the IEEE 118-bus system, the WECC 240-bus test system, the 1354-bus PEGASE system, and the 2383-bus Polish winter-peak
test system.
Keywords— power network resilience, network interdiction, (N-k) vulnerability, nonlinear optimization, convex relaxation, cuttingplane algorithm, network flow, DC power flow, AC power flow

1.

Stackelberg game (see [3]) with an attacker and a defender
where the attacker’s and defender’s actions are sequential
and the attacker has a perfect model of how the defender
will optimally operate the transmission system. The objective of the attacker is to identify exactly k components
to maximize the impact to the transmission system, measured in terms of the minimum load that must be shed by
the defender weighted by the probability of the simultaneous or near-simultaneous failure of the k components.
Hence, N and k refer to the total number of components
and the number of them that can be interdicted in the system, respectively. The number of possible N-k contingencies, even for small values of k, makes total enumeration
computationally intractable.

Introduction

In modern society, the electric power transmission grid
plays a crucial role in sustaining the socioeconomic systems that civilization depends on. While this dependence
is widely observed and discussed, large-scale failures (e.g.,
Superstorm Sandy, Hurricane Matthew, cyber-attacks) remain a threat. These historical events highlight the need
for methods that identify sets of components whose failure leads to significant impacts. One model that is used
for these identifications is the N-k interdiction or the N-k
failure-identification problem [1], which focuses on identifying an N-k contingency, i.e., a set of k critical components of the transmission system whose simultaneous or
near-simultaneous failure would maximize the disruption,
measured in terms of the amount of load shedding in the
system. This paper considers a probabilistic generalization of the N-k problem, which we will refer to as the
PNK, where the probability of failure of each component
is known a priori and the objective is to find a set of probable k components that causes the maximum disruption
to the transmission system under a simultaneous failure.
The probability of failure of each component can be generated using fragility models based on exposure to extreme
events [2], among other things.
It can be useful to consider the problem as a bilevel
∗ Corresponding

1.1.

Related work

The literature contains many variants of the N-k problem, and algorithms that can obtain optimal solutions and
heuristics have been extensively studied for these variants.
For ease of exposition, we will analyze the work done
in the literature using the following three categories: (1)
deterministic or stochastic variant, (2) linear or convex
nonlinear representation of the physics governing the flow
of power in electric transmission systems, and (3) heuristic
or exact approach to solving the corresponding problem.
The physics governing the flow of power through elec-

author: rbent@lanl.gov

1

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

tric transmission networks is nonlinear and is given by
Kirchoff’s laws; the equations that govern the physics
of power flow are referred to as the alternating current
(AC) power flow equations. Typically, problems relating
to the failure of transmission lines and buses in an electric
transmission system use a linear approximation of the AC
power flow equations; the equations obtained after linear
approximation are called the linearized direct current (DC)
power flow equations.
We will start by reviewing the literature on the deterministic N-k interdiction problem that uses the DC power
flow equations to develop exact approaches. We remark
that the deterministic N-k interdiction problem and the
deterministic N-k failure-identification problem are the
the same; hence, we shall use the two names interchangeably. To the best of our knowledge, the first work in the
literature to develop mixed-integer programming-based
formulations for the deterministic variant of the PNK using the DC power flow model is by Salmeron et al. [4].
The authors in [4] develop a mixed-integer bilevel program to model the deterministic interdiction of lines and
buses in the transmission system. A Benders decomposition algorithm, based on the work in [4], is developed
to solve the same problem in [5] on small test instances.
The authors in [5] also perform AC analysis using PowerWorld [6] on the interdiction plan computed using the
DC power flow models; their AC analysis involves computing the percentage deviation of the real power flow
on the lines provided by the DC power flow model and
PowerWorld. The lower-level problem in [4] is replaced
by its dual in [7] and is approximated using KKT conditions in [8]. The first work to systematically develop a
decomposition algorithm to solve a deterministic bilevel
N-k problem based on Benders decomposition on large
test instances is [1]. The algorithm developed in [1], an
extension of the work in [4], is tested on a “U.S. Regional
Grid” with 5000 buses, 5000 lines, and 500 generators. The
works in [1, 5] measure disruption in terms of long-term
power shedding, ignoring short-term shedding resulting
from cascading outages immediately after the attack, and
use a DC power flow model to compute this long-term
power shedding. More recently, [9] develops computationally efficient algorithms to solve a minimum cardinality
variant of the N-k problem, where the objective is to find
a minimum cardinality attack with a throughput less than
a prespecified bound; the DC power flow model is used
to model the physics of the transmission system. They
2

also formulate and solve a nonlinear continuous version
of the problem where the attacker is allowed to change the
transmission line parameters to disrupt the system instead
of removing transmission lines. The authors in [10] study
the deterministic N-k problem where the system operator
is allowed to use both load shedding and line switching
as defensive operations via a Benders decomposition algorithm. In the context of this literature, our approach
makes two key contributions. We analyze the theoretical
and computational effects of using better representations
of the AC power flow equations, viz., convex nonlinear
relaxations, and we consider a probabilistic generalization.
The first work in the literature to use the full AC power
flow model to analyze vulnerability in power systems is
[11]; the work develops two continuous optimization models where the attacker is allowed to increase the impedance
of transmission lines, similar to [9]. They develop a Frank–
Wolfe decomposition algorithm to compute an optimal
solution to the problem. Unlike this paper, the authors
in [11] do not model the removal of any component in
the transmission network; modeling component removal
would require the introduction of binary variables and
possibly disjunctions that require different solution techniques. A different line of work in the literature (see
[12, 13]) focuses on vulnerability assessment in power
grids using a nonlinear approximation of the AC power
flow equations where the problems are formulated as
bilevel mixed-integer programs. These optimization problems either maximize unmet demands or minimize the
number of lines to attack. These two approaches were
shown to be equivalent in [14]. The authors in [13] also
develop an algorithm based on graph partitioning to compute the best N-k attack based solely on the generation
capacities, the demands, and the transmission line limits
while ignoring the power flow constraints; they develop
solution methodologies that can compute the best N-k
attack on a transmission system in the western United
States.
Other heuristic approaches in the literature for the deterministic interdiction problem include [15]; they present
a greedy algorithm that uses a DC power flow model. Although the heuristic approaches produced favorable scaling, the solutions were not optimal. Heuristic approaches
based on genetic algorithms can be found in [16].
Another related problem of contingency identification
has received considerable attention in the literature (see
[17, 18, 19, 20] and references therein). In particular, [18]

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

develops a heuristic approach to identify multiple contingencies that can initiate cascades on large transmission
systems, and [19] uses current injection-based methods
and introduces “line outage distribution factors” to identify contingencies. This work in the literature focuses on
using a variety of “criticality” measures, based on the DC
power flow model, for the different components in the
system that aid in identifying these contingencies.
A different line of work for general flow networks is
based on a generalized network performance measure.
This measure is in turn used to rank the nodes and edges
of the network (see [21]); the downside of using these
ranking schemes on electric transmission systems is that
they are not validated using power flow models.
As for the literature in stochastic network interdiction
problems, the authors in [22] provide a thorough introduction to many variations of stochastic network interdiction.
In this context, we remark that the PNK is very different
from stochastic variants of the N-k interdiction problem
that compute an N-k contingency that maximizes the expected damage given a set of component-outage scenarios.
In contrast, the PNK uses the probability of occurrence
of an N-k contingency to weight the damage incurred by
that corresponding contingency; it seeks to identify an
N-k failure that maximizes the weighted damage. The
PNK generalizes the deterministic N-k interdiction problem, i.e., when the likelihood of occurrence of every N-k
contingency is the same, it exactly reduces to the deterministic problem. Stochastic network interdiction has also
been studied in the context of identifying locations for
installing smuggled nuclear material detectors [23]. The
authors in [23] formulate the problem of identifying locations for installing smuggled nuclear material detectors
as a two-stage stochastic mixed-integer program with recourse, and showed that the problem is NP-Hard. To the
best of our knowledge, the stochastic variants of the network interdiction problem, previously addressed in operations research literature, are only for general capacitated
flow networks [22, 24] and not for electric transmission
systems. Capacitated flow networks can be represented
by linear constraints and in some cases have a total unimodularity structure that can be exploited for developing
computationally efficient algorithms [25]. Furthermore,
using the probability of occurrence of an N-k failure as
weights to the damage incurred by that failure is a natural
generalization and extension of failure modeling in transmission systems. Transmission system operators identify

a N-k failure as a failure with exactly k components failing,
even in a probabilistic sense. To the best of our knowledge,
this is the first work in the literature to analyze a probabilistic generalization of the network interdiction problem
for electric transmission systems.

1.2.

Contributions

The PNK of this paper can be considered a generalization
of the work in [1] in two ways:
1. We consider a probabilistic outage model rather than
a deterministic model.
2. We use the AC power flow equations and their convex
relaxations instead of the DC power flow approximation.
The first generalization makes the objective function
concave, after suitable transformations, and the second
generalization introduces a set of convex-quadratic constraints instead of a set of linear constraints. The main
contributions of this paper are as follows.
1. We develop the first formulation of the PNK as a
mixed-integer, nonlinear program (MINLP) using the
AC power flow equations.
2. We develop one convex relaxation of the PNK that is
based on existing convex relaxations of the AC power
flow equations. We develop two approximations of
the PNK with linear constraints.
3. We develop a general cutting-plane algorithm for all
formulations of the PNK problem and prove its exactness for certain cases.
4. We perform extensive computational experimentation that demonstrates our approach is scalable to
networks with thousands of nodes.

2.

Nomenclature

This section presents nomenclature and terminology, including those that are well understood and used routinely
by the power systems community. An interested reader is
referred to [26] for a detailed description of the terminology.
Sets:
N - set of buses (nodes) in the network
N (i ) - set of nodes connected to bus i by an edge
3

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

E - set of from1 edges (lines) in the network
E r - set of to edges in the network
S - set of N-k contingency scenarios
Ẽs - set of damaged edges in N-k contingency scenario s
Es - set of operational from edges in the network during
N-k contingency scenario s
Esr - set of operational to edges in the network during N-k
contingency scenario s
Variables:
Vis = vis ∠θis - AC voltage at bus i during scenario s ∈ S
gs
gs
gs
Si = pi + iqi - AC power generated at bus i during
scenario s ∈ S
`is - percent AC power (load) shed at bus i during scenario
s∈S
Sijs = pijs + iqijs - AC power flow on line (i, j) during
scenario s ∈ S
xij - binary interdiction variable for line (i, j) ∈ E
zs - a continuous variable representing the minimum load
shed for scenario s ∈ S
Constants:
Sid = pid + iqid - AC power demand at bus i
Yij = gij + ibij - admittance of line (i, j)
(vil , viu ) - bounds for voltage magnitude at bus i
gl

gu

gl

gu

( pi , pi ) - bounds for active power generated at bus i
(qi , qi ) - bounds for reactive power generated at bus i
θ∆
ij - maximum phase angle difference across line (i, j )
tij - thermal limit of line (i, j)
Prij - i.i.d. probability of failure of line (i, j)
Prs = ∏(i,j)∈Ẽs Prij - probability of N-k contingency
scenario s
k - number of damaged/interdicted lines in each N-k
contingency scenario
Other notations:
Re(·) - real part of a complex number
Im(·) - imaginary part of a complex number
(·)∗ - conjugate of a complex number
|·| - magnitude of a complex number
∠ - angle of a complex number
Unless otherwise stated, all the values are in per-unit
(pu). A pu system is the expression of system quantities
as fractions of a defined base unit quantity; this is ex1 Although

a transmission line is modeled as an undirected edge, the
power flowing on the line is asymmetric; from and to edges are used to
model this asymmetric power flowing on a line.

4

tensively used in power system engineering to normalize
the magnitude of different system quantities and improve
numerical stability. We also note that the probability of
failure of a line (i, j) is assumed to be i.i.d. This is a reasonable assumption in the power system context because
although extreme events such as hurricanes, earthquakes
and geomagnetic disturbances have a lot of structure, that
structure is determined via simulation, and is not considered as an uncertainty in typical fragility assessments (see
for instance [2]). Once the extreme event is fixed, individual component failures are conditionally independent of
one another during the particular event.

3.

A bilevel MINLP formulation

In this section, we formulate the PNK as a bilevel MINLP
using the nomenclature introduced in the previous section.
Although the approach presented in this paper holds for
any type of component (generator, transformer, substation,
transmission line), for ease of exposition, we assume that
only the transmission lines in the system can be disrupted
or interdicted; hence, throughout the rest of the article, N
and k represent the total number of transmission lines and
the number of them that can be interdicted in the system,
respectively. The problem is formulated as a bilevel maxmin problem. From here on, we will refer to the max and
min problem as the outer and inner problem, respectively.
The formulation is as follows:

(Fac ) :

max
s∈S

Prs · zs

(1)

where zs is defined by the following optimization problem
for every scenario s ∈ S :
zs = min

∑ Re(Sid )`is

subject to:

(2)

i ∈N
vil 6 |Vis | 6 viu ∀i ∈ N ,
gl
gs
gu
pi 6 Re(Si ) 6 pi
∀i ∈ N ,
gl
gs
gu
qi 6 Im(Si ) 6 qi
∀i ∈ N ,
[
s
|Sij | 6 tij ∀(i, j) ∈ Es Esr ,

− θ∆ij 6 ∠Vis Vjs∗ 6 θ∆ij
0 6 `is 6 1
gs
Si

(3)
(4)
(5)
(6)

∀(i, j) ∈ Es ,

(7)

∀i ∈ N ,

− (1 − `is )Sid

=

(8)

∑S

(i,j)∈Es

Sijs = Yij∗ Vis Vis∗ − Yij∗ Vis Vjs∗

Esr

Sijs

∀i ∈ N , and

∀(i, j) ∈ Es

[

Esr .

(9)
(10)

In the above formulation, the outer problem is analogous to the attacker’s problem, where the attacker is

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

interested in computing the N-k scenario s that maximizes
the damage to the system (see Eq. (1)); here, the attacker
measures the damage in terms of the minimum load shed
incurred by the defender (the system operator) during
scenario s weighted by the probability Prs of scenario s.
This objective is shown to be an effective severity indicator
for an N-k scenario via extensive simulations in the dissertation by Nedic [27]. Computing the value of Prs for every
N-k scenario s, presents a mathematical challenge. In the
following section, we will show how the independence
of the line failure probabilities Prij can be leveraged to
address this mathematical challenge via the use of bilinear
terms which are convexified using logarithms.

The inner problem is the defender’s problem, which is
to operate the power network for contingency s in such a
way as to minimize the active load shed in the network.
Although the defender seeks to minimize the active load
shed in the system, constraint (9) ensures that a certain
quantity of reactive load is also shed in the system. This
is achieved in constraint (9) by forcing the load-shedding
factor `is at each bus i during scenario s ∈ S to act on
both the real and reactive loads. This type of load shedding is also referred to as constant power-factor shedding.
The constraints for the inner problem can be classified
into two categories, viz., the physical limits of the network (operational side constraints) and the AC power flow
equations that govern the physics of power flowing in the
network. Constraints (3) and (7) bound the voltage magnitude at each bus and the phase angle difference between
the buses connecting a line, respectively. Constraints (4)
and (5) limit the real and reactive power generated at a bus
during contingency s, respectively. Constraint (8) gives
the restrictions on the load-shedding variable `is . Eqs. (9)
and (10) together define how power flows in the network
and form the core building block of many power system
applications; they are also the source of nonconvexity in
the defender’s problem. The PNK is a difficult problem
(in fact, the inner problem is NP-hard [28, 29]). Hence,
in the next section, we present a convex relaxation of the
inner problem and a convexification procedure to handle
the nonlinear terms in the outer problem to convert the
MINLP to a mixed-integer convex program. We note that
the convex relaxations of the AC power flow equations
in (9), presented in the next section, are established and
broadly accepted in the power systems and optimization
literature (see [26] and references therein).

4.

Convex reformulation of the outer
problem

As mentioned in the previous section, the outer problem
has a bilinear term in the objective function. Also, the
value of Prs that defines the probability of an N-k scenario
s ∈ S depends on the k transmission lines that are being
interdicted in scenario s. Generally, it is computationally
inefficient to maintain a look-up table that calculates the
value of Prs for each scenario. Instead, we model the
relationship between Prs and lines interdicted in the N-k
scenario s explicitly by using binary interdiction decision
variables xij for each line (i, j). We then use logarithms to
address both these issues and convert the outer problem
to a maximization of a concave function. After taking
logarithms, the outer problem is as follows:
log Prs + log zs

max
s∈S

log Prs =

∑

subject to:

log Prij · xij ,

(11)

(i,j)∈E

∑

xij = k,

(12)

(i,j)∈E

xij ∈ {0, 1}

∀(i, j) ∈ E ,

∀s ∈ S .

zs > 0

(13)
(14)

The decision variables in this problem are xij , zs , and
Prs . Once the N-k scenario s ∈ S is fixed, the values of
all the decision variables are completely determined. In
particular, the value of zs is determined by the the inner
problem, given by (2) – (10). Hence, the maximization
is over the set of all scenarios s ∈ S . Constraints (11),
(12), and (13) define the value of Prs based on the values
taken by the interdiction variables xij corresponding to the
N-k scenario s ∈ S . The assumption that the line-failure
probabilities Prij are independent is used to express Prs
using Prij s in constraint (11). We now let p = log Prs and
the optimization problem is equivalently written as:
max
p=

p + log zs

∑

subject to:

log Prij · xij , (12), (13), and (14).

(15)

(i,j)∈E

The function p + log zs in the objective is jointly concave
when p, zs > 0. Hence, the outer problem is a convex
optimization problem, provided p, zs > 0.
In the next section, we present convex relaxations and
linear approximations for the inner problem in the PNK.
The inner problem, as presented in Section 3, is nonlinear
5

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

because of the AC power flow equations in (10). This
makes computing a global optimum, even for reasonably
sized test instances, computationally intractable [30]. As
a result, there has been recent work developing convex
relaxations in the literature. The relaxations include the
semidefinite programming [31], the second-order cone
(SOC) [32], the quadratic convex [33], and the momentbased relaxations [34]. Each of these relaxations has an
associated solution quality/computation time trade-off. In
this paper we use the SOC relaxation because it has a favorable balance of computation time and solution quality in
this application domain. We also present two approximations to the AC power flow equations: the linearized DC
power flow formulation and a traditional network flow
formulation. It is important to note that the algorithm
presented in this paper is generic and can be applied to
any convex relaxations of the AC power flow equations.

5. Convex relaxation and linear
approximations of the inner problem
5.1.

The SOC relaxation

First we present the SOC relaxation of the AC power flow
equations in (10) to obtain a SOC program for the inner
problem for a fixed N-k scenario. We begin by rewriting
the constraints for the inner problem in the following
equivalent lifted form:
Wijs = Vis Vjs∗

∀(i, j) ∈ Es ,

(vil )2 6 Wiis 6 (viu )2
tan(−θ∆
ij ) 6

Im(Wijs )
Re(Wijs )

Sijs = Yij∗ Wiis − Yij∗ Wijs

(16)

∀i ∈ N ,

(17)

6 tan(θ∆ij ) ∀(i, j) ∈ Es ,

(18)

∀(i, j) ∈ Es

[

Esr ,

(19)

and constraints (4), (5), (6), (8) and (9).
The introduction of these auxiliary variables ensures
that the nonlinearity in the inner problem is isolated in
constraint (16). Given the above set of lifted constraints
for the inner problem, the SOC relaxation is obtained by
replacing constraint (16) with |Wijs |2 6 Wiis Wjjs for every
(i, j) ∈ Es . For the sake of completeness, the full formulation for the PNK, with the convex reformulation of
the outer problem and the SOC relaxation of the inner
problem, is shown below:

subject to: p =

6

max

p + log zs


log Prij · xij , (12), and (13)

(i,j)∈E

where zs for every N-k scenario s ∈ S is defined by

∑ Re(Sid )`is

zs = min

i ∈N

subject to:

|Wijs |2 6 Wiis Wjjs

∀(i, j) ∈ Es ,

(20)

(17), (18), (19), (4), (5), (6), (8), and (9).

5.2.

The DC approximation

In this section, we present the linearized DC approximation of the AC power flow equations in (10); we note that
the DC approximation results in the inner problem being
a linear program for a fixed N-k scenario s ∈ S . To make
the presentation of the DC approximation simple, it is
useful to represent the AC power flow equations in their
polar form. Eq. (10), in polar coordinates, is given by the
following set of two equations for every (i, j) ∈ Ēs :


pijs = (vis )2 gij − vis vsj gij cos(θis − θ js ) + bij sin(θis − θ js ) ,


qijs = −(vis )2 bij − vis vsj gij sin(θis − θ js ) − bij cos(θis − θ js ) .
Three basic assumptions are used to derive the DC approximation of the AC power flow constraints in Cartesian
form. They are as follows:
(A1) |bij |  | gij | for every line (i, j).
(A2) The voltage magnitude at each node or bus is 1 pu,
i.e., vis = 1 pu for each bus i.
(A3) The voltage angle difference (θis − θ js ) across any line
(i, j) ∈ E is small enough that cos(θis − θ js ) ≈ 1 and
sin(θis − θ js ) ≈ (θis − θ js ).
Using these three assumptions on the AC power flow
equations yields
pijs = −bij (θis − θ js ) and qijs = 0

∀(i, j) ∈ Es

[

Esr .

(21)

Applying Eq. (21) in the inner problem for every scenario
s ∈ S and the convex reformulation of the outer problem
yields the following DC approximation of the PNK:

(Fdc ) :

max

subject to: p =

(Fsoc ) :

∑

p + log zs

∑

(i,j)∈E


log Prij · xij , (12), and (13)

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

where

∑

zs = min
gl

pi 6

| pijs |

i ∈N
gs
gu
pi 6 pi

subject to:

(22)

∀i ∈ N ,

(23)

∀(i, j) ∈ Es ,

6 tij

0 6 `is 6 1
gs
pi

pid `is

(24)

∀i ∈ N ,

− (1 − `is ) pid

(25)

∑S

=

(i,j)∈Es

pijs = −bij (θis − θ js )

Esr

pijs

∀(i, j) ∈ Es

∀i ∈ N , and
[

Esr .

(26)
(27)

In the next section, we present another linear approximation of the inner problem of the PNK: a network flow (NF)
approximation. We note that the NF approximation is
actually a relaxation of the DC approximation presented
in this section.

5.3.

NF approximation

In this section, we present a relaxation of the DC approximation using capacitated flows in networks. The NF
approximation for the PNK is obtained by dropping constraint (27) from the DC approximation presented in the
previous section. It is trivial to observe that for a given N-k
scenario, the inner problem without (27) is a maximum
flow problem. For the sake of completeness, the PNK with
the NF approximation is shown below:

(Fnf ) :

max

subject to: p =

p + log zs

∑


log Prij · xij , (12), and (13)

(i,j)∈E

where
zs = min

∑

pid `is subject to: (23), (24), (25), and (26).

i ∈N

6.

A cutting-plane algorithm

In this section, we present a generic cutting-plane algorithm that is common to the MINLP, the SOC relaxation,
and the DC and NF approximation for the PNK. The
cutting-plane algorithm uses the Stackelberg game structure that is inherent to all the formulations. The main
difficulty in developing any algorithm for the PNK arises
from the nonconvex max-min nature of the problem. A
number of techniques have been studied and analyzed in
the literature to convert bilevel max-min problems like the
PNK to a single mixed-integer program (see [5, 7]). These
techniques are known to have scaling issues for large-scale

test instances. To understand why scaling is a problem, it
is useful to think about an N-k contingency scenario s as
modifying the original set of lines E to be Es . Representing this behavior as constraints is key to converting the
bilevel program to a single-level mixed-integer program
and can be achieved by using disjunctions. Furthermore, if
the N-k scenario includes node or substation failures, this
further complicates the model with additional multilinear
terms [1]. Instead, we present a cutting-plane algorithm
that, instead of converting the bilevel PNK to a single-level
problem, generates cutting planes (using the solution of
the inner problem) that are added sequentially to the outer
problem.
The algorithm relies on constructing a sequence of piecewise linear functions that upper bound the total active load
shed zs , given by the solution to the inner problem or its
relaxation and approximations, for any s ∈ S . For any N-k
scenario ŝ ∈ S and the associated solution xij for (i, j) ∈ E ,
let zŝ denote the minimum load shed that is obtained by
solving the inner problem given by (2) – (10) or one of its
convex relaxation or approximations in Section 5. Then,
the algorithm computes coefficients αij (ŝ) for each line
(i, j) ∈ E such that

∑

zs 6 zŝ +

αij (ŝ) · xij

∀s ∈ S .

(28)

(i,j)∈E

The cut in (28) has been developed and found to be effective for the deterministic variant of formulation Fdc (see
[1]). In this section, we show that this cut generalizes to
the PNK and to different relaxations and approximations
of the AC power flow equations. The inequality in (28)
provides an upper bound for zs , the objective function of
the inner problem, for every feasible N-k scenario s ∈ S .
The linear cut in (28) is very general, and there are many
choices for the cut coefficients αij (ŝ). The key challenge
is to choose tight values for each coefficient that do not
remove any feasible N-k scenario. For the PNK, these
coefficients αij (ŝ) are computed using a combination of
the subproblem solution for scenario s and the physics of
power flow in the network. Using the inequality in (28),
the PNK is equivalently written as:

(Fcp )

max

zs 6 zŝ +

∑

p + log zs

subject to:

αij (ŝ) · xij

∀ŝ ∈ S ,

(i,j)∈E

p=

∑


log Prij · xij , (12), and (13).

(i,j)∈E

Pseudocode of the cutting-plane algorithm for the PNK
7

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

using formulation Fcp is given in Algorithm 1. The pseudocode assumes a black box for solving the inner problem
and computing the cut coefficients in (28).
Algorithm 1 Cutting-plane algorithm: pseudocode
Input: optimality tolerance, ε > 0
Output: s∗ ∈ S , an ε-optimal N-k scenario to the PNK
1: initial problem: Fcp without constraint (28)
2: f ∗ ← − ∞
. lower bound on the optimal obj. value
3: f¯ ← + ∞
. upper bound on the optimal obj. value
4: ŝ ← any initial N-k scenario
5: p (ŝ ) ← ∑(i,j)∈Ẽ log Prij
ŝ
6: solve subproblem for ŝ and let zŝ be the objective value
7: f (ŝ ) ← p (ŝ ) + log zŝ
8: if f (ŝ ) > f ∗ then f ∗ ← f (ŝ ) and s∗ ← ŝ
compute αij (ŝ) for every (i, j) ∈ E satisfying (28)
add zs 6 zŝ + ∑(i,j)∈E αij (ŝ) · xij to Fcp and resolve
11: update ŝ, and set f¯ using solution from Step 10
12: if f¯ − f ∗ 6 ε f ∗ then (s∗ , f ∗ ) is the ε-optimal solution
to the PNK, stop
9:

10:

13:

return to step: 6

6.1.

Computing cut coefficients

Here we present an algorithm to compute the cut coefficients αij (ŝ) for a given N-k scenario ŝ and for every
(i, j) ∈ E using the solution to the inner problem. First
we present an expression for the cut coefficients, and then
we present a proof of validity of the cut when the NF
approximation is used for the inner problem followed by
the necessary conditions under which the cut is still valid
for the DC approximation and the SOC relaxation of the
inner problem. Regardless of the formulation that is used
for the inner problem (Fsoc , Fdc , or Fnf ), the solution to
any of these problems for a fixed N-k scenario ŝ provides
S
a value of pijŝ for each line (i, j) ∈ Eŝ Eŝr . Each pijŝ repS
resents the active power flowing on line (i, j) ∈ Eŝ Eŝr
during the N-k scenario ŝ. Given these values of pijŝ , the
cut coefficients are defined by the following equation:

max{| pŝ |, | pŝ |} if (i, j) ∈ E \ Ẽ ,
ŝ
ij
ji
αijŝ =
(29)
0
otherwise.
Recall that Ẽŝ is the set of lines that is removed for the
N-k scenario ŝ. When the coefficient αijŝ provides an upper
bound on the total amount of active load that must be
8

shed as the line (i, j) is removed from the network during
an N-k scenario ŝ, the cut given by (28) is valid. The
following theorem proves that the cut in (28) with the cut
coefficients in (29) is valid for formulation Fnf , which uses
the NF approximation of the inner problem.
Theorem 1. For any N-k scenario ŝ ∈ S , the cut zs 6 zŝ +
∑(i,j)∈E αij (ŝ) · xij , with αij (ŝ) taking the values in Eq. (29), is
valid for formulation Fnf .
Proof. As discussed in Section 5.3, the inner problem in
formulation Fnf is a maximum flow problem because
minimizing the active load shed in the system is equivalent to maximizing the active load satisfied in the system.
Maximizing the amount of active load satisfied in the
system subject to the NF constraints in Fnf is simply the
maximum flow problem in a directed graph; pijŝ represents the amount of active power flowing through line
S
(i, j) ∈ Eŝ Eŝr . At the optimal solution of the inner problem of Fnf for a fixed s ∈ S , the maximum amount of
load that can be shed when line (i, j) is removed from the
network is given exactly by the expression in (29). This
follows from the max-flow min-cut theorem [25].
We note that a trivial value of αij (ŝ) = ∑i∈N pid for every
(i, j) ∈ E \ Ẽŝ (total active power demand in the system)
results in a valid cut for the linearized DC approximation
and the SOC relaxation in Sections 5.2 and 5.1, respectively.
But, this value for cut coefficients is very weak and leads
to poor convergence of the cutting-plane algorithm. This
motivates the use of heuristic choices for cut-coefficients
that have better convergence properties. As a result, the
value of the cut coefficients, as given by Eq. (29), is valid
for the formulation Fnf . For the DC approximation of
the PNK presented in Section 5.2, the cut given in (28)
using the coefficients in (29) is not necessarily valid and
can remove feasible solutions to the PNK. An interested
reader is referred to [35, 1] for a counterexample. Nevertheless, the following theorem mathematically formalizes
an empirical statement in [1] concerning the validity of
the cut in (28) using the coefficients in (29). This theorem
specifies a necessary condition under which the cut given
by (28) is valid for the formulations Fdc and Fsoc .
Theorem 2. For any N-k scenario ŝ ∈ S , the cut zs 6 zŝ +
∑(i,j)∈E αij (ŝ) · xij , with αij (ŝ) taking the values in Eq. (29),
is valid for the formulations Fdc and Fsoc if the following
condition is satisfied:

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

(B1) The removal of the set of lines in (i, j) ∈ Eŝ , each carrying
pijŝ and pŝji units of real power during contingency ŝ from i
to j and j to i, respectively, leads to an additional total real
power load shedding of at most ∑(i,j)∈Eŝ max{| pijŝ |, | pŝji |}.
Proof. The cut in (28) can be interpreted as follows: given
a N-k scenario ŝ ∈ S with an inner problem objective
value zŝ , removal of a set of lines (i, j) ∈ Eŝ would result in
an increase in load shed by a value of at most ∑(i,j)∈Eŝ αij .
Eq. (29) and the assumption (B1) together indicate that

∑

αij =

(i,j)∈Eŝ

∑

max{| pijŝ |, | pŝji |},

7.

Computational results

(i,j)∈Eŝ

completing the proof.
In general, assumption (B1) need not hold (see [35])
because of Braess’s paradox as exhibited in electric transmission systems. The validity of assumption (B1) has been
empirically observed for many available power system test
cases via extensive computational experiments (see [1]) for
the deterministic variant of the DC approximation. As for
the SOC relaxation of the PNK given by formulation Fsoc ,
we make assumption (B1); in fact, the sufficiency condition
for the validity of the cut given in (28) and (29) for models
Fsoc and Fdc is still an open question. This makes the
cutting-plane algorithm a heuristic for formulations Fsoc
and Fdc .

6.2.

makes the cutting-plane algorithm a complete enumeration algorithm. This itself ensures finite convergence
because there are a finite number of N-k scenarios. (28), if
valid (for Fnf or for Fdc and Fsoc under the necessary conditions given by Theorem 2), ensures that not all feasible
solutions are enumerated and accelerates the convergence.
If the cut is not valid, then using results in a suboptimal
solution does not affect the finite convergence property,
completing the proof.

Logical constraints

We now present an additional logical constraint that is
added to speed up the convergence of the cutting-plane
algorithm. Given a feasible N-k scenario ŝ ∈ S , the constraint is given by

∑

xij 6 k − 1.

(30)

(i,j)∈Ẽŝ

This constraint forces any new N-k scenario to differ from
ŝ in at least one component. This type of cut has been
observed to be effective in similar interdiction problems
[36, 1]. In the following theorem, we prove that the cuttingplane algorithm, along with the logical constraints in (30),
converges in a finite number of iterations.
Theorem 3. The cutting-plane algorithm presented in Algorithm 1, along with the logical constraints in (30), converges in
finite number of iterations.
Proof. We start by observing that adding just the logical
constraints in (30) (and neglecting the cuts in (28)) to Fcp

In this section, we demonstrate the computational effectiveness of the cutting-plane algorithm and heuristics for
computing an optimal/feasible solution to the linear approximations and relaxation of the PNK. We use five test
instances for all our simulations, classified into the following four categories:
(C1) Small test instances: IEEE 14-bus and IEEE single-area
RTS96 (with 24 buses) test systems.
(C2) Medium test instances: IEEE three-area RTS96 (with
73 buses) and IEEE 118-bus test systems.
(C3) Large test instances: WECC 240 test system. We note
that the WECC 240 is an aggregated model of a real
transmission system in the western United States.
(C4) Very large test instances: PEGASE 1354-bus and the
Polish 2383-bus Winter peak test systems.
All computational experiments were run on an Intel
Haswell 2.6 GHz, 62 GB, 20-core supercomputing machine at Los Alamos National Laboratory. The algorithms were implemented using PowerModels [37], an
open-source framework for implementing power flow formulations. In the interest of using state-of-the-art commercial solvers, which support only convex quadratically
constrained quadratic programs, a linear outer approximation of the outer problem objective function is used
in place of its concave counterpart. CPLEX 12.7 is used
for solving all the mixed-integer linear and convex optimization problems presented in this paper. The data on
the probability of line failures are available as a part of
reliability data for the IEEE single-area and three-area test
cases [38]. For the remaining three test cases, the probability of line failure for each line is generated using the
following procedure. We compute the maximum and minimum probability (pmax and pmin , respectively) among all
9

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

the failure probabilities for the IEEE single-area RTS96 test
case and generate the probability of line failures uniformly
in the interval [ pmin , pmax ] for the other test cases. This
data and the code are made available in [39].
In order to ensure that the solutions obtained using
the cutting plane algorithm on Fsoc , Fdc , and Fnf are
meaningful for the PNK, the nonlinear inner problem in
Fac is solved using the N-k scenario obtained from the
cutting-plane algorithm using Ipopt [40] as the local solver.
This process of obtaining a feasible solution to the PNK
using its MINLP formulation, Fac , is also called primal
recovery or feasible solution recovery.

# buses

k

cutting-plane
obj. value (MW)

complete enumeration
obj. value (MW)

14
14
14
24
24
24

2
3
4
2
3
4

23.99
11.51
7.47
28.75
15.52
20.48

23.99
11.51
7.47
28.75
15.52
20.48

Table 2: Comparison of the objective values obtained by using the
cutting-plane algorithm on Fdc (DC approximation) and the
optimal solution obtained via complete enumeration.

7.1. Comparison of solution with total enumeration

# buses

k

cutting-plane
obj. value (MW)

complete enumeration
obj. value (MW)

In this section, we concentrate on only the small instances
in category (C1) where an optimal solution to formulations Fsoc , Fdc , and Fnf can be obtained by complete
enumeration of all the N-k scenarios. Here, we compare
the solutions obtained by the three formulations with their
corresponding optimal solutions obtained via complete
enumeration of all possible feasible solutions. We restrict
our attention to values of k in {2, 3, 4}.

14
14
14
24
24
24

2
3
4
2
3
4

24.02
11.53
7.47
28.75
19.18
20.97

24.02
11.53
7.47
28.75
19.18
20.97

# buses

k

cutting-plane
obj. value (MW)

complete enumeration
obj. value (MW)

14
14
14
24
24
24

2
3
4
2
3
4

23.99
11.51
7.47
28.75
15.52
20.48

23.99
11.51
7.47
28.75
15.52
20.48

Table 1: Comparison of the objective values obtained by using the
cutting-plane algorithm on Fnf (NF approximation) and the
optimal solution obtained via complete enumeration.

It can be observed from Tables 2 and 3 that the cuttingplane algorithm is able to compute an optimal solution to
formulations Fdc and Fsoc , despite being a heuristic, for
the instances in category (C1) for values of k ∈ {2, 3, 4}.
The results in Table 1 are not surprising because the
cutting-plane algorithm is indeed an exact algorithm for
the PNK with the NF approximation. Figs. 1 and 2 present
the number of iterations of the cutting-plane algorithm for
the SOC relaxation, the DC approximation, and the NF
approximation in the IEEE 14-bus system and the IEEE
10

Table 3: Comparison of the objective values obtained by using the
cutting-plane algorithm on Fsoc (SOC relaxation) and the
optimal solution obtained via complete enumeration.

single-area RTS96 system with 24 buses, respectively. A
consistent trend that emerges from Figs. 1 and 2 is that
the number of iterations of the cutting-plane algorithm is
highest for the NF approximation and lowest for the SOC
relaxation of the PNK, except for in the IEEE 14-bus system when k = 2. The computation time is not reported for
these small instances because the cutting-plane algorithm
converged to an optimal solution within 5 seconds for all
the runs.

7.2. Performance of the cutting-plane algorithm
on categories (C2) and (C3)
In this section, we present the computational results for the
test cases in categories (C2) and (C3). Unlike the instances
in category (C1), computing the optimal solution using
complete enumeration is time consuming for the instances
in categories (C2) and (C3), so we present the results of the
cutting-plane algorithm applied to formulations Fsoc , Fdc ,
and Fnf for instances in category (C1) only. A tolerance of
ε = 1% is set for every run of the cutting-plane algorithm.

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

10

# of iterations

Fnf

Fdc

Fsoc

8
6
4
2
k=2

k=3

k=4

k=5

Figure 1: Number of iterations of the cutting-plane algorithm applied
to formulations Fsoc , Fdc , and Fnf in the IEEE 14-bus
system.

# of iterations

Fnf

Fdc

Fsoc

100

50

0
k=2

k=3

timed out, the gap in % is reported instead of the objective
value; the gaps are enclosed in parentheses to differentiate
them from the ε-optimal solutions). We note that this
column reports the cost of the heuristic solution, obtained
using the cutting-plane algorithm, for formulations Fsoc
and Fdc .
ac: the objective value of the feasible solution to the PNK
obtained via primal recovery. This objective value is not
reported for instances that timed out.

k=4

k=5

Figure 2: Number of iterations of the cutting-plane algorithm applied
to formulations Fsoc , Fdc , and Fnf in the IEEE single-area
RTS96 system with 24 buses.

A computation time limit of 24 hours was imposed on
all the runs of the algorithm. Table 4 summarizes the
computational behavior of the cutting-plane algorithm
and the heuristic for the instances in category (C2) and
(C3), respectively. The column headings used in the table
are defined as follows:
# buses: number of buses (nodes) in the system.
k: number of interdicted lines (edges).
Fi : cutting-plane algorithm applied to formulation Fi ,
i ∈ {nf, dc, soc}.
time: time in seconds for the cutting-plane algorithm to
converge to ε-optimal solution; if the algorithm timed out
at 24 hours, then “TO” (timed-out) status is reported.
iter: number of iterations of the cutting-plane algorithm.
sol: ε-optimal objective value (in MW) if the cutting-plane
algorithm terminated within 24 hours or best feasible
solution obtained at the end of 24 hours (if algorithm

We observe from Table 4 that the SOC relaxation of
the PNK given by formulation Fsoc performs the best
from a computational perspective. The cutting-plane algorithm consistently takes the least number of iterations
when applied to the SOC relaxation. The cutting-plane
algorithm performs the worst when applied to the NF
approximation of the PNK. This can be attributed to the
fact that the NF approximation has a larger feasible region of the PNK, whereas the SOC relaxation is a smaller
representation of the feasible set of the PNK. Despite the
fact that the cut in (28) is not necessarily valid for the SOC
relaxation and the DC approximation, the primal recovery procedure to obtain a feasible solution to the MINLP
formulation for the PNK is observed to provide similar
objective values to both these formulations. Hence, it is
hard to distinguish between or identify any trends in the
actual solutions produced by the cutting-plane algorithm
when applied to formulations Fdc and Fsoc . Fig. 3 shows
the Hamming distance between the solutions obtained by
the cutting-plane algorithm on the three formulations for
the WECC 240 test system; it can be observed from the
figure that the Hamming distance between the interdiction
plans obtained using formulations Fsoc and Fdc is zero
for k ∈ {2, 3, 4, 5}.
There are outliers to the fact that the cutting-plane algorithm on Fdc and Fsoc results in the same choice of
interdiction plans; for instance consider the instance IEEE
three-area RTS96 with 73 buses: when k = 3, the Hamming distance between the interdiction plans obtained
from the two formulations is 2 and the difference between
the objective values obtained using the primal recovery
procedure is 15%. In fact, there are many other problems
in power systems [41, 42] where the SOC relaxation and
the DC approximation of the power flow constraints produce contrasting results. We delegate to future work the
need for further study to determine if the sameness of
the solutions between the two formulations holds more
11

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

Hamming distance

2
1.5
1
0.5
0
k=2
NF vs. DC

k=3

k=4

NF vs. SOC

k=5
DC vs. SOC

from category (C3).
For ease of exposition, we will denote the line-failure
probability of line (i, j) ∈ E generated using the reliability
data by Prijr [38]. We perform two sets of computational
experiments, referred to as (D1) and (D2). The first set
of experiments is aimed at examining how an increase
in line-failure probability in a certain geographical area
affects the objective value obtained using the cutting-plane
algorithm. This experiment is intended to simulate a
severe event that could potentially increase the line-failure
probabilities in a geographical region. To that end, we

Figure 3: Hamming distance between the interdiction plans obtained
by the cutting-plane algorithm on the three formulations for
the WECC 240 test system.

generally for the PNK.
However, taking into account the computation time reveals the SOC relaxation in formulation Fsoc to be the
clear winner between the three formulations, from both
a solution quality and computation time standpoint, for
the PNK. A likely reason for the relative computational
efficiency of the cutting-plane algorithm applied on Fsoc
is that the nonlinear line losses captured by the SOC relaxation has the property of removing symmetric solutions.
In contrast, the DC approximation can potentially allow
solutions (N-k failures) which produce the same amount
of load shed. The number of iterations of the cutting-plane
algorithm on Fsoc and Fdc is also consistent with this hypothesis. This type of behaviour was also observed in the
literature on similar problems such as Optimal Transmission Switching (see [33, 41]). However, it is uncommon.
Hence, we will use the SOC relaxation of the PNK for the
remaining computational experiments.

Figure 4: Buses that are assumed to be affected by a severe event on
the coast of California are represented using squares and the
other buses in the WECC 240 test system are represented
using circles.

first choose a set of 117 transmission lines, L ⊂ E , in
the state of California using the geolocated WECC 240
test instance. The set of buses on which these lines are
incident is shown in Fig. 4; this subset of buses is assumed
to be affected by a severe event on the coast of California.
The line-failure probability for each line in E is generated
using the following equation:

r
Pr r + n · 1−Prij if (i, j) ∈ L,
ij
5
Prij =
(31)
Pr r
otherwise.
ij

7.3. Effect of the failure probability distributions
For all the computational experiments thus far we either
used the line-failure probabilities from publicly available
data (reliability data) or generated data from a uniform
distribution in the range defined by the reliability data.
Now, we study the effect of a change in distribution of
line-failure probabilities on the solution obtained by the
cutting-plane algorithm. For the purposes of this study,
we restrict our attention to the SOC relaxation of the PNK
given by the formulation Fsoc and the instance WECC 240
12

We use n = {1, 2, 3} to generate three sets of line-failure
probabilities. We also note that the higher the value of
n, the higher the probability of line failure for the lines
in set L. Hence, n = 3 can correspond to a more severe
event than n ∈ {1, 2}. Also, when n = 0, the line-failure
probabilities reduce to the probabilities that were obtained
using the reliability data, studied in the previous section.
Fig. 5 shows the objective value obtained by the cuttingplane algorithm for formulation Fsoc for the WECC 240
test system for k ∈ {2, . . . , 15} using the line-failure probabilities generated for the computational experiments in

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

# buses

Fnf

k

Fdc

Fsoc

time

iter

sol

ac

time

iter

sol

ac

time

iter

sol

ac

73
73
73
73
73
73
73
73
73

2
3
4
5
6
7
8
9
10

235.06
45804.10
1715.02
18438.44
TO
TO
TO
TO
TO

242
1256
336
641
911
794
716
687
673

28.75
15.15
20.48
11.06
(26.97)
(31.19)
(26.17)
(23.68)
(22.23)

28.75
15.15
20.96
11.32
—–
—–
—–
—–
—–

31.09
580.29
16.66
61.12
118.49
335.76
511.23
640.90
435.06

107
382
47
99
127
167
202
194
217

28.75
15.52
20.48
11.06
5.97
3.22
1.71
0.91
0.47

28.75
15.52
20.96
11.32
6.11
3.30
1.73
0.91
0.47

29.61
66.27
6.30
9.69
19.18
29.87
40.43
47.84
57.70

117
144
22
35
60
77
86
93
111

28.75
19.19
20.96
11.32
6.11
3.30
1.74
0.91
0.47

28.75
19.18
20.96
11.32
6.11
3.30
1.74
0.91
0.47

118
118
118
118
118
118
118
118
118

2
3
4
5
6
7
8
9
10

208.81
13565.87
TO
TO
TO
TO
TO
TO
TO

260
1199
1809
1497
1293
1095
951
866
840

21.22
11.40
(16.99)
(16.61)
(26.97)
(31.19)
(26.17)
(23.68)
(22.23)

21.11
11.40
—–
—–
—–
—–
—–
—–
—–

9.90
24.83
61.30
348.21
1205.63
3161.65
8125.70
15615.67
20622.50

33
83
124
245
348
458
588
677
708

21.12
11.40
6.83
3.67
2.02
1.09
0.58
0.30
0.16

21.11
11.40
6.83
4.57
2.02
1.09
0.58
0.30
0.16

8.72
26.87
43.81
115.36
319.29
508.70
1543.57
1777.69
2867.04

28
75
101
167
252
297
447
444
507

21.12
11.40
6.92
3.73
2.02
1.09
0.58
0.31
0.16

21.11
11.40
8.52
4.60
2.02
1.09
0.58
0.36
0.19

240
240
240
240
240
240
240
240
240

2
3
4
5
6
7
8
9
10

9.60
14.13
15.71
39.54
127.70
526.72
1293.63
2362.78
3939.01

21
41
46
88
169
312
409
473
523

2227.27
1329.31
778.54
430.66
235.24
127.03
69.04
37.28
20.09

2287.13
1418.38
831.09
451.17
242.83
131.13
71.04
38.36
20.70

1.66
2.73
2.87
3.09
4.65
4.88
8.94
11.08
8.02

6
13
13
15
22
23
36
44
43

2331.26
1381.54
796.09
437.94
236.49
129.09
70.05
37.88
20.30

2482.68
1471.39
831.09
458.66
247.68
134.86
72.19
38.98
21.17

7.35
7.70
10.17
11.04
18.16
19.83
20.50
20.46
17.62

8
8
12
13
23
25
27
27
23

2393.18
1422.93
808.91
447.24
241.77
130.55
70.82
38.32
20.65

2482.68
1471.39
831.09
458.66
248.79
134.35
72.83
39.39
21.24

Table 4: Computational results for instances in categories (C2) and (C3).

13

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

obj. value (MW)

3,000

3,000

=0
=1
=2
=3

obj. value (MW)

n
n
n
n

4,000

2,000
1,000

2,000

1,000

0

0
2

4

6

8

10

12

14

2

16

k

Figure 5: Plot of the objective value vs. k using the probability data
using Eq. (31). n = 0 shows the plot for the case when
line-failure probabilities are obtained using the reliability
data, Prijr .

set (D1). It is clear from the figure that the objective value
increases with increasing line-failure probabilities; this is
expected given the objective function of the PNK.
The second set of experiments, (D2), is aimed at examining the effect of the line-failure probability distribution on
the solution obtained by the cutting-plane algorithm. To
make the comparison fair, we assume a probability budget,
which is defined as the sum of the line-failure probabilities
of all the lines using the probabilities generated using the
reliability data, Prijr . The probability budget is defined as
B := ∑(i,j)∈E Prijr . We will refer to the line-failure probabilities Prijr as “rel”, where rel is an abbreviation of reliability.
We consider three cases for the line-failure probabilities:
(1) we assume each Prij is exactly 0.5 (this is equivalent
to the deterministic variant of the PNK and we will refer
to this case as “det”), (2) we assume Prij ∼ U (0, 1) (a
uniform distribution in the interval (0, 1)), and (3) we assume Prij ∼ texp(1), where “texp(1)” denotes a truncated
exponential distribution in (0, 1) with rate 1. Once we generate the line-failure probabilities for all three cases, these
values are normalized to sum to the probability budget B .
Fig. 6 shows the objective value obtained by the cuttingplane algorithm decreasing for increasing values of k; this
is because the larger the value of k, the lower the probability of an N-k failure (this observation is also valid for
Fig. 5). The sum of the individual line-failure probabilities
is normalized to a constant value for all four cases; this
can be interpreted as moving the line-failure probabilities
from one line to another. This change does not seem to
affect the value of k at which the objective becomes comparatively small. This leads to a hypothesis that for a given
14

rel
det
U (0, 1)
texp(1)

4

6

8

10

12

14

16

k

Figure 6: Plot of the objective value vs. k for different distributions of
line-failure probability.

system and for a constant value of probability budget B ,
there exists a value of k greater than which the likelihood
of the system sustaining damage is low. Also, the objective value for a fixed B is observed to be greatest when
Prij ∼ texp(1) and least for the “det” case, when all the
probabilities take a value of 0.5 before normalization. The
fact that this objective value is lowest in the “det” case and
increases significantly in the “exp” case suggests that the
interdiction plan obtained by solving the PNK is sensitive
to the nature of the failure distribution. This observation
is also corroborated by Fig. 7. Fig. 7 shows the regions in
n=3
rel
det

Figure 7: Regions that are affected for the different experiments on the
WECC 240 instance.

which the lines are interdicted for the case when a severe
event occurs on the coast of California, with line-failure
probabilities generated using Eq. (31) with the value of n
set to 3, the “det” case, and the “rel” case.

7.4.

Scalability of the cutting-plane algorithm

In this section, we demonstrate the scalability of the
cutting-plane algorithm on the very large test instances

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

from category (C4). Given the computational and solution
quality superiority of Fsoc on instances in categories (C2)
and (C3), for the purposes of this study, we restrict our
attention to Fsoc . Furthermore, we increase the tolerance
from ε = 1% to ε = 5% and set a computation time limit of
24 hours for every run of the cutting-plane algorithm. An
ε value of 5% was chosen for this set of experiments as this
value of ε provided a reasonable trade-off between computation time and the solution quality. Table 5 summarizes
the computational behavior of the cutting-plane algorithm
applied to the formulation Fsoc for the “PEGASE” and
“Polish” test instances. The additional column headings,
used in the Table 5 are defined as follows:

natural or manmade severe events. Future work should
focus on modeling the dependence on component-failure
probabilities, dependence on severe event probabilities, as
well an in depth comparison of solutions obtained using
a deterministic N-k model. Further work should also focus on stochastic variants of the N-k interdiction problem
on power systems and comparing the solutions with the
PNK. Finally, including models of cascading failures in
the context on the N-k problems remains a computational
challenge and is also an important avenue of future work.

feas: objective value (MW) of the heuristic solution provided by the cutting-plane algorithm for ε = 5%.

[1] Javier Salmeron, Kevin Wood, and Ross Baldick.
Worst-case interdiction analysis of large-scale electric power grids. IEEE Transactions on Power Systems,
24(1):96–104, 2009.

( f¯− f ∗ )

gap: relative gap in % i.e., f ∗ · 100% after the termination of the cutting-plane algorithm.
The results in Table 5 indicate that the cutting-plane
algorithm can scale to very large test instances of the order
of 2000 buses. Even for the Polish instances that timed
out in the k = 6, . . . , 10 cases, the relative gap is below 6%
further corroborating the effectiveness of the algorithm.

8.

Conclusion

In this paper, we have presented an MINLP formulation, a
convex (SOC) relaxation, and linear approximations to the
probabilistic interdiction problem in power transmission
networks. To the best of our knowledge, this is the first
study of the N-k problem on power transmission systems
that incorporates the convex relaxation of the AC power
flow constraints. In summary, we have shown the following: (1) the computational effectiveness of the generic
cutting-plane algorithm based on work in the literature
for the deterministic variant of the PNK with the DC approximation [1] extends to the probabilistic generalization
and also to convex relaxations of the AC power flow equations; (2) the nonlinear convex relaxation (SOC) of the AC
power flow equations accelerates the convergence of the
cutting-plane algorithm in comparison to linear DC/NF
approximations; (3) the interdiction plan obtained using
the cutting-plane algorithm is sensitive to changes in the
component-failure probabilities, and this sensitivity can be
used to compute interdiction plans in smaller subnetworks
of the original system by spatially biasing the probabilities
in the regions of interest; and (4) the PNK is very effective at identifying high-risk failure plans in the case of

References

[2] Hazus, Dept. of Homeland Security (DHS). (Accessed: December, 2017).
[3] Gerald Brown, Matthew Carlyle, Javier Salmerón,
and Kevin Wood. Defending critical infrastructure.
Interfaces, 36(6):530–544, 2006.
[4] Javier Salmeron, Kevin Wood, and Ross Baldick. Analysis of electric grid security under terrorist threat.
IEEE Transactions on Power Systems, 19(2):905–912,
2004.
[5] Rogelio E Alvarez. Interdicting electrical power grids.
PhD thesis, Naval Postgraduate School, Monterey,
California, 2004.
[6] PowerWorld Corporation, Champaign, IL. PowerWorld Simulator.
[7] Alexis L Motto, José M Arroyo, and Francisco D
Galiana. A mixed-integer LP procedure for the analysis of electric grid security under disruptive threat.
IEEE Transactions on Power Systems, 20(3):1357–1365,
2005.
[8] José M Arroyo and Francisco D Galiana. On the
solution of the bilevel programming formulation of
the terrorist threat problem. IEEE Transactions on
Power Systems, 20(2):789–797, 2005.
[9] Daniel Bienstock and Abhinav Verma. The n-k problem in power grids: New models, formulations, and
numerical experiments. SIAM Journal on Optimization,
20(5):2352–2380, 2010.
15

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

k

2
3
4
5
6
7
8
9
10

PEGASE 1354-bus test system

Polish 2383-bus test system

time

iter

feas

gap

ac

time

iter

feas

gap

ac

133.25
159.27
398.36
954.42
2861.25
5469.59
28372.89
28401.17
29980.16

57
66
125
215
352
425
721
624
545

255.22
172.07
102.94
60.27
32.55
18.19
10.32
5.43
3.09

0.12
0.53
2.10
0.02
5.00
5.00
2.59
4.96
4.06

255.23
173.74
103.82
60.73
32.79
19.46
10.57
5.75
3.44

1046.78
2937.39
7520.70
17126.15
TO
TO
TO
TO
TO

64
164
335
557
1106
1015
913
832
750

100.72
58.88
34.75
18.76
9.95
5.87
3.10
1.70
0.89

0.83
1.31
0.64
5.00
5.98
5.02
6.75
6.77
7.39

109.77
63.02
36.85
19.90
10.74
7.23
3.82
2.08
1.10

Table 5: Computational results for instances in category (C4).

[10] Andrés Delgadillo, José Manuel Arroyo, and Natalia
Alguacil. Analysis of electric grid interdiction with
line switching. IEEE Transactions on Power Systems,
25(2):633–641, 2010.
[11] Taedong Kim, Stephen J Wright, Daniel Bienstock,
and Sean Harnett. Analyzing vulnerability of power
systems with continuous optimization formulations.
IEEE Transactions on Network Science and Engineering,
3(3):132–146, 2016.
[12] Vaibhav Donde, Vanessa López, Bernard Lesieutre,
Ali Pinar, Chao Yang, and Juan Meza. Severe multiple
contingency screening in electric power systems. IEEE
Transactions on Power Systems, 23(2):406–417, 2008.
[13] Ali Pinar, Juan Meza, Vaibhav Donde, and Bernard
Lesieutre. Optimization strategies for the vulnerability analysis of the electric power grid. SIAM Journal
on Optimization, 20(4):1786–1810, 2010.
[14] José Manuel Arroyo. Bilevel programming applied
to power system vulnerability analysis under multiple contingencies. IET Generation, Transmission &
Distribution, 4(2):178–190, 2010.
[15] Vicki M Bier, Eli R Gratz, Naraphorn J Haphuriwat,
Wairimu Magua, and Kevin R Wierzbicki. Methodology for identifying near-optimal interdiction strategies for a power transmission system. Reliability Engineering & System Safety, 92(9):1155–1161, 2007.
[16] José M Arroyo and Francisco J Fernández. Application of a genetic algorithm to n-k power system
security assessment. International Journal of Electrical
Power & Energy Systems, 49:114–121, 2013.
16

[17] Mark K Enns, John J Quada, and Bert Sackett. Fast linear contingency analysis. IEEE Transactions on Power
Apparatus and Systems, PAS-101(4):783–791, 1982.
[18] Margaret J Eppstein and Paul DH Hines. A “random
chemistry” algorithm for identifying collections of
multiple contingencies that initiate cascading failure.
IEEE Transactions on Power Systems, 27(3):1698–1705,
2012.
[19] C Matthew Davis and Thomas J Overbye. Multiple
element contingency screening. IEEE Transactions on
Power Systems, 26(3):1294–1301, 2011.
[20] Saleh Soltan and Gil Zussman. Quantifying the effect
of k-line failures in power grids. In Power and Energy
Society General Meeting (PESGM), 2016, pages 1–5.
IEEE, 2016.
[21] Qiang Qiang and Anna Nagurney. A unified network
performance measure with importance identification
and the ranking of network components. Optimization
Letters, 2(1):127–142, 2008.
[22] Kelly J Cormican, David P Morton, and R Kevin
Wood. Stochastic network interdiction. Operations
Research, 46(2):184–197, 1998.
[23] Feng Pan, William S. Charlton, and David P. Morton.
A stochastic program for interdicting smuggled nuclear
material, pages 1–19. Springer US, 2003.
[24] Udom Janjarassuk and Jeff Linderoth. Reformulation
and sampling to solve a stochastic network interdiction problem. Networks, 52(3):120–132, 2008.

Probabilistic N-k Failure-Identification for Power Systems • Sundar et.al • October 2017 • Vol. X, No. Y

[25] Thomas H Cormen. Introduction to Algorithms. MIT
Press, 2009.

open-source framework for exploring power flow formulations. 2017.

[26] Carleton Coffrin, Hassan L Hijazi, and Pascal
Van Hentenryck. The QC relaxation: A theoretical and computational study on optimal power flow.
IEEE Transactions on Power Systems, 31(4):3008–3018,
2016.

[38] Cliff Grigg, Peter Wong, Paul Albrecht, Ron Allan,
Murty Bhavaraju, Roy Billinton, Quan Chen, Clement
Fong, Suheil Haddad, Sastry Kuruganty, et al. The
IEEE reliability test system-1996. a report prepared
by the reliability test system task force of the application of probability methods subcommittee. IEEE
Transactions on Power Systems, 14(3):1010–1020, 1999.

[27] Dusko P Nedic. Simulation of large system disturbances.
PhD thesis, University of Manchester, 2003.
[28] Abhinav Verma. Power grid security analysis: An optimization approach. PhD thesis, Columbia University,
2009.
[29] Karsten Lehmann, Alban Grastien, and Pascal
Van Hentenryck. AC-feasibility on tree networks is
np-hard. IEEE Transactions on Power Systems, 31(1):798–
801, 2016.
[30] Chen Chen, Alper Atamtürk, and Shmuel S Oren.
Bound tightening for the alternating current optimal
power flow problem. IEEE Transactions on Power Systems, 31(5):3729–3736, 2016.
[31] Xiaoqing Bai, Hua Wei, Katsuki Fujisawa, and Yong
Wang. Semidefinite programming for optimal power
flow problems. International Journal of Electrical Power
& Energy Systems, 30(6):383–392, 2008.
[32] Rabih A Jabr. Radial distribution load flow using
conic programming. IEEE Transactions on Power Systems, 21(3):1458–1459, 2006.

[39] Kaarthik Sundar. GitHub repository for the probabilistic N-k failure-identification problem for power
systems.
https://github.com/kaarthiksundar/
n-k, April 2017.
[40] Andreas Wächter and Lorenz T Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1):25–57, 2006.
[41] Carleton Coffrin, Hassan L Hijazi, Karsten Lehmann,
and Pascal Van Hentenryck. Primal and dual bounds
for optimal transmission switching. In Power Systems
Computation Conference (PSCC), 2014, pages 1–8. IEEE,
2014.
[42] Harsha Nagarajan, Russell Bent, Pascal Van Hentenryck, Scott Backhaus, and Emre Yamangil. Resilient
transmission grid design: AC relaxation vs. DC approximation. 2017.

[33] Hassan Hijazi, Carleton Coffrin, and Pascal Van Hentenryck. Convex quadratic relaxations for mixedinteger nonlinear programs in power systems. Mathematical Programming Computation, 9(3):321–367, 2017.
[34] Daniel K Molzahn and Ian A Hiskens. Momentbased relaxation of the optimal power flow problem.
In Power Systems Computation Conference (PSCC), 2014,
pages 1–7. IEEE, 2014.
[35] Ross Baldick. Course notes for ee394v – restructured
electricity markets: Locational marginal pricing, 2006.
[36] Eitan Israeli and R Kevin Wood. Shortest-path network interdiction. Networks, 40(2):97–111, 2002.
[37] Carleton Coffrin, Russell Bent, Kaarthik Sundar,
Yeesian Ng, and Miles Lubin. Powermodels.jl: An
17

