1

Management of Cascading Outage Risk Based on
Risk Gradient and Markovian Tree Search

arXiv:1704.02556v2 [cs.SY] 6 Nov 2017

Rui Yao, Member, IEEE, Kai Sun, Senior Member, IEEE, Feng Liu, Member, IEEE, Shengwei Mei, Fellow, IEEE

Abstract—Since cascading outages are major threats to power
systems, it is important to reduce the risk of potential cascading
outages. In this paper, a risk management method of cascading
outages based on Markovian tree search is proposed. With the
tree expansion on the cascading outage risk, risk gradient is
computed efficiently by a forward-backward tree search scheme
with good convergence, and it is then employed in an optimization
model to minimize control cost while effectively reducing the
cascading outage risk. To overcome the limitation with linearization in computing risk gradient, an iterative risk management
(IRM) approach is further developed. Tests on the RTS-96 3area system verify the accuracy of the computed risk gradient
and its effectiveness for risk reduction. Time performance of the
proposed IRM approach is tested on the RTS-96 system, a 410bus US-Canada northeast system and a 1354-bus Mid-European
system, and demonstrates its potentials for decision support on
practical power systems online or on hourly basis.
Index Terms—cascading outage, tree search, risk management,
risk gradient, control cost, trade-off, multi-objective optimization

I. I NTRODUCTION
ASCADING outages are major threats to power system
operations. To study and prevent cascading outages,
people have proposed various methods on the simulation, risk
assessment and risk management of cascading outages [1]. The
mitigation methods [2]–[4] based on power-law analysis can
provide rough mitigation strategies in the planning horizon, but
lack accurate tactics for online operations. In grid operations,
an extensively used approach is to check deterministic N-1 or
N-k criteria to guarantee no limit being violated under given
contingencies. Another methodology is robust optimization
[5], which ensures the strategy to cover all foreseen high-risk
scenarios. However, the prevention of all potential cascading
outages could be too expensive to achieve for both planning
and operations of the system. Also it is practically infeasible
to completely wipe out the risk of any cascading outages [6].
Although the deterministic methods are more computationally
efficient, they cannot distinguish the actual probabilities and
consequences of different events, so the control cost could be
unnecessarily spent on preventing some low-risk events [7].
Therefore, it is more desirable to introduce the index of risk
(i.e. the expectation of consequences) to provide more effective
guidance for decision-making in power system operations.

C

This work was supported by the CURENT Engineering Research Center.
R. Yao and K. Sun are with the Department of EECS, University
of Tennessee, Knoxville, TN 37996, USA (emails: yaorui.thu@gmail.com,
kaisun@utk.edu).
F. Liu and S. Mei are with the State Key Laboratory of Power Systems,
Department of Electrical Engineering, Tsinghua University, Beijing 100084,
China. (email: lfeng@mail.tsinghua.edu.cn)

Several risk-based approaches have been proposed in the
recent years. Some methods apply risk indices based on
pre-defined severity functions to optimization models [8]–
[11]. The risk indices in these methods quantify the expected
violation of operation constraints as indirect representations
of the risk faced by the end users [7]. While some methods
establish the risk management under a given contingency
set as a chance-constrained optimization problem [12], [13].
Also some methods incorporate reliability indices into the
economic dispatch or unit commitment models [14]–[16],
and realizes coordination between the economic profit and
operation risk [17]. However, these methods did not consider
the dependency among outages and hence may underestimate
the risk when assessing cascading outages. Moreover, many
efficient reliability assessment and risk management methods
for independent contingencies are no longer effective when
dealing with dependent outages. In summary, there is not yet
a method designed for the risk management against dependent
cascading outages. Such risk management should be able to
provide direct risk measurement in the amount of load shed,
energy shortage or economic loss. Moreover, to avoid omitting
unaware risky scenarios, the effective risk management should
not be confined in a pre-set contingency set [7]. It is rather,
desirable to efficiently search the most risky set of cascading
outage patterns for risk management.
To realize effective and efficient risk management, all the
following features are desirable:
1) Reasonable simulation and efficient risk assessment of
multi-timescale dependent cascading outages: firstly, a reasonable modeling and simulation method should reflect the primary characteristics of cascading outages, e.g. the dependency
among outages, timescales of related processes, etc. [6], [18];
then, from numerous simulated outages, a risk assessment
method needs to identify the most risky cascade paths; finally,
the assessed risk is used in a risk management method to
reduce the risks of cascading outages.
2) Risk metrics that directly evaluate the risk faced by the
end users for effective risk management, such as load or energy
loss, and economic loss.
3) A risk management formulation allowing direct adjustment on reduction of risk, e.g. the amount of desired decrease
in the load loss, energy loss or economic loss, which will
facilitate to determining the trade-off between risk reduction
and control cost.
4) An efficient risk management algorithm to derive strategies for risk reduction on cascading outages.
Note that failing to realize any of the above desirable features would downgrade the effectiveness and efficiency in risk

2

management of cascading outages. This paper proposes a risk
management method having all the above features. To meet
the feature 1, the quasi-dynamic simulation method in [19]
and the Markovian tree (MT) based risk assessment method in
[20] are applied in this paper. The former can better reflect the
multi-timescale nature and dependencies in cascading outages,
and latter can efficiently identify the most risky cascading
outage patterns and provide a variety of practical risk metrics.
To meet the features 2 and 3, this paper introduces risk
gradient and risk constraint, which allow direct adjustment on
how much risk to be reduced on load loss, energy shortage
and economic loss. Moreover, to overcome the limitation of
linearization in the risk gradient computation, an iterative risk
management (IRM) approach is proposed for more effective
reduction of risk to meet feature 4. Although risk gradient
can be calculated in multiple ways, e.g. the cumbersome
perturbation method, this paper proposes an efficient forwardbackward algorithm to calculate risk gradient using the MT in
order for seamless integration into the MT search algorithm
for risk assessment. Thus, risk assessment and risk gradient
computation are performed at the same time to improve the
overall efficiency of the proposed approach.
The rest of the paper is organized as follows. Section
II first briefly introduces our previous work selected as a
foundation of the proposed approach (Section II.A). With the
MT expansion of risk, the generic model of risk management is
proposed based on the concept of risk gradient (Section II.B).
Section III is the analytical derivation of the risk gradient and
the efficient algorithm for calculating the risk gradient. The
risk gradient is derived in Section III.A-C. The computation of
risk gradient is realized as an efficient forward-backward treesearch algorithm, which is integrated with the MT-based risk
assessment algorithm (III.D). After deriving the risk gradient,
the Section IV implements the risk management based on
the generic form introduced in Section II: Section IV.A proposes the full optimization model of risk management (RM);
Section IV.B further presents the iterative risk management
(IRM) model to overcome the limitation of linearization; and
then Section IV.C discusses the application potentials of the
proposed approach. Section V is the test cases. Section V.A
verifies the accuracy and the efficiency of the risk gradient
computation on the RTS-96 test system. And then the RM and
IRM are tested on RTS-96 test system, a US-Canada 410-bus
system, and a 1354-bus Mid-European system, respectively.
II. G ENERIC M ODEL OF R ISK M ANAGEMENT BASED
MT AND R ISK G RADIENT

ON

A. Retro of risk assessment using Markovian tree search
The risk management needs quantification of risk. This
paper selects our previous work on efficient risk assessment
method with Markovian tree search [20], and the quasidynamic simulation method [19]. All the possible cascading
outage paths are organized in a MT, as shown in Fig. 1.
The time elapse is incorporated in the MT model, where
each level on the MT corresponds to a time interval τD , and
each node is labeled with the outage sequence from the root,
as (ik1 · · · ikn ), where ikr is either a positive integer denoting

Fig. 1. MT representation of cascading outage paths

the index of the element failed on the rth level or 0 if no
outage occurs on this level, n is the ending level of the
cascading outage. The cost corresponding to state (ik1 · · · ikn )
is C(ik1 · · · ikn ), and the conditional probability of outage
event ikr+1 after the state (ik1 · · · ikr ) is Pr(ikr+1 |ik1 · · · ikr ).
Then with the MT structure, the risk is expressed as the
following expansion [20]:
R = C0 +

X

Pr(ik1 )C(ik1 ) +

X
k1

Pr(ik1 )

Pr(ik1 )

k1

k1

+

X

X
k2

Pr(ik2 |ik1 )

X

X

Pr(ik2 |ik1 )C(ik1 ik2 )

k2

Pr(ik3 |ik1 ik2 )C(ik1 ik2 ik3 ) + · · ·

k3

(1)
Suppose control actions are taken to reduce the cascading
outage risk after the initial outages. Then the expansion (1) can
be divided into two parts: C0 is the cost of the control, and all
the other terms are the risks of subsequent cascading outages,
whose sum is denoted by R′ . The risk is assessed by searching
the MT and adding the risk terms into (1) corresponding to the
newly-visited states. To further accelerate the risk assessment,
a forward-backward tree search algorithm based on a risk
estimation index (REI) was proposed in [20].
B. Generic optimization model of risk management
Cascading outages in the early stage usually develop slowly,
so there is some time to adjust system states to reduce the
risk of potential cascading outages after the initial outage.
Control measures after the initial outages result in cost C0
while reducing the subsequent risk R′ , so a compromise
between the effect and cost of risk management should be
concerned. It is desirable that the risk of subsequent cascading
outages is confined below a certain level RS′ , while the cost
of control is minimized. Therefore, the basic formulation of
risk management can be written as
min f = C0 (x∗ )
s.t. R′ (x∗ ) ≤ RS′
∗

g(x ) ≤ 0

(2)

3

where x∗ is the target system state, C0 (x∗ ) is the cost of
control, and R′ (x∗ ) is the subsequent risk. The last constraint
represents the constraints in operations, e.g. load, generation
and transmission capacity constraints. In real power systems,
the commonly used control measures for overloading relief
or power balance include re-dispatch of generation and load
shedding initiated by operators or emergency control facilities. Moreover, to prevent voltage instability, reactive power
compensation devices can also be adjusted. In emergency, the
under-voltage load shedding, under-frequency load shedding
and controlled islanding can also be triggered to prevent
system-wide collapse. In this paper, cascading outages are
modeled and simulated with DC power flow model, and the the
re-dispatch operations for overloading relief are considered,
which includes the re-dispatch of generators and curtailment
of loads controlled by the operators. The re-dispatch is very
commonly used in the operations of power systems, and it can
be described with an optimal power flow model. The cost of
re-dispatch C0 (x∗ ) can be described as follows:
∗

C0 (x ) =

∗
−cT
D (Pd

− Pd ) +

∗
cT
G |Pg

− Pg |

(3)

Here x∗ = [P∗d T , P∗g T ]T is a vector of target re-dispatch
T T
is the pre-control system state. cD
state. x = [PT
d , Pg ]
and cG are per unit costs of load shedding and generation
adjustment, respectively. It should be noted that the proposed
risk management approach is not limited to the re-dispatch of
active power of generators and loads. For example, since the
simulation and risk assessment can also be implemented in
an AC power flow model, in that case the control of reactive
power compensators can also be realized.
However, quantifying the R′ as a function of x∗ in (2) is
not straightforward. Since cascading outages involve complex
dependent events, any change in the system state will affect all
the following states, and thus the risk terms on all the levels of
the MT are changed. Therefore, it is infeasible to analytically
quantify R′ as a function of x∗ , but the risk at the original
target state x∗0 can be linearized to obtain the risk gradient:
Γ=

∂R′
∂x∗

(4)
x∗ =x∗
0

then the risk management model (2) becomes
min f = C0 (x∗ )
s.t. Γ · (x∗ − x∗0 ) ≤ RS′ − R0′
g(x∗ ) ≤ 0

(5)

The next section will address the calculation of the risk
gradient by using the result of risk assessment.
III. C ALCULATION

OF RISK GRADIENT

A. Derivative chain of states on the MT
From (1) and (5), the risk gradient depends on the derivatives of conditional probability Pr(ikr |ik1 · · · ikr−1 ) and cost
C(ik1 · · · ikr ) of states on each level of the MT. Such calculation requires the analysis of the chain of states on the cascading
outage path. As shown in Fig.2, denote the post-outage state
on the rth level as x(r)′ , and the state after re-dispatch as x(r)∗ .

The costs of the short-timescale process and re-dispatch of all
(r)
(r)
ikr (r = 1, · · · , n) are briefly denoted as vectors CF and CR ,
(r)
(r)
respectively, and the total cost on level r is C(r) = CF +CR .
(r)
Briefly denote outage probabilities as a vector Pr .

Fig. 2. State dependencies on one level of cascading outage

In this paper, assume that the probabilities and the costs
are differentiable. If in reality these quantities are non-smooth
or discontinuous functions of states, then they need to be
treated in segments, or the sub-derivatives are used as an
approximation. And if there are discrete variables, then they
need to be treated as continuous ones temporarily. To cal(r)
(r)
culate the gradient of risk, the terms ∂Pr
and ∂C
are
∂x(0)
∂x(0)
necessary according to (1). In fact, according to Fig.3, if the
(r)
∂Pr(r)
partial derivatives on levels up to r (i.e. ∂x∂C
(r−1) , ∂x(r−1) and
(r)
(r)
∂x(r)
) are obtained, then ∂Pr
and ∂C
can be calculated
∂x(r−1)
∂x(0)
∂x(0)
iteratively. The derivation of the partial derivatives on each
level depends on the analysis of the cascading process, which
will be elucidated below.
1) Mid-term Random Outage: The probability of element
i outage on the MT is [20]:

P
λi 
− j λj τD
P
(6)
1
−
e
PrMT
=
i
j λj

where λi is the failure rate of branch i, which is assumed to
be a function of its branch flow Fi [3], [21], [22]. And Fi is
a function of the system state x. So the partial derivative of
branch outage probability on level r + 1 to x(r) is


∂Pr(r+1)
∂Pr(r+1) ∂λ 
·
· −yD MY+ , yD MY+
(7)
=
(r)
∂λ
∂F
∂x
where λ and F are vectors of λi and Fi , respectively. yD is
a diagonal matrix of branch admittances. M is a |V | × |E|
matrix and each of its column Mik corresponding to a branch
ik = {u, v} satisfies Mik ,u = 1, Mik ,v = −1 and all the
other entries are 0. Y+ is the Penrose-Moore pseudo-inverse
of admittance matrix Y.
2) Short-timescale process: As Fig. 3 shows, a shorttimescale process may comprise of several outages, and each
outage may directly lead to cost on the loss of load due to
load shedding [23]. The cost of short-timescale outages on
level r + 1 is the sum of costs caused by all outages.
nr+1
(r+1)

CF

=

X

k=1

(k|r+1)

CF

(8)

4

B. Iterative calculation of terms in risk gradient

Fig. 3. Illustration of short-timescale outage events and costs

where nr+1 is the number of short-timescale outages on level
(k|r+1)
r + 1, and the cost of the kth outage is CF
. x(k|r+1)
is the state after the kth outage on level r + 1 (note that
(k|r+1)
x(0|r+1) = x(r) and x(nr+1 |r+1) = x(r+1)′ ). Also ∂x∂x(k−1|r+1)

With the analysis in III.A and the chain rule of derivatives,
(r)
(r)
the terms ∂Pr
and ∂C
of each level r can be calculated.
∂x(0)
∂x(0)
(m) ′
(m)
, ∂x
Assume for any level m (1 ≤ m ≤ r), the terms ∂x
∂x(0)
∂x(0)
(m)∗
(r+1)
and ∂x
are obtained, then for level r + 1 the terms ∂x∂x(0) ,
∂x(0)
(r+1)∗
∂x(r+1) ′
and ∂x∂x(0) are obtained as follows
∂x(0)
∂x(r+1)′
∂x(r+1)′ ∂x(r)
=
∂x(0)
∂x(r) ∂x(0)

(15)

∂x(r+1)∗
∂x(r+1)∗ ∂x(r+1)′
=
∂x(0)
∂x(r+1)′ ∂x(0)

(16)

(k|r+1)

∂CF
∂x(k−1|r+1)

can be derived by sensitivity analysis of load
and
shedding [24], [25]. As the simulation of each outage on level
(k|r+1)
r + 1 continues, ∂x∂x(r) can be derived with:
∂x(k|r+1)
∂x(k|r+1) ∂x(k−1|r+1)
=
∂x(r)
∂x(k−1|r+1)
∂x(r)

(9)

(r+1)

∂CF
∂x(r)

(1)∗

(1) ′

X ∂C(k|r+1)
F
=
(r)
∂x
k=1

∂Pr(r) ∂x(r−1)
∂Pr(r)
=
∂x(0)
∂x(r−1) ∂x(0)

(11)
(r)

(r)

s.t. 1T (Pg − Pd ) = 0

(19)

C. Recursive form of risk gradient

(12)

Pmin
≤ Pg ≤ Pmax
g
g
where
and
are the target load and generation given by
dispatch center, P′d and P′g are the load and generation before
dispatch. Pd and Pg are the states to solve, i.e. the states after
dispatch at time τD . rg is the vector of ramping rates of all
the generation buses. (12) can be briefly denoted as follows:

X

(r+1)′

(r+1)∗

,x

, τD )

(13)

(r+1)

where p
is the parameters on the r + 1th level, such
as network topology, branch parameters, branch flow limit,
etc. x(r+1)∗ is the re-dispatch target state to fulfill, which is
determined by solving another LP problem (the conventional
DC-OPF in this paper) [20], [21] as:
x(r+1)∗ = LPa (p(r+1) , x(r+1)′ )
(r+1)

∂C

(r+1)

Pr(ikr+1 |ik1 · · · ikr )C ′ (ik1 · · · ikr+1 )

(20)
(20) shows a recursive relationship, so C ′ (ik1 · · · ikr ) could
be calculated and updated reversely from the terminal back to
the root of the MT. Similarly, define

R′ (ik1 · · · ikr+1 ) ,

∂C

Pr(ikr+1 |ik1 · · · ikr )C(ik1 · · · ikr+1 ) + · · ·

ikr+1

P∗g

,x

X

ikr+1

= C(ik1 · · · ikr ) +

P∗d ≤ Pd ≤ P′d

= LPe (p

(r)

∂CF ∂x(r−1)
∂CR ∂x(r−1)′
∂CR ∂x(r−1)∗
∂C(r)
= (r−1)
+ (r−1)′
+ (r−1)∗
(0)
(0)
(0)
∂x
∂x
∂x
∂x
∂x
∂x
∂x(0)

C ′ (ik1 · · · ikr ) , C(ik1 · · · ikr ) +

− τD rg ≤ Pg − P′g ≤ τD rg

x

(18)

Define equivalent cascading outage cost C ′ (ik1 · · · ikr ) as

∗
∗
T
min f = cT
D (Pd − Pd ) + cG Pg − Pg

(r+1)

(1)

Since ∂x
, ∂x
and ∂x
can be obtained from III.A, so
∂x(0)
∂x(0)
∂x(0)
according to (15)-(17), for any r (1 ≤ r ≤ n, where n is the
(r)∗
(r) ′
(r)
final level of cascading outage), ∂x
, ∂x
and ∂x
are
∂x(0)
∂x(0)
∂x(0)
obtained iteratively in the process of cascading outage path
(r)
(r)
simulation. So ∂Pr
and ∂C
are calculated with
∂x(0)
∂x(0)

nr+1

3) Re-dispatch: Re-dispatch is usually modeled as an optimization problem. Under the DC power flow assumption,
the execution of re-dispatch can be modeled as a linear
programming (LP) problem [20], which aims to minimize the
distance between the actual post-dispatch state and the target
state in a given time interval τD .

(r+1)

(17)

(k|r+1)

(k|r+1)

∂CF
∂CF
∂x(k−1|r+1)
= (k−1|r+1)
(10)
(r)
∂x
∂x
∂x(r)
From Fig. 3, the partial derivative of states in short-timescale
(r+1) ′
outages ∂x∂x(r) is obtained by applying (9) iteratively. And
from (8) and (10), the partial derivative of the short-timescale
outage cost is derived as

P∗d

∂x(r+1) ∂x(r+1)∗
∂x(r+1) ∂x(r+1)′
∂x(r+1)
= (r+1)∗
+ (r+1)′
(0)
(0)
∂x
∂x
∂x
∂x
∂x(0)

(r+1)

(14)
(r+1)

∂x
∂x
R
R
From (13), ∂x(r+1)
and ∂x
can
′ , ∂x(r+1)∗ , ∂x(r+1) ′
(r+1)∗
be obtained by means of Lagrange multiplier and sensitivity
(r+1)∗
analysis [26]. Similarly, ∂x
can be calculated from (14).
∂x(r+1) ′

X

Pr(ikr+1 |ik1 · · · ikr )C ′ (ik1 · · · ikr+1 ) (21)

ikr+1

With given ik1 · · · ikr , abbreviate all R′ (ik1 · · · ikr+1 ) as
vector R(r) ′ , and C ′ (ik1 · · · ikr ) as C(r) ′ , then
(r)
T
T ∂Pr
∂R′ (ik1 · · · ikr )
= C(r)′
+ Pr(r)
(0)
(0)
∂x
∂x

∂C(r)
∂R(r) ′
+
(0)
∂x
∂x(0)

!

(22)

(22) is also a recursive form. Note that R(0)′ = R′ , so
the gradient of risk can be computed with forward-backward
scheme in the risk assessment based on MT search.

5

D. Forward-backward scheme of risk gradient calculation
With (15)-(19), the partial derivatives are calculated in the
process of forward searching with the risk assessment [20], as
demonstrated in Algorithm 1.
Algorithm 1. Forward calculation of partial derivatives
Step 1. Initialize level on MT r = 0.
Step 2. Sample mid-timescale outages and calculate

∂Pr(r+1)
∂x(r)

Step 3. Simulate short-timescale outages [20], and calculate
(r+1)

∂CF

∂x(r)

with (7).
∂x(r+1) ′
,
∂x(r)

with (8)-(11).
(r+1)

(r+1)

(r+1)

Step 4. Simulate re-dispatch. Calculate ∂Pr (0) , ∂x (0) and ∂C (0) .
∂x
∂x
∂x
Step 5. If the cascading outage path ends, exit and start simulation of a
new path. Otherwise assign r = r + 1 and jump to Step 2.

After searching a cascading outage path on the MT, the risk
gradient is updated reversely, as Algorithm 2 shows.
Algorithm 2. Backward update of risk gradient.
Assume the cascade path is ik1 · · · ikn , the variables with superscript (r)
correspond to (ik1 · · · ikr ), and (0) correspond to the root state on MT.
Step 1. Define b(r) , r = 1 · · · n. If state (ik1 · · · ikr ) has been reached
before, then b(r) = 0, otherwise b(r) = 1.
(r)

Step 2. Assign r = n. Assign S(r) = b(r) ∂C(0) , ∆C (r) = b(r) C (r) .
∂x
Step 3. Reverse to the previous state on the cascade path. Assign r = r−1.
(r+1)

Step 4. Let S(r) = Pr(r+1) S(r+1) + ∆C (r+1) ∂Pr (0) .
∂x
Step 5. Let ∆C (r) = b(r) C (r) + b(r+1) Pr(r+1) ∆C (r+1) .
Step 6. If r = 0 then exit. Otherwise jump to Step 3.

established based on the generic form of (5) as
∗
∗
T
min f = −cT
D (Pd − Pd ) + cG Pg − Pg
#
" ∗
Pd − P∗d0
(0)∗
≤ RE − R′ (x0 )
s.t. − Γ ·
P∗g − P∗g0

1T P∗g − P∗d = 0

− Fmax ≤ yD MY+ P∗g − P∗d ≤ Fmax

(24)

Pmin
≤ P∗g ≤ Pmax
g
g
0 ≤ P∗d ≤ Pd

where the first constraint is the risk constraint using the
(0)∗
T
T
risk gradient. x0 = [P∗d0 , P∗g0 ]T is the target state from
(0)∗
the original dispatch strategy, R′ (x0 ) is the risk of subsequent cascades of original dispatch strategy, and RE is the
expected risk after the RM. The other constraints are the
limits of branches, generators and loads. Here the variables
are continuous, so the RM is an LP problem. If there are
discrete variables, then RM will become a mixed-integer linear
programming (MILP) problem.
The (24) reduces cascading outages risk by setting its
solution as the target state of re-dispatch. The extent of
reduction of cascading outage risk is adjusted by changing
the expected subsequent cascade risk RE . The constraint will
force the solution to a less risky state and the control cost is
expected to be higher. Denote the expected risk decrease as
(0)∗
∆R = R′ (x0 ) − RE ≥ 0, then the bigger ∆R, the more
reduction of risk is expected.
B. Iterative risk management (IRM)
Risk assessment on original
control strategy

Obtain risk gradient

Solve RM model

Risk assessment on renewed
dispatch strategy

Fig. 4. Time series of the derivation and execution of control strategy

By applying Algorithms 1 and 2 repeatedly along with the
forward searching and backward updating procedure of risk
(0) ′
assessment [20], the S(0) will converge to ∂R
. Note that in
∂x(0)
(0)
operations, x is indirectly changed by altering the dispatch
target state x∗ , as Fig. 4 shows. The gradient of risk in the
space of control variables is
Γ = S(0)

(0)

∂x
∂x∗

(23)

IV. I MPLEMENTATION OF R ISK M ANAGEMENT
A. Full optimization model of risk management (RM)
After obtaining the risk gradient in the space of control
variables, the risk management (RM) optimization model is

Y

Risk decreased?
N
End

Fig. 5. The procedure of IRM

The RM model is based on the linearization of risk at the
original operating point. When ∆R goes outside an effective region of linearization, then there will be considerable
linearization errors. To overcome such limitation, consider
iterating the procedure of RM so as to accumulate the effect of
linearization-based RM step by step. The procedure of iterative
RM (IRM) is shown as Fig. 5.
The IRM first assesses risk on the original control strategy,
and solves the RM problem (24). The new dispatch strategy is
then evaluated with risk assessment. If the risk is decreased,

6

C. Framework of RM/IRM application
Physical System

Interface/GUI

Generating Risk Management Strategy

Periodic Analysis/
Online Invocation

1.0

0.20

0.8

0.16

0.6

0.12

δdir
m

δm

then the strategy is expected to be effective. Such a procedure
is iterated until the risk does not decrease.

0.4

0.08

0.2

0.04

0.0

0.00
0

2000 4000 6000 8000 10000
Search Attempts

0

2000 4000 6000 8000 10000
Search Attempts

(b)

(a)

Risk assessment

Fig. 7. Convergence of (a) risk gradient (b) the direction of risk gradient
Risk gradient calculation

250
Solve RM strategy

Interface for
Decision Support
& Control

Real-time Strategy
Lookup

200

IRM

Risk (MW)

Power
System

Strategy Set

Strategy Database

150
100
50
0
0

Fig. 6. The framework of RM/IRM application

The RM and IRM can be used off-line to generate control
strategies on given set of system working conditions. The
generated strategies are stored in a database and can be
extracted when corresponding events occur. Moreover, the
RM and IRM also have potential of online assessment and
decision support. The framework of RM/IRM application is
demonstrated in Fig. 6.
V. C ASE S TUDIES

reduction, so in situations that require fast computation, the
number of tree search attempts can be significantly reduced.
Nevertheless, in such a case, the accuracy of estimating the
extent of risk reduction will be lower.

kΓm − Γ∗ k
kΓ1 − Γ∗ k

(25)

Γ∗
Γm
−
kΓm k kΓ∗ k

(26)

δm reflects the convergence of the vector of risk gradient,
dir
and δm
evaluates the convergence of the direction of risk
gradient. The test is conducted on the RTS-96 3-area system,
in which the parameters are set as Tmax = 150min, τD =
15min. After 10000 times of tree search, the risk, δm and
dir
δm
are all considered as converged. Fig. 7 demonstrates the
convergence profile of risk gradient.
The convergence of risk gradient is slower than that of risk
[20] (Fig. 8). This is caused by the more complicated form
of risk gradient, so the derivation of risk gradient may cost
more computation time than risk assessment. However, the
convergence of the direction of risk gradient only requires
several hundreds of search attempts, which is much faster
than the convergence of the vector of risk gradient. Actually,
only obtaining the direction of risk gradient is enough for risk

1) RTS-96 System: After validating the accuracy and the
convergence of the calculation of risk gradient, the effectiveness of RM in the RTS-96 3-area system model is tested.
Here all the cost and risk are converted into economic metrics.
Assume that adjusting 1MW of generation in an interval τD
costs $100, and 1MW load loss in an interval corresponds to
the loss of $10000. Set the initial failure on branches 22, 23
and 24. When utilizing conventional re-dispatch, the total risk
(i.e. the cost of control plus the risk of subsequent cascading
outages) is $696775. Then use the RM to reduce cascading
outage risk and evaluate the performance with risk assessment.
The effect of RM under different values of ∆R is shown in
Fig. 9. The RM effectively decreases risk within a certain
range, but the risk stops decreasing when ∆R reaches around
$600000, which means the linearization is no longer effective.
8x10 5

Total risk by RM
Expected risk

7x10 5
6x10 5

Risk ($)

According to III.C, the calculation of risk gradient can
be integrated into the forward-backward scheme of risk assessment. To efficiently derive risk management strategy, the
convergence of risk gradient is of great significance.
In the computation process, the risk gradient after the mth
search attempt is denoted as Γm , and the converged risk
gradient is denoted as Γ∗ . To evaluate the convergence profile
of risk gradient, propose the following convergence indices:

dir
δm
=

Fig. 8. Convergence profile of risk [20]

B. Effectiveness of RM

A. Convergence of risk gradient in RTS-96 system

δm =

2000 4000 6000 8000 10000
Search Attempts

5x10 5
4x10 5
3x10 5
2x10 5
1x10 5
0
0

200000

400000

600000

R($)

Fig. 9. Convergence profile of risk under RM

Next, the performance of RM with several existing methods
are compared. The risk-based OPF method (RB-OPF) [8] and

7

Risk max =5 0.9 0.8 0.7

RM
RM-1
RB-OPF

0.65
6x10 5

106

0.6

105
104
103
102
101
100
20

5x10 5

40

60

80

100

120

Line Number
0.5

4x10 5

0.1

0.05

0

5000

6000

0.2

3x10 5

Fig. 11. Risk of subsequent cascading outages by level-1 line outage.

2000

3000

4000

Control Cost C0($)

Fig. 10. Risk reduction performance of RM, RM-1 and RB-OPF

The mechanisms of RB-OPF and RM are further analyzed
and compared. In RB-OPF, as Riskmax reduces, the RB-OPF
optimization model adjusts the loading of lines to reduce the
severity of post-outage states of all the lines in the system.
However, as RB-OPF does not have accurate knowledge of
the most risky cascading outage patterns, it is not as effective
as RM for cascading outages. The risk in RB-OPF is described
as the outage probability times post-contingency severity, but
the severity does not very well reflect the actual loss (load or
economic loss) of the potential subsequent cascading outages,
so the RB-OPF tends to lower the loading of lines systemwide, disregarding that the risk under many subsequent out-

Line loading rate

2x10 5

10

10

08

08

Line loading rate

Subsequent Risk R’($)

7x10

5

ages are actually very low. As Fig. 11 shows, the distribution
of subsequent cascading outage risk in terms of level-1 line
outage is very uneven: most risk is concentrated on the level-1
outages of lines 21, 25, 49, 87 and 121 (note that there are
only 120 lines in RTS-96 system, and the "line 121" stands
for no outage in the studied time interval). Particularly, the
risk after line 25 outage accounts for about 95% of the total
risk. We select the lines whose loading rates are significantly
changed by these methods, and the relationship between line
loading rates and control cost. As shown in Fig. 12(a), the
RB-OPF reduces the risk by suppressing the loading level on
many lines even if the actual risk of the following cascading
outages is relatively trivial. In contrast, Fig. 12(b) shows
the line loading rates under RM. Since the risk assessment
identifies that the outage of line 25 is the most risky, the
RM mainly reduces the loading of line 25, and thus more
effectively reduces the overall subsequent cascading outage
risk with lower control cost. Fig. 13 compares the distribution
of risk before control, after RB-OPF and after RM. The results
show that RM identifies the "bottleneck" of cascading outage
risk and better concentrates the control effort on the risky
cascading outage patterns. As a result, the overall risk under
RM is lower at the same level of control cost. This is attributed
to the feature of the RM optimization model that makes use
of the information of risk distribution obtained from the risk
assessment and risk gradient.
Subsequent risk by line ($)

a variant of ref. [16] are tested in this case. In the RBOPF, the risk can be reduced by lowering the upper limit of
severity-based risk index Riskmax . Decreasing the Riskmax
is expected to reduce the risk. Ref. [16] changed the risk
constraint of RB-OPF/SCOPF as the expected load shedding
after N-1 outages. Since [16] is similar to the risk management
considering only the first level of cascading outages, the risk
assessment and calculating the risk gradient is conducted only
to level-1 outages, as a variant of [16]. Such a variant of [16]
is named as RM-1. Fig. 10 compares the costs of control and
subsequent cascading outage risks under the RM, RM-1 and
RB-OPF. The RM achieves a lower risk than the other two
approaches at the same level of cost, so it is verified that the
RM is more effective in reducing the cascading outage risk.
The RM-1, regarded as an simplified version of RM, does
not consider the multi-level cascading outages, so the accuracy
of reducing the cascading outage risk is not as satisfactory as
RM. It should be noted that although the result of RM-1 in Fig.
10 still seems fine, the effectiveness of RM-1 actually highly
depends on the distribution of the risk on different levels of
cascading outages. In this test case, the risk mainly exists on
the level-1 outages, so the RM-1 does not have significant
difference from RM. But if substantial risk is on other levels
of outages, the accuracy of RM-1 will be significantly affected.
Fig. 10 further annotates the cost-risk relationship with different values of Riskmax in the RB-OPF. The result shows that
the relationship between cascading outage risk and Riskmax
is non-smooth, which will cause difficulties in managing the
cascading outage risk. By comparison, the RM proposed in
this paper achieves a better performance in reducing risk, and
the amount of risk reduction is more controllable.

06
04
02

06
04
02
00

00
3000

4000

5000

Control cost C0 ($)
(a) RB-OPF

6000

3000

4000

5000

Control cost C0 ($)
(b) RM

Line No.
30
31
25
27
26
15
20
19
116
32
33
11
29
59

Fig. 12. Loading rates of selected lines under (a) RB-OPF, (b) RM.

2) US-Canada Northeast system: Next, the RM is tested
on a real-world system model, i.e. the US-Canada Northeast
power grid model [19] retaining all generator buses and the
buses at 230 kV or above. The system contains 410 buses
(287 load buses and 233 generator buses, note that a bus may
be both load bus and generator bus) and 882 branches. In
this case, increase the system load to 1.5 times of the base

8

Subsequent risk ($)

than RB-OPF, and the average computation time is less than
1s. However, since the optimization model of RM requires
risk gradient, the computation time of risk assessment and risk
gradient calculation should also be included. In this case, each
risk assessment uses 200 Markovian tree search attempts, and
the average computation time is longer than that of RB-OPF.

Before control
RB-OPF
RM

7x10 5
6x10 5
5x10 5
4x10 5
3x10 5
2x10 5
1x10 5

TABLE I
T IME CONSUMPTION COMPARISON IN US-C ANADA N ORTHEAST SYSTEM

0
21

25

49

87

121

Line number

Time consumption (s)
RM
RB-OPF

Subprograms
Fig. 13. Cascading outage risks of selected level-1 outages.

Subsequent cascade risk ($)

load and set initial outages on branches 320 and 321. Fig.
14 compares the cost-risk profiles obtained by RM and RBOPF, respectively. It can be seen that the RM reduces the
risk by over 70%, while the RB-OPF only reduces the risk
by about 10%. The effect of RB-OPF is limited due to its
mechanism. According to the definition of severity function in
RB-OPF, if the loading rates of a line under all N-1 states are
below 90%, then the risk of the line is 0. So when setting the
Riskmax = 0, the RB-OPF will guarantee the loading rates
of all the lines under all N-1 states under 90%. However, in
this case, substantial risk may still exist under some cascading
outage patterns, which cannot be reduced by RB-OPF. In
contrast, since RM uses the information obtained from the risk
assessment of cascading outages, it more effectively reduces
the risk of subsequent cascading outages.
1.0x10 5
0.0
1

Ri
sk
ma
x

=0

3 1.5 0.8
2 1

95000

8.0x10 4

90000

6.0x10 4

85000

0.5
0.4
0

4.0x10 4

2000

0.2
0.3

0.05
0.1

4000

0.0
0

10000

20000

30000

40000

79.41
25.40
0.36

–
–
72.29

Total

105.17

72.29

C. Risk-cost coordination realized by IRM
1) RTS-96 System: The performance of the RM is limited
only within a range where linearization is effective. Fig. 9
shows that with the RM, the risk stops decreasing at around
$2.5×105. IRM keeps updating risk gradient at new operating
points and the risk is further reduced. Table II and Fig. 15
demonstrates the cost-risk relationship derived by IRM.
TABLE II
C OST- RISK IN THE ITERATION PROCESS IN RTS-96 TEST SYSTEM
Round

∆R ($)

Control
cost($)

Subsequent
risk($)

Total
risk($)

0
1
2
3
4
5

–
600000
100000
100000
20000
1000

2475
3789.8
4141.6
5034.9
5456.4
5590

694300
249590
174840
40738
12434
13132

696775
253379.8
178981.6
45772.9
17890.4
18722

6000

RM
RB-OPF

2.0x10 4

Risk assessment
Computation of risk gradient
Solving RM/RB-OPF model

50000

The results indicate that the IRM can effectively overcome
the limitation of linearization with the RM and further reduce
risk of cascading outages. After 4 rounds of iterations, the
subsequent risk of cascading outages has reduced by 97.3%.

Control cost ($)

The computational efficiency of RB-OPF and RM is also
tested. Table I shows the average time of RM and RB-OPF
in US-Canada northeast system. The methods are developed
with MATLAB and are tested on a PC with 2.6GHz CPU
and 32GB RAM. Similar with the computational strategies
proposed by [10], the RB-OPF is realized with the technique of
benders decomposition to handle large-system case. Since the
RB-OPF considers the severity function of all the lines under
all the N-1 outages, the scale of the optimization problem
grows in proportion to the square of the number of branches
in the system. In the US-Canada northeast system with 882
branches, the average computation time of RB-OPF is 72.29s.
While the optimization model of RM has much smaller scale

Subsequent Risk R’($)

Fig. 14. Risk reduction by RM and RB-OPF in US-Canada northeast system

7x10 5

0

6x10 5
5x10 5
4x10 5
3x10 5

1
2

2x10 5
1x10

5

0
2000

3
3000

4000

4* 5

5000

6000

Control Cost C0($)

Fig. 15. Cost-risk trajectory of IRM in RTS-96 test system

2) US-Canada Northeast system: Next, the IRM is tested
on the US-Canada Northeast power grid model. The cost-risk
characteristics derived by IRM are demonstrated in Fig. 16.

9

Risk ($)

8x10

0

Subsequent cascade risk
Total risk

4

1

6x10 4

2 3

4 5

4x10 4

6*

7

2x10 4
0
0

10000

20000

30000

40000

Control cost C0($)

Fig. 16. Cost-risk trajectory of IRM in US-Canada Northeast system

The results indicate that the total risk at the 6th iteration
reaches the lowest among all iterations, where the total risk
is expected to decrease by 54.5%, and the subsequent cascade
risk drops significantly by 93.6%. The subsequent cascade risk
is even lower in the 7th iteration, but the drop in subsequent
cascade risk is offset by a substantially high cost, which causes
a higher total risk than that of the 6th iteration. In practice,
the adopted strategy for risk management may vary depending
on the risk preference. Therefore, the strategies in Fig. 16 can
also be regarded as the results of multi-objective optimization
of risk and cost, and the control strategy for risk management
can be selected depending on the risk preference.
Regarding the computational efficiency, the average performance of a single IRM run in the US-Canada Northeast system
model is shown in Table III. The result shows that solving the
RM models takes only several seconds, but much more time is
consumed in cascading outage simulation, risk assessment and
calculating the risk gradient. The speed of computation can be
further enhanced with parallel computing on high-performance
computation platform, and this method also has potentials for
on-line application on a period of 5-15 minutes for operators’
decision support to prevent cascading outages.
TABLE III
T IME CONSUMPTION OF IRM IN US-C ANADA N ORTHEAST SYSTEM
Subprograms of IRM

Time consumption (s)

Cascading outage simulation & risk assessment
Computation of risk gradient
Solving RM model
Total

660.41
237.42
3.48
901.31

As for the drawbacks of the proposed approach, it is
observed that the accuracy of risk gradient calculation may
decrease as the size of the power system grows. This is because
as the system size grows, the number of nonlinear behaviors
(e.g. switched active constraints in dispatch model, etc.) also
grows. To maintain desired accuracy, more iterations in IRM
may be necessary, and the computation speed will be adversely
affected. The estimation of effective linearization region and
the derivation of desirable step size in IRM will be studied
in the future. Moreover, since the calculation of risk gradient
(r) ′
(r)∗
, ∂x
and
requires to store all the sensitivity matrices ∂x
∂x(0)
∂x(0)
(r)
∂x
, the memory usage in large-scale systems will be high.
∂x(0)
This problem can be alleviated by using compressed storage

(CS) of the sensitivity matrices, since most elements in the
sensitivity matrices have very low absolute values (in large
systems, generally less than 1% of the elements have absolute
values larger than 10−3 , and less than 10% have absolute
value larger than 10−5 ). Compressing the matrices by dropping
low-absolute-value elements practically does not affect the
accuracy, but can significantly decrease memory usage.
3) 1354-bus Mid-European System [27]: To further demonstrate the performance, the proposed approach is tested on
a larger Mid-European backbone system model that appears
in [27], [28]. The system has 1354 buses (including 260
generation buses and 688 load buses) and 1991 branches. The
branch flow limits were modified to secure N-1. We maintain
the same number of MT search attempts as the US-Canada
Northeast System case, and run the IRM. The effectiveness
of IRM is shown in Fig. 17, and the computational time is
shown in Table IV. The computation can be finished within
1 hour with our desktop computer, which shows promising
potential for online decision support in practical systems. Also,
the speed can be further enhanced with parallel computing
on high-performance platforms (currently the approach is
developed and tested on Matlab and has not been further
optimized for performance).
Subsequent cascade risk
Total risk

0
4x10 4

1
2

Risk ($)

1x10 5

3x10 4

3
4

2x10 4

5*

6

1x10 4
0
5000

10000

15000

20000

Control cost C0($)

Fig. 17. Cost-risk trajectory of IRM in Mid-European system

TABLE IV
T IME CONSUMPTION OF IRM IN M ID -E UROPEAN SYSTEM
Subprograms of IRM

Time consumption (s)

Cascading outage simulation & risk assessment
Computation of risk gradient
Solving RM model
Total

1806.17
1449.69
6.079
3261.94

For large scale systems, the CS can be adopted to avoid
memory overflow. Table V compares the performance with
and without CS. The threshold of absolute sensitivity value
is set as 10−5 . Although the CS causes some computational
overhead in compressing and indexing matrices, the memory
usage is significantly reduced. Per the observation from the
test results on more systems, the computation time complexity
when using the CS is approximately O(N 2 ), where N is the
number of buses. Moreover, it is observed that the number of
remaining sensitivity elements after CS grows approximately
linearly with system size, while the full sensitivity matrix

10

grows quadratically with system size. Such a desirable spatial
complexity of CS contributes to the less overall computation
time in large-scale systems cases, as Table V shows. Therefore,
the CS is recommended for large systems.
TABLE V
T IME & M EMORY USAGE OF RISK GRADIENT COMPUTATION
Systems

Without CS
Time(s)
Memory(GB)

With CS
Time(s)
Memory(GB)

US-Canada
Mid-European

237.42
2184.60

670.10
1449.69

4.87
30.25

0.76
3.32

VI. C ONCLUSION
This paper proposes an approach for risk management of
cascading outages based on risk gradient and MT search. The
expansion of risk corresponding to the MT structure is linearized in the space of control variables and then the gradient
of risk is obtained. The computation of risk gradient adopts an
efficient forward-backward algorithm, which can be combined
with the procedure of risk assessment based on Markovian tree
search, and the risk gradient has good convergence profile.
With the risk gradient, the constraint of risk is established
and incorporated into a dispatch model to formulate the risk
management (RM) optimization model. The RM minimizes
control cost while limiting the cascading outage risk under a
given level for a desirable trade-off between cost and risk. The
risk and cost in the RM model have clear physical meanings
and thus can practically give insights for operators’ decision
support. The effectiveness of RM is verified on the RTS-96
3-area system and a 410-bus US-Canada system model.
To overcome the limitation of linearization in the RM, the
iterative RM (IRM) approach is proposed to achieve more
effective reduction of cascading outages risk. The test cases on
the US-Canada northeast system and 1354-bus Mid-European
system verify the effectiveness of the RM and IRM. Moreover,
compressed storage (CS) of sensitivity matrix technique is
introduced to significantly reduce memory usage without compromising the accuracy, which further improves the practicality
of the proposed approach in large-scale system applications.
The RM and IRM may be utilized in on-line decision support
for preventive control against potential cascading outages
following possible initial contingency scenarios.
R EFERENCES
[1] P. Hines, J. Apt, and S. Talukdar, “Large blackouts in north america:
Historical trends and policy implications,” Energy Policy, vol. 37, no. 12,
pp. 5249–5259, 2009.
[2] B. A. Carreras, V. E. Lynch, D. E. Newman, and I. Dobson, “Blackout
mitigation assessment in power transmission systems,” in Proceedings of
the 36th Annual Hawaii International Conference on System Sciences.
IEEE, 2003, pp. 10–pp.
[3] J. Chen, J. S. Thorp, and I. Dobson, “Cascading dynamics and mitigation
assessment in power system disturbances via a hidden failure model,”
International Journal of Electrical Power & Energy Systems, vol. 27,
no. 4, pp. 318–326, 2005.
[4] D. E. Newman, B. A. Carreras, V. E. Lynch, and I. Dobson, “Exploring
complex systems aspects of blackout risk and mitigation,” IEEE Trans.
Reliab., vol. 60, no. 1, pp. 134–143, 2011.

[5] D. Bertsimas, E. Litvinov, X. A. Sun, J. Zhao, and T. Zheng, “Adaptive robust optimization for the security constrained unit commitment
problem,” IEEE Trans. Power Syst., vol. 28, no. 1, pp. 52–63, 2013.
[6] M. Vaiman, K. Bell, Y. Chen, B. Chowdhury, I. Dobson, P. Hines,
M. Papic, S. Miller, and P. Zhang, “Risk assessment of cascading
outages: Methodologies and challenges,” IEEE Trans. Power Syst.,
vol. 27, no. 2, pp. 631–641, 2012.
[7] D. Kirschen and D. Jayaweera, “Comparison of risk-based and deterministic security assessments,” IET Generation, Transmission & Distribution, vol. 1, no. 4, pp. 527–533, 2007.
[8] F. Xiao, Risk based multi-objective security control and congestion
management. ProQuest, 2007.
[9] Q. Wang, J. D. McCalley, T. Zheng, and E. Litvinov, “A computational
strategy to solve preventive risk-based security-constrained opf,” IEEE
Trans. Power Syst., vol. 28, no. 2, pp. 1666–1675, 2013.
[10] Q. Wang and J. D. McCalley, “Risk and n-1 criteria coordination for
real-time operations,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 3505–
3506, 2013.
[11] Q. Wang, G. Zhang, J. D. McCalley, T. Zheng, and E. Litvinov, “Riskbased locational marginal pricing and congestion management,” IEEE
Trans. Power Syst., vol. 29, no. 5, pp. 2518–2528, 2014.
[12] E. Karangelos and L. Wehenkel, “Probabilistic reliability management
approach and criteria for power system real-time operation,” in Power
Systems Computation Conference (PSCC), 2016. IEEE, 2016, pp. 1–9.
[13] C. Hamon, M. Perninge, and L. Söder, “A computational framework for
risk-based power systems operations under uncertainty. part i: Theory,”
Electric Power Systems Research, vol. 119, pp. 45–53, 2015.
[14] D. N. Simopoulos, S. D. Kavatza, and C. D. Vournas, “Reliability
constrained unit commitment using simulated annealing,” IEEE Trans.
Power Syst., vol. 21, no. 4, pp. 1699–1706, 2006.
[15] M. A. Ortega-Vazquez, “Assessment of nk contingencies in a probabilistic security-constrained optimal power flow,” in Power and Energy
Society General Meeting (PESGM), 2016. IEEE, 2016, pp. 1–5.
[16] F. Capitanescu, “Enhanced risk-based scopf formulation balancing operation cost and expected voluntary load shedding,” Electric Power
Systems Research, vol. 128, pp. 151–155, 2015.
[17] J. He, L. Cheng, D. S. Kirschen, and Y. Sun, “Optimising the balance
between security and economy on a probabilistic basis,” IET Generation,
Transmission & Distribution, vol. 4, no. 12, pp. 1275–1287, 2010.
[18] R. Yao, X. Zhang, S. Huang, S. Mei, Z. Zhang, X. Li, and Q. Zhu, “Cascading outage preventive control for large-scale ac-dc interconnected
power grid,” in IEEE PES General Meeting. IEEE, 2014, pp. 1–5.
[19] R. Yao, S. Huang, K. Sun, F. Liu, X. Zhang, and S. Mei, “A multitimescale quasi-dynamic model for simulation of cascading outages,”
IEEE Trans. Power Syst., vol. 31, no. 4, pp. 3189–3201, July 2016.
[20] R. Yao, S. Huang, K. Sun, F. Liu, X. Zhang, S. Mei, W. Wei, and
L. Ding, “Risk assessment of multi-timescale cascading outages based
on markovian tree search,” IEEE Trans. Power Syst., vol. 32, no. 4, pp.
2887–2900, July 2017.
[21] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, “An initial
model for complex dynamics in electric power system blackouts,” in
Proceedings of the Annual Hawaii International Conference on System
Sciences. Citeseer, 2001, pp. 51–51.
[22] M. Yang, J. Wang, H. Diao, J. Qi, and X. Han, “Interval estimation
for conditional failure rates of transmission lines with limited samples,”
IEEE Trans. Smart Grid, 2016.
[23] H. Ren, I. Dobson, and B. A. Carreras, “Long-term effect of the n-1
criterion on cascading line outages in an evolving power transmission
grid,” IEEE Trans. Power Syst., vol. 23, no. 3, pp. 1217–1225, 2008.
[24] C. Reddy, S. Chakrabarti, and S. Srivastava, “A sensitivity-based method
for under-frequency load-shedding,” IEEE Trans. Power Syst., vol. 29,
no. 2, pp. 984–985, 2014.
[25] A. A. Girgis and S. Mathure, “Application of active power sensitivity
to frequency and voltage variations on load shedding,” Electric Power
Systems Research, vol. 80, no. 3, pp. 306–310, 2010.
[26] P. R. Gribik, D. Shirmohammadi, S. Hao, and C. L. Thomas, “Optimal
power flow sensitivity analysis,” IEEE Trans. Power Syst., vol. 5, no. 3,
pp. 969–976, 1990.
[27] C. Josz, S. Fliscounakis, J. Maeght, and P. Panciatici, “Ac power flow
data in matpower and qcqp format: itesla, rte snapshots, and pegase,”
arXiv preprint arXiv:1603.01533, 2016.
[28] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power
systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1,
pp. 12–19, 2011.

