This article has been accepted for inclusion in a future issue of this journal.

Content is final as presented, with the exception of pagination.

IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS 1

Modeling the Dynamics of Cascading

Failures in Power Systems
Xi Zhang, Choujun Zhan, and Chi K. Tse, Fellow, IEEE

Abstract— In this paper, we use a circuit-based power flow and control systems [1]. Protection equipment is responsible
model to study the cascading failure propagation process, and for maintaining reliability through applying switching actions
combine it with a stochastic model to describe the uncertain of relays and circuit breakers. Relays are essential auxiliary
failure time instants, producing a model that gives a complete
dynamic profile of the cascading failure propagation beginning components in transmission lines, generators, transformers,
from a dysfunctioned component and developing eventually to and other kinds of power apparatus. When a relay detects an
a large-scale blackout. The sequence of failures is determined abnormal operating condition such as over-current and voltage
by voltage and current stresses of individual elements, which dip, it will switch off the affected component to remove the
are governed by deterministic circuit equations, while the time fault from the network, thereby preventing further damages
durations between failures are described by stochastic processes.
The use of stochastic processes here addresses the uncertainties and hence ensuring the normal operation of the rest of the
in individual components’ physical failure mechanisms, which system. The on/off states of relays determine the structure of
may depend on manufacturing quality and environmental factors. the power network, thus influencing the overall operational
The element failure rate is related to the extent of overloading. state of a power system. Normally the power grid is designed
A network-based stochastic model is developed to study the to maintain its power distribution function even when a few
failure propagation dynamics of the entire power network. Sim-
ulation results show that our model generates dynamic profiles elements are removed [2]. However, when the power grid is
of cascading failures that contain all salient features displayed under stressed conditions, for instance, due to heavy loads
in historical blackout data. The proposed model thus offers and outages of equipment, the removal of some elements may
predictive information about occurrences of large-scale blackouts. lead to huge disturbances and subsequent tripping of other
We further plot cumulative distribution of the blackout size to elements, causing a possible severe blackout [1].
assess the overall system’s robustness. We show that heavier loads
increase the likelihood of large blackouts and that small-world The dynamic cascading failure process in a power grid can
network structure would make cascading failure propagate more be viewed as a sequence of tripping events, leading eventually
widely and rapidly than a regular network structure. to power outage affecting a very large area. It has been
Index Terms— Complex network, power system, dynamics of observed that the 1996 Western North America blackouts [3],
cascading failure, power flow study, stochastic process. the 2003 Northeastern America and Canadian blackouts [4]
and other historical blackout data all display a typical profile
I. I NTRODUCTION characterized by a relatively slow initial phase followed by a
sharp escalation of cascading failures. Such a universal form

E LECTRICITY supply network is an essential part of the

infrastructure of modern society. Large power blackouts
cause inconvenience to residents in the affected areas as well
of dynamic profiles strongly suggests that a common model
can be used to describe the dynamic cascading failure process.
The study of the dynamic propagation of cascading fail-
as considerable economic loss to the community at-large. ures provides useful hints for system vulnerability detec-
Power system’s security has always been an issue of serious tion, robustness assessment and network control. Recently,
concern among utility providers, infrastructure developers, Chen et al. [5] used a generalized Poisson model, neg-
and policy makers, and is also a topic attracting significant ative binomial model and exponentially accelerated model
attention of electrical engineers and researchers. The power to generate the probabilities of the propagation of trans-
distribution network is a complex and highly interconnected mission outages which fit the observed historical data.
network, consisting of power apparatus, protection equipment Dobson et al. [6], [7] used branching processes to analyze the
Manuscript received October 12, 2016; revised December 22, 2016; propagation of cascading failures in power grids. Much of the
accepted January 19, 2017. This work was supported by the Hong Kong previous work primarily applied data fitting methods to investi-
Research Grants Council GRF under Project PolyU 5258/13E. This paper gate the statistical characteristics of power systems’ blackouts,
was recommended by Guest Editor H. Ho-Ching Iu.
X. Zhang and C. K. Tse are with the Department of Electronic and but fell short of considering the essential electrical circuit
Information Engineering, Hong Kong Polytechnic University, Hong Kong operations or the impact of the network structure. Moreover,
(e-mail: xi.r.zhang@connect.polyu.hk; encktse@polyu.edu.hk). cascading failures in power grids have also been studied in
C. Zhan was with The Hong Kong Polytechnic University, Hong Kong.
And now he is with the Department of Electronics Communication and terms of the sequential trippings of electrical elements in real
Software Engineering, Nanfang College of Sun Yat-sen University, Guangzhou networks. Among the many switching mechanisms of relays,
510900, China (e-mail: zchoujun2@gmail.com). overloading is the most prominent one and has been widely
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org. studied [8]–[13]. In a cascading failure process, the failure of
Digital Object Identifier 10.1109/JETCAS.2017.2671354 one element leads to power flow redistribution in the grid,
2156-3357 © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

2 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

which can cause some other elements to be overloaded. These

overloaded electrical elements can then be tripped by their
relays, causing another round of failures until the remaining
elements are all within their respective operating limits.
Various of tripping sequence settings of the overloaded
elements have been explored [9], [10], [14]. Typically, in each
round of the cascading simulations, the power flow distrib-
ution in the network is computed, and overloaded electrical
elements are removed at the same time. The actual time delays Fig. 1. Dynamic description of failure in terms of state transitions. State “0”
is the normal connected state; state “1” is the removed or tripped state.
and dynamical profiles of the process are not considered in
this kind of models, making them unable to simulate the
dynamic propagation of a cascading failure. In order to show
the network is dependent on the sum of individual extents
the dynamic profile, some previous studies made the simple
of overloading of all elements in the network. Simulation
deterministic assumption that the duration for an overloaded
results show that the cumulative number of failed elements
 t +t to be tripped is equal to t which is given by
element
triggered by some initial failures shows a universal growing
t ( f j (τ ) − f¯j )dτ = o j , where f j is the power flow
pattern which is consistent with historical blackout data. Thus,
of overloaded element j , f¯j is the flow limit and o j is a
our model can offer insights into the mechanism of cascading
specific threshold of that element [15], [16].
propagation in a power system as well as provide predictive
Considering the complexities and uncertainties in real power
information for the failure spreading in the network. Our study
grids [17], a few researchers turned to use probabilistic models
also includes the effects of loading conditions and network
to characterize the tripping events of the elements in power
structure on the extent and rapidity of blackouts in power
grids [11], [12], [18]. For instance, Wang et al. [11] used
systems. The UIUC 150 bus system with different consumer
a Markov model to study cascading failures, where the trip-
load distributions and several types of network structure are
pings of the elements are regarded as state transitions, which
studied for comparison purposes. It is shown that heavy load
are memoryless and probabilistic. In Wang et al.’s work,
conditions increase the risk of large blackouts in the same
the overloaded elements share one same tripping rate, which
power system, and that small-world network structure is more
is much larger than the natural failure rate of the electrical
prone to rapid propagation of cascading failures than the
equipment. Alternatively, an overall state transition probability
regular structure.
can be determined by considering the maximum capacity
of the failed elements and a random tripping process [12]. II. FAILURE M ECHANISMS OF C OMPONENTS
It is also shown [19] that a component will experience more
failures under heavy load conditions. The varying tripping A power system is composed by various electrical stations
rates for elements under different extents of overloading stress connected by transmission lines, and each station or trans-
have not been thoroughly considered in the aforementioned mission line is protected by protective equipment. In this
stochastic models. Study of essential collective behavior of paper, we model electrical stations as nodes and transmission
a power network must be pursued according to the gov- lines as links, with nodes being connected by links forming a
erning physical laws which in the case of power systems power network [10]. Deterministic power flow equations are
should involve circuit-based power flow equations [10] (see used to generate the sequence of failures and their locations.
Appendix). By suitably combining the power flow study A node or link is a basic element of a power network.
with probabilistic methods for describing inevitable uncertain- We refer to an element’s tripping event as an element state
ties, the dynamic profile of cascading failure processes can transition (EST). The cascading failure propagation in a power
be realistically revealed, hence offering important predictive network can be viewed as a sequence of ESTs in the network.
information about the occurrence of large-scale blackouts. In this section, we investigate the state transition behavior of
In this paper, we study the dynamics of cascading failure a basic element, and in the next section, we apply probabilis-
propagations in power systems. The key contributions are as tic theory to study the collective transition behavior of the
follows. First, we apply circuit-based power flow equations to network.
determine the sequence of failures in accordance to the extent
of overloadings of individual components. In order to describe A. Time to Failure of a Basic Element
the complete dynamic profile, we need to determine the time Let si (t) be the state of element i of a given network,
durations between failures in the propagation sequence. Due and si (t) ∈ [0, 1], with si (t) = 0 corresponding to a
to the complexities and uncertainties of the involving physical connected element i at time t, and si (t) = 1 corresponding
failure mechanisms of the components (e.g., manufacturing to a removed (tripped or open-circuited) element i at time t,
quality, environmental factors, etc.), stochastic processes are as shown in Fig. 1. Here, λi (t) is the rate of transition of node i
used to model the dynamic changes. Then, to study the going from state “0” to “1”, and μi (t) is the transition rate
collective behavior of the entire system in terms of failure from “1” to “0”. Then, the future state of an element is solely
propagation in the whole network, an extended chemical determined by its present state and the transition rule. Suppose
master equation (CME) model is used [20], [21]. Based on the present time is t, and dt is an infinitesimal time interval.
the CME model, we show that the failure propagation rate of As si (t) ∈ {0, 1}, P{si (t + dt) = 1} and P{si (t + dt) = 0}
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 3

can be written separately as where L i (t) is the power loading of element i that can be
found from the power flow calculation, Ci is the capacity
P[si (t + dt) = 1]
of that element, and ai is the basic unit rate (trippings
= P[si (t + dt) = 1 |si (t) = 0 ]P[si (t) = 0]
per second). For normal operating condition, λ1i = 0. In a
+ P[si (t + dt) = 1 |si (t) = 1 ]P[si (t) = 1] cascading failure process, λ1i  λ0i [24]. Without loss of
P[si (t + dt) = 0] generality, we assume that λi (t) ≈ λ1i (t) in our analysis of
= P[si (t + dt) = 0 |si (t) = 0 ]P[si (t) = 0] cascading failures in power systems.
+ P[si (t + dt) = 0 |si (t) = 1 ]P[si (t) = 1] (1) For the sake of completeness, we also allow a
tripped or removed element to be repaired, and hence be
where P[si (t) = 1] and P[si (t) = 0] denote the probability
restored to its normal connected state. Thus, we define μi (t)
that node i is in state “1” and “0” at time t, respectively;
as the transition rate of element i going from state “1” to “0”
P[si (t + dt) = 1 |si (t) = 0) ] is the conditional probability
as a result of repair actions or self-healing ability of the power
that given si (t) = 0 element i transits to state “1” in the time
system. In practice, an element’s state cannot be switched
interval (t, t +dt); and P[si (t +dt) = 0 |si (t) = 1) ] is defined
arbitrarily. Also, the time delay for recovering a tripped
in a likewise manner. Using the state transition rates shown in
element should be considered and can be included in the actual
Fig. 1, P[si (t + dt) = 1 |si (t) = 0) ] can be written as
representation of μ (t). This recovery process can be used to
P[si (t + dt) = 1 |si (t) = 0 ] = λi (t)dt. (2) study the power restoration process after the power blackout.
In this paper, we focus on analyzing the cascading failure
Also, P[si (t +dt) = 0 |si (t) = 0) ] is the probability that given
process. Thus, considering that not all elements could be
si (t) = 0, element i remains in state “0” in time interval
repaired in a short time and an element cannot keep changing
(t, t + dt) (i.e., no state transition occurs). Thus, we have
its status frequently, we take μi (t) as 0 for a fast cascading
P[si (t + dt) = 0 |si (t) = 0 ] = 1 − λi (t)dt. (3) process.
Likewise, we have
C. Power Flow Calculation
P[si (t + dt) = 0 |si (t) = 1 ] = μi (t)dt, (4)
In addition to equation (7), power flow calculation is
P[si (t + dt) = 1 |si (t) = 0 ] = 1 − μi (t)dt. (5)
still needed for the analysis of cascading failures. Several
algorithms and tools are available for computing power
B. State Transition Rates of Basic Elements flows [25], [26]. The actual power system is a high-order
In this section, we discuss the physical meanings of element complex nonlinear network, and any abrupt change of network
state transition rates λi (t) and μi (t) in a fast cascading failure structure can change the power flow distribution, and at the
process. In statistical terms, an event rate refers to the number same time cause large transients, oscillations, and bifurca-
of events per unit time. Specifically, λi (t) is the rate of tions [27]. Using our definition of state transition of elements,
element i becoming disconnected in the network which is the tripping probability of each element is an integration of the
caused by either a natural equipment malfunction or tripping tripping rate (extent of overloading) with time. In this study,
by its protective equipment, i.e., we assume that the system can always reach a steady state
λi (t) = λ0i (t) + λ1i (i ) (6) when tripping occurs and that the transient before the system
reaches the next steady state is sufficiently short, making
where λ0i (t) is the equipment malfunctioning rate in the accumulative effects negligible. As far as the propagation
absence of loading stress and its value is constant and derivable of cascading failures is concerned, it suffices to consider
from past statistics [16]; and λ1i (t) is the removal or tripping blackouts caused by overloading, ignoring the nonlinear char-
rate by protective relays and is determined by the (over)- acteristics of the circuit elements and possible oscillatory
loading condition and the capacity of element i . behavior. In our previous work, a model that can accurately
Among the many tripping mechanisms of relays [22], [23], track the load change in a power network during a cascading
power overloading is a dominant one. In this study, we focus failure has been developed [10]. This model is adopted in our
on switching actions caused by overloading. When the load study here. Given information of power consumption, power
of element i is within its capacity, it is assumed to work in generation and grid topology, the voltage of each node can be
the normal condition and will not be removed or tripped by found using
the protective relay, namely λ1i (t) = 0. However, when the
element exceeds its capacity, there will be a short delay before A∗V = B (8)
it is finally removed. The tripping rate is relevant to the extent
of overloading. In other words, if there is a large overloading where A is a matrix describing the power network; B =
T
of element i , it will be tripped more rapidly compared to the · · · Ii 0 v k 0 · · · with Ii and v k representing the sink
case of a light overloading [19]. Based on this assumption, current and voltage of nodes i and k, respectively. A brief
we can write λ1i (t) as description is given in the Appendix and more details can be
⎧   found in [10]. The current flowing through the transmission
⎨a L i (t) − Ci , if L (t) > C line (i, j ) can be calculated as
i i i
λ1i (t) = Ci (7)
⎩
0, if L i (t) ≤ Ci Ii j = (v i − v j ) ∗ Yi j (9)
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

4 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

where v i and v j are voltages of nodes i and j , and Yi j is

the admittance of line (i, j ). Equation (8) is derived from
consideration of circuit laws and thus rigorously describes the
behavior of the power network.
Fig. 2. Time line of network state transitions.
III. FAILURE P ROPAGATION IN THE N ETWORK
A power network is represented as an undirected graph G
3) The probability that two or more element state transitions
consisting of m elements. The state of G is defined as
occur after dt is given by
S = {s1 , s2 , ..., sm }, which is a vector containing the states
of all m elements. Network G can have 2m possible network P[S(t + dt) = R S |S(t) = N S ] = 0 (12)
states, and any state transition of an element will lead to a
network state transition of G. where R S denotes the network state that two or more of the
The dynamic propagation of cascading failures in G is “0”-state elements in N S become “1”. From equation (12),
equivalent to the dynamic evolution of S(t). Given the current there is at most one element state transition at a time.
state of the network, the network state transition can be
described by (i) the time of the next state transition; and B. Extended Gillespie Method
(ii) identification of the next element that will transit (be In this section, we derive S(t) using an extended Gillespie
tripped). method [28], which was used for analyzing coupled chemical
reactions [20], [21].
A. Basics As shown in Fig. 2, the state of the power system at t1
First, we consider the network state transitions in an infin- is N S , i.e., S(t1 ) = N S . Let Q(τ ) denote the probability that
itesimal time interval dt. Suppose S(t) = N S , which is a given S(t1 ) = N S , no transition occurs in (t1 , t1 + τ ), i.e.,
specific network state among the 2m possible states. Thus,
Q(τ ) = P[S(t1 + τ ) = N S |S(t1 ) = N S ] (13)
S(t + dt) is the network state after a duration of dt. Only
those elements in state “0” may transit, leading to a network Similarly, Q(τ + dt) can be written as
state transition. Let 0 be the set of elements in state “0”, and
1 be the set of removed (tripped) elements. From elementary Q(τ + dt)
probability theory, we have the following basic results: = P[S(t1 + τ + dt) = N S |S(t1 ) = N S ]
1) Omitting O(dt), the probability that no element under- = P[S(t1 + τ + dt) = N S |S(t1 + τ ) = N S ]Q(τ ) (14)
goes a state transition after dt can be written as
Given S(t1 ) = N S , power flow calculation can be performed,
P[S(t + dt) = N S |S(t) = N S ] as described in Section II-C, and λi (t1 ) can be derived based
= [1 − λi (t)dt)] on the settings in Section II-B. If no state transition occurs
i∈0 during time interval (t1 , t1 + τ ), we have S(t) = S(t1 ) and
= 1− λi (t)dt + λx1 (t)λx2 (t)(dt)2 λi (t) = λi (t1 ) for t ∈ (t1 , t1 + τ ). From (10), we get
i∈0 x 1 ,x 2 ∈0
P[S(t1 + τ + dt) = N S |S(t1 + τ ) = N S ] = 1− λi (t1 )dt
− λx1 (t)λx2 (t)λx3 (t)(dt)3 + · · · i∈0
x 1 ,x 2 ,x 3 ∈0 (15)
= 1− λi (t)dt + O(dt) ≈ 1 − λi (t)dt (10)
i∈0 i∈0 Thus, by putting (15) in (14), we get

where x 1 , x 2 , · · · are the elements in 0 . Q(τ + dt) = Q(τ )(1 − λ∗ (t1 )dt) (16)
2) The probability that only one element state transition (say
where λ∗ (t1 ) = i∈0 λi (t1 ). Furthermore, re-arranging (16)
element k) occurs after dt, i.e., only element k transits, can
and taking the limit dt → 0, we get
be written as
d Q(τ ) Q(τ + dt) − Q(τ )
P[S(t + dt) = M S |S(t) = N S ] = lim = −λ∗ (t1 )Q(τ ),
dτ dt →0 dt
= λk (t)dt [1 − λi (t)dt] ⇒ Q (τ ) = −λ∗ (t1 )Q(τ ). (17)
i∈0 /{k}
= λk (t)dt − λk (t)λx1 (t)(dt)2 The probability that nothing happens in zero time is one,
i.e., Q(0) = P{S(t1 ) = N S |S(t1 ) = N S } = 1. Then,
x 1 ∈0 /{k}
the analytical solution of (17) is
+ λk (t)λx1 (t)λx2 (t)(dt)3 + · · · ∗ (t
x 1 ,x 2 ∈0 /{k} Q(τ ) = e−λ 1 )τ . (18)
= λk (t)dt + O(dt) ≈ λk (t)dt (11)
Let h i (τ, dt) denote the probability of the event that given
where x 1 , x 2 , · · · are the elements in 0 /{k} and M S denotes S(t1 ) = N S , the next transition occurs in the interval (t1 +
the network state that only one of the “0”-state elements in N S τ, t1 + τ + dt) in element i . There are two conditions for
becomes “1”. this event to occur. The first condition is that there is no state
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 5

transition during (t1 , t1 + τ ). The second condition is that a

state transition occurs in element i during (t1 + τ, t1 + τ + dt).
Thus, h i (τ, dt) can be written as
h i (τ, dt) = P[S(t1 + τ + dt) = M S |S(t1 + τ ) = N S ]Q(τ )
Fig. 3. Relative probability for elements in 0 to be first tripped given
(19) S(t1 ) = N S .

Putting (11) and (18) in (19), we get

∗ (t In our stochastic model, any overloaded element in 0 may
h i (τ, dt) = e−λ 1 )τ λi (t1 )dt (20)
be tripped first. From a probabilistic viewpoint, the element
Let H (τ, dt) denote the probability that the next transition with a higher h i (τ, dt) will more likely be tripped first. Thus,
occurs in the time interval (t1 + τ, t1 + τ + dt), given we define the relative probability for element i (i ∈ 0 ) to be
S(t1 ) = N S . It is readily shown that tripped first as:
H (τ, dt) = h i (τ, dt) = λ∗ (t1 )e−λ
∗ (t
1 )τ
dτ (21) h i (τ, dt) λi (t1 )
r fi = = ∗ . (26)
i∈0 H (τ, dt) λ (t1 )
Further, let τ denote the time interval between two adjacent where λ∗ (t1 ) = i∈0 λi (t1 ). Our model can incorporate this
network state transitions, and f (τ ) denote the state transition tripping order using the following steps:
probability density function (PDF): 1) A random number z 2 is generated uniformly in (0,1).
H (τ, dt) − H (τ, 0) ∗
2) Suppose there are l overloaded elements in 0 . With no
f (τ ) = li m = λ∗ (t1 )e−λ (t1 )τ (22) loss of generality and for ease of referral, let these over-
dt →0 dt
loaded elements be elements L1, L2, . . . , L j, . . . , Ll.
i.e., Fig. 3 shows the relative probability of an overloaded
f (τ ) = λ∗ (t1 )e−λ
∗ (t
1 )τ (23) element in 0 to be first tripped, given that S(t1 ) = N S .
3) The j th element in 0 is selected to be tripped according
The accumulative probability density function that the next to
transition occurs before time t1 + τ , given S(t1 ) = N S , can be j −1 j
written as λ Lk λ Lk
 z2 < (27)
τ λ∗ λ∗
∗ −λ∗ (t1 )t −λ∗ (t1 )τ k=0 k=0
F(τ ) = λ (t1 ) e dt = 1 − e (24)
0 where λ L0 = 0.
Note that one can also get F(τ ) from F(τ ) = 1 − Q(τ ).
Equations (23) and (24) show that τ follows an exponential IV. C ASCADING FAILURE S IMULATIONS
AND PARAMETERS
distribution and that the network transition rate is λ∗ (t1 ). Here,
λ∗ (t1 ) is the sum of the element state transition rates of all In this section we describe the simulation algorithm and
the working elements in the network, and is determined by some important characterizing parameters of our model that
the sum of the extents of overloading of all the overloaded are relevant to predicting the occurrence of power blackouts.
elements. The time interval τ is expected to be short when
λ∗ (t1 ) is large, i.e., the network state transition (cascading A. Simulation Algorithm
process) occurs very rapidly. Thus, the physical meaning of Fig. 4 shows the flow chart for simulating the cascading
λ∗ (t1 ) can be interpreted as the overloading stress of the entire failure process which can be summarized as follows:
power system. 1) Initial Settings: At the start of the simulation, all volt-
In order to include this characteristic in our model, we take ages at the power generation stations, currents flowing
the following steps to determine the time of the next network into the consumer nodes, admittances of the transmission
state transition, given S(t1 ) = N S : lines, and capacities of elements are set.
1) A random number z 1 is generated uniformly in (0,1). 2) Initial Failure: An initial failure is planted by remov-
2) Let F(τ ) = z 1 , and τ is derived as ing one element from the network, which triggers the
cascading failure process.
ln(1 − z 1 )
τ= . (25) 3) Iterative Process: Based on S(t), we remove the tripped
−λ∗ (t1 ) elements from the network, and keep all elements
whose states is “0”. The remaining network may be
C. Order of State Transition disconnected, forming so-called islands, due to the
A number of working elements (elements in 0 ) can removal of the tripped elements. For a disconnected
possibly undergo state transition. In our analysis pre- sub-network (island) containing no generator node, all
sented in Section III-A, we allow only one element to be elements within it would have no access to power and all
removed (tripped) at a time. Pfitzner et al. [29] pointed out power flows become zero. All nodes in this sub-network
that the order in which overloaded lines are tripped influences are unserved. Note that these elements are not tripped,
the cascade propagation significantly. In this section, we study and their states are still “0”. Moreover, for a sub-network
the order in which element state transitions take place. containing at least one generator node, equation (8) can
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

6 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

Fig. 5. Typical propagation profile and onset time tonset . Maximum

propagation rate gm occurs at t = tm .

elements at time t (NoTE(t)) is used. Here, we take NoTE(t)

as the number of “1”-state elements in S(t). Also, the number
of elements not served is another important metric used to
measure the blackout size. As the power grid’s operation
depends on the connection of the elements, the tripping of
some elements in the network can disconnect the grid and
island the consumer nodes from the power sources. These
nodes are being deprived of power, and are labelled as
unserved nodes in our analysis. We use NoUN(t) to denote the
number of cumulative unserved nodes at t. To find NoUN(t),
Fig. 4. Flow chart for simulating cascading failure. we remove the tripped elements in G, and identify all sub-
networks in the remaining part of G. The consumer nodes
that are isolated from generators are all unserved nodes.
be used to compute the power flow distribution in this
Furthermore, during a cascading failure process, it is par-
sub-network. Power flows of all the “0”-state elements
ticularly important to track the growing rate of the num-
in G can be computed, and the tripping rate of each
ber of failed or tripped elements, i.e., the frequency of
element λi can be obtained using (7). If all tripping
removal or tripping of overloaded elements. Specifically, any
rates are positive, we determine the next network state.
rapid increase in the frequency of removal of overloaded
Specifically, we first determine the time of the next net-
elements is a precursor to an onset of a large blackout. Thus,
work state transition using (25), and determine the ele-
a metric that effectively gives the critical time from which
ment in 0 that will be tripped next. The network state
tripping begins to take place more rapidly is extremely relevant
transition is determined using (27). Then, we update
to prevention of power blackouts. This metric, called onset
S(t) = S(t + τ ), and iterate the process until all the
time (tonset ) here, can simply be defined as the time after which
transition rates are found to be zero (i.e., no overloaded
the propagation rate of the cascading failure increases rapidly,
elements). With no more overloaded elements in the
as depicted in Fig. 5. In other words, tonset is a critical time
network, no state transition will occur and S(t) is a
point before which remedial control and protection actions
stable state. We can then end the simulation and get
should be applied to the power grid. After tonset , the power grid
the final network.
undergoes a short phase of very rapid tripping of overloaded
elements leading to large power blackout within a very short
B. Parameter Settings and Metrics time. To compute the onset time, we identify the maximum
The time of the initial failure is set as zero, and the time of rate of the growing profile by solving d 2 NoTE(t)/dt 2 = 0,
the final network state transition (after which there are no over- which gives t = tm as the time point where the growing rate is
loaded elements in the network, and the network state enters highest. Assuming that the value of dNoTE(t)/dt at t = tm is
a stable point) is tfinal . Using the above algorithm, we can gm and the initial phase has a very slow growing rate, the onset
simulate the dynamic profile of S(t) for power network G, time is simply given by
from t = 0 to t = tfinal . For t > tfinal , S(t) remains unchanged. NoTE(tm )
The dynamic profile of S(t) is thus the dynamic propagation of tonset = tm − (28)
gm
cascading failures in the network. In order to better represent
and visualize the characteristics of the dynamics of a cascading In practice, we can use any handy algorithm to find tonset ,
failure, we use the following metrics, which are extracted for instance, by locating the time instant where dNoTE(t)/dt
from S(t). starts to increase rapidly. Section V-B offers one simple
We propose several metrics to investigate the cascading algorithm.
failure in a power system. Finally, to characterize the general severity in the event of
First, to characterize the propagation profile of a cascading a possible blackout, we use the following statistical metric.
failure in a power system, the cumulative number of tripped Suppose a large number of cascading failure cases, initialized
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 7

by failures of different elements in the network, are simulated.

The probability that the blackout size of a randomly picked
case is larger than a chosen threshold BS (Blackout Size) is
given by
n[x(t) ≥ BS]
P[x(t) ≥ BS] = (29)
n
where x(t) can be NoTE(t) or NoUN(t), n is number of
the total blackout cases simulated, and n[x(t) ≥ BS] is the
number of cases whose blackout size at t is larger than BS.
Based on (29), we can evaluate the cumulative blackout size
distribution of a network to reveal the probability (risk) of
having a blackout of a specific level of severity in a given
network.

V. A PPLICATION C ASE S TUDY

In this section, we simulate cascading failures in the
UIUC 150 Bus System using the model proposed above.
The UIUC 150 Bus is a power test case offered by Illinois
Center for a Smarter Electric Grid at UIUC [30]. It contains
150 buses and 217 links that operate in 3 different voltage base
values. We merge the parallel lines that connect the same two
buses into one link, resulting in 203 links in our simulation.
We assume that the current sinks of the consumer buses given
in the UIUC 150 Bus are the normal load demands of these
consumers. The voltages of generators are all 1.04 p.u. based
on the data in the test case.
From the historical blackout reports [3], [4], one can find
that the tripped elements are mostly generators, transmission
lines and transformers. Thus, in our simulation, we set current
limits for the transmission lines, transformers and the genera-
tors according to Ci = (1 + α) ∗ Ii (normal), where Ii (normal)
is the current flowing through a transformer or a transmission Fig. 6. Propagation profile of the Western North America power blackout
line, or the total current flowing out of a generator under in (a) July 1996; (b) August 1996.
normal load demand condition; α is the safety margin and
is set to 0.2. The current limits of other elements (consumer the same fashion, the cascading failure propagated slowly at
buses and distribution buses) are set to values that are large the beginning, but 6000 seconds later, the propagation rate
enough to avoid tripping during a cascading failure. accelerated sharply. The main propagation finished in 120 sec-
onds. The 2003 American-Canada blackout [4] also showed
A. Dynamics of Cascading Failure Propagation a similar growing pattern, with NoTE growing very slowly
We first study the failure spreading during a blackout for the initial 4 hours and then accelerating rapidly to its
process. Fig. 6(a) shows the profile of cumulative tripped final state. The rapid increase in NoTE occurred in a few
elements of the blackout in the Western North American minutes, which was a small fraction of the whole cascading
system in July 1996 [3]. The blackout started from the failure period (0, tfinal ).
of the 345 kV Jim Bridger-Kinport line (the time of that In the following, we use the proposed stochastic model
initial failure is 0 in the figure). As shown in the Fig. 6(a), to simulate the dynamic propagation of cascading failures
NoTE grew very slowly, at the initial phase, until the failing triggered by the failure of one single line. First, we simulate
of the 230 kV Brownlee-Boise Bench line at 1600 seconds 100 different propagation profiles of the cascading failure
after the initial failure. Then, the cascading failure speeded up process triggered by the initial failure of line (2, 21), under
abruptly, and within 380 seconds, NoTE reached 33 from 6 at a normal load demand condition and the condition with 5%
the end of the initial phase. Finally, the cascading failure increase in load demands. Note that when the loading of
settled at a final state where 34 major elements were tripped, the power system is increased by 5% and S(t = 0) = 0,
depriving 10% consumers of the Western interconnection area there are no overloaded elements, i.e., the 5% increase in
from access to electrical power. Fig. 6(b) shows the profile load demands will not cause any outage in the power grid.
of NoTE in another blackout of the same power system that From the data of historical blackouts, the time duration
occurred in August 1996. The cascading failure was triggered of the failure propagation is usually between 1 hour and
by the failure of the 500 kV Big Eddy-Ostrander line, and then 4 hours. In this simulation, we use a uniform ai for all the
continued with a sequence of tripping of elements. In almost elements in the UIUC 150 Bus system, and fit ai to make the
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

8 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

TABLE II
S EQUENCE OF E LEMENT T RIPPING E VENTS

growing pattern as the historical blackout data. Using

Fig. 7. Simulation of the dynamics of a cascading failure event in the equation (23), the growth rate of the cascading failure is
UIUC 150 Power System caused by an initial failure of line (2, 21). determined by λ∗ (t). Fig. 7(b) shows the values of λ∗ (t)
(a) NoTE and NoUN; (b) λ∗ .
throughout the simulated cascading process. Initially, the value
of λ∗ (t) is relatively small, until the breakdown of some
TABLE I critical elements, its value increases very rapidly. This means
S IMULATION R ESULTS FOR THE C ASCADING FAILURE that the power network operates under a high overloading
T RIGGERED BY THE FAILURE OF L INE (2, 21) stress. When the stress comes down again, the propagation
slows down. The tripping of elements ceases when λ∗ (t)
reduces to 0, and the network reaches its final condition.
The consistency of our simulated cascading failure process
with the historical data verifies the validity of our model in
describing realistic blackout processes. The following two key
issues should be noted.
averaged tfinal of the 100 simulated results under the condition (1) We use stochastic methods to investigate the cascading
of 5% increase in load demands to be 3600 s. Thus, ai is failure propagation. Our model takes into consideration the
set as 0.035 s−1 in our simulation. Table I lists the averaged high complexities as well as uncertainties of the involving
simulated values of NoTE(tfinal), NoUN(tfinal ) and tfinal . From mechanisms which can be investigated with probabilistic
Table I, we see that under a normal load demand condition, methods. The tripping rates of the elements are related to
the failure of line (2, 21) will not cause severe disturbance to the overloading extents of the corresponding elements and the
the power system. When the network is stressed by heavier more heavily overloaded ones will be more likely to be tripped
loads, the failure of the same transmission line can lead to a first. Equations (25) and (27) incorporate these considerations.
large blackout in the power system. Thus, for the same system and same initial failure, different
Table II lists the sequence of the element tripping events in simulations may yield different results due to the stochastic
one simulated cascading failure process under the condition nature of the model. Fig. 7 is one particular simulation run
of 5% increase in load demand. We plot the profiles of of the UIUC 150 Bus with initial failure of line (2, 21).
NoTE and NoUN in Fig. 7(a), which display the same typical Furthermore, Figs. 8 (a)-(f) show results derived from 6 other
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 9

Fig. 8. Simulation of failure propagations in UIUC 150 Bus power system with initial tripping of line (2, 21) using the proposed model. (a)-(f) Six separate
simulation runs.

Fig. 9. Simulation of failure propagation initiated by (a) failure of line (2, 14) in UIUC 150 Bus; (b) failure of line (34, 35) in UIUC 150 Bus; (c) failure
of line (103, 105) in IEEE 118 Bus.

simulation runs. From Fig. 8, we can see that these 6 sets of failure arises from a generator which is the only generating
results share the same characteristic profile, where the growing unit in its power network. Moreover, the failure propagation
rate of the blackout size is uneven and a relatively slow profile shown in Fig. 7 is unique for power systems, which
initial phase is followed by a sharp escalation of cascading is determined by the specific failure spreading mechanism.
failures. Such profile is not normally observed in other failure spreading
(2) From our simulations, we observe cascading failure pat- mechanisms, such as disease propagations in human networks,
terns as shown in Fig. 7. Figs. 9 (a) and (b) show the simulated rumour spreading on the Internet, and so on.
cascading failure propagation in the UIUC 150 Bus initiated
by the failure of line (2, 14) and the failure of line (34, 35), B. Blackout Onset Time
respectively. Fig. 9 (c) shows the simulated cascading failure To evaluate tonset , we adopt an intuitive algorithm that
propagation in the IEEE 118 Bus initiated by the failure of locates the time point at which NoTE begins to escalate
line (103, 105). It should be noted that not all initial failures rapidly. Suppose this time point is tk which corresponds to the
will generate such cascading failure profiles. In fact, the initial time when the kth element is tripped. The gradient of NoTE
failure of some elements will not cause further cascading before this time point is gi = (k − 1)/tk , and the gradient
failures at all. Another extreme case is that initial failure of of NoTE after this time point is gm = w/(tk+w − tk ), where
some crucially important element in the network will make all w is an arbitrary additional number of elements tripped after
other elements to be unserved instantly. For instance, the initial the kth time point for the purpose of computing the gradient.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

10 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

of NoTE(tfinal ) and NoUN(tfinal ) are plotted using equa-

tion (29). These simulations are repeated for the conditions
that the load demands are increased by various percentages.
Figs. 11(a) and (b) show the cumulative distributions for
NoTE(tfinal ) and NoUN(tfinal ) of the UIUC 150 Bus System
under several different load demand conditions. Specifically,
all simulated cascading failures will result in NoTE(t) > 0,
as any cascading failure simulation has a tripped element as
initial failure. Thus, for BS threshold = 0, P[x(t) ≥ 0] = 1
where x(t) is NoTE(t). We see that under a normal load
condition, the probability of large blackouts of the UIUC
150 Bus System is relatively low. However, when the load
demands are increased, the probability of large blackouts
increase significantly. We also observe that the probability
of having severe blackouts does not grow in linearly with
Fig. 10. Probability density function of tonset . the growth of load demands. From Figs. 11(a) and (b),
P[x(tfinal ) ≥ BS] grows relatively slowly relatively when the
TABLE III load demand increases by less than 5%, but more rapidly when
C ONFIDENCE I NTERVALS OF tonset the load demand increases by 5% to 10%.

D. Effects of Network Structure

It has been shown that the network topology plays a
significant role in determining the dynamics of propagation
and spreading of disease or information in networks [31]. For
example, infectious disease spreads more readily in small-
world networks than in regular ones [31]. It is shown that the
topological characteristics of many real-world power systems
We compare the two gradients, and if gm > γ gi , where γ > 0, are not uniform [32]. Thus, it is meaningful to investigate
we accept this kth time point as the onset time. the relationship between network structure and functional
We perform 10,000 simulations of cascading failures trig- properties of power systems, and to identify better connec-
gered by removal of line (2, 21). In our algorithm, we use tivity styles. However, since the mechanism of infectious
w = 10 and γ = 50 to find tonset of these 10,000 simulation disease or information spreading is totally different from that
runs and analyze the probability density function of tonset . of failure cascading in power systems, the conclusion derived
From Fig. 10, we see that there is a peak in the inter- in prior studies cannot be applied directly to the power systems
val (1200 s,1400 s), implying that tonset is more likely to be directly.
around 1300 s. Also, Monte Carlo method is applied to derive Regular networks and small-world networks are used as test
the confidence interval of tonset for three different confidence power systems for comparison purposes. A regular network is
levels, as listed in Table III. It should be noted that the tonset generated with 150 nodes, each node’s degree being 4. We
distribution shown in Fig. 10 is only valid for the cascading allocate 30 generators in the regular network, whose voltages
failures in the UIUC 150 Bus power system with initial failure are set as 1.04 p.u. The remaining nodes are consumer nodes,
of line (2, 21). Different systems should have different tonset each sinking 0.3 p.u. of current, and the admittances of all
distributions, which should be derived from computation on links in this network are set as 2 × 103 p.u. Then, we generate
the specific power systems. 3 small-world networks by rewiring the links in the regular
network with rewiring probabilities 0.1, 0.2 and 0.3 [31].
For each test network, we increase the consumers’ load
C. Effects of Heavy Load Demands demands by 5%, and then simulate 10 profiles of the cascading
Another common characteristic of the three historical black- failure initiated by the failure (removal) of each line. We plot
outs is that they all took place in the hot summer when the the cumulative blackout size distributions for the 4 networks
power demand is high. In this section, we investigate the based on equation (29). Fig. 12(d) shows the cumulative
overall influence of load demands on blackout risk of a power distribution of NoTN(tfinal ), and we see that the risk of large
network. The cumulative blackout size distribution is used final blackouts is higher for small-world networks than for
to indicate the risk of severe power blackouts of the power regular ones. In order to show the speed of the propagations
system. in these networks, we plot the accumulative distributions of
Under a normal load demand condition, we simulate NoTE at different time points. Figs. 12(a), (b) and (c) show
10 profiles of cascading failure triggered by the failure of the cumulative distributions of NoTE at 100 s, 500 s and
each line. Thus, we have altogether 2030 blackout cases for 1000 s, respectively. For the same duration [0, t0 ], a higher
the UIUC 150 Bus System. Then, we analyze the profile value of P{NoTE(t0 ) ≥ BS} for the same BS indicates a faster
data of these 2030 blackouts. The cumulative distributions speed of cascading. From Fig. 12, we can conclusion that the
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 11

Fig. 11. Cumulative blackout size distributions of UIUC 150 Power System. Blackout size measured in (a) NoTE; (b) NoUN.

Fig. 12. Cumulative blackout size distributions. (a) t = 100 s; (b) t = 500 s; (c) t = 1000 s; (d) t = tfinal . “pbeta” is rewiring probability for generating
small-world networks.

cascading failure propagates faster in small-world networks combining deterministic power flow equations and stochastic
than in regular networks. time duration descriptions. An extended chemical master equa-
tion method is adopted to analyze the network failure dynam-
VI. C ONCLUSIONS ics. It has been verified that the model produces propagation
In this paper, we develop a model to investigate the profiles that contain the key features displayed in historical
dynamics of the cascading failure processes in power systems, blackout data. We studied the UIUC 150 Bus system and
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

12 IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS

a few important representative network structures with the R EFERENCES

model, and identified the effects of heavy load demands and
[1] S. H. Horowitz and A. G. Phadke, Power System Relaying. New York,
network structure on the rapidity of propagation of possible NY, USA: Wiley, 2008.
blackouts. We also develop metrics to evaluate the risk of [2] M. Rausand and A. Høyland, System Reliability Theory: Models, Statis-
large-scale blackouts in terms of cumulative blackout size tical Methods, and Applications, 2nd ed. New York, NY, USA: Wiley,
2004.
distributions. The model described in this paper thus provides [3] E. C. Eakeley et al., “1996 system disturbances,” North Amer. Elect.
predictive information for possible power blackout events in Rel. Council, Disturb. Anal. Working Group, Atlanta, GA, USA,
power systems. Tech. Rep., 2002. [Online]. Available: http://www.nerc.com/pa/rrm/ea/
System%20Disturbance%20Reports%20DL/1996SystemDisturbance.pdf
[4] U.S.-Canada Power System Outage Task Force, “Final report
A PPENDIX on the August 14, 2003 blackout in the United States and
P OWER F LOW C ALCULATION Canada: Causes and recommendations,” U.S. Dept. Energy Nat.
Res. Canada, Tech. Rep., 2004. [Online]. Available: https://
Four kinds of nodes are considered in a power network: www3.epa.gov/region1/npdes/merrimackstation/pdfs/ar/AR-1165.pdf
consumer nodes, distribution nodes, transformer nodes and [5] Q. Chen, C. Jiang, W. Qiu, and J. D. McCalley, “Probability mod-
generation nodes. Voltage and currents in the network are els for estimating the probabilities of cascading outages in high-
voltage transmission network,” IEEE Trans. Power Syst., vol. 21, no. 3,
computed according to the kinds of nodes as described below. pp. 1423–1431, Aug. 2006.
[6] I. Dobson, B. A. Carreras, and D. E. Newman, “Branching process
A. Consumer Nodes (Loads) models for the exponentially increasing portions of cascading failure
blackouts,” in Proc. 39th Annu. IEEE Hawaii Int. Conf. Syst. Sci.,
A consumer node i dissipates power, and at the circuit level, Jan. 2005, p. 64.
it sinks current Ii and is unregulated. The current value is [7] I. Dobson, “Estimating the propagation and extent of cascading line
outages from utility data with a branching process,” IEEE Trans. Power
negative as the node consumes power, i.e., Syst., vol. 27, no. 4, pp. 2146–2155, Feb. 2012.
 [8] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, “An initial
−Yi1 · · · Yii · · · −Yin ∗ V = Ii (30) model for complex dynamics in electric power system blackouts,”
 T in Proc. 34th Annu. IEEE Hawaii Int. Conf. Syst. Sci., Jan. 2001,
where V = · · · v i v j v k v h · · · ; v n and In are the pp. 710–718.
voltage and external injected current at node n, respectively; [9] S. Pahwa, C. Scoglio, and A. Scala, “Abruptness of cascade failures in
Yi j is the admittance of the transmission line connecting nodes power grids,” Sci. Rep., vol. 4, Jan. 2014, Art. no. 3694.
[10] X. Zhang and C. K. Tse, “Assessment of robustness of power systems
i and j; and Yii = j =i Yi j . If there is no transmission line from a network perspective,” IEEE Trans. Emerg. Sel. Topics Circuits
between nodes i and j, Yi j = 0. Syst., vol. 5, no. 3, pp. 456–464, Sep. 2015.
[11] Z. Wang, A. Scaglione, and R. J. Thomas, “A Markov-transition model
for cascading failures in power grids,” in Proc. 45th IEEE Hawaii Int.
B. Distribution Nodes Conf. Syst. Sci., Jan. 2012, pp. 2115–2124.
A distribution node j is modelled as a transfer node without [12] M. Rahnamay-Naeini, Z. Wang, N. Ghani, A. Mammoli, and
M. M. Hayat, “Stochastic analysis of cascading-failure dynamics in
any external injected current. power grids,” IEEE Trans. Power Syst., vol. 29, no. 4, pp. 1767–1779,
 Jul. 2014.
−Y j 1 · · · Y j j · · · −Y j n ∗ V = 0 (31) [13] O. Yağan, “Robustness of power systems under a democratic-fiber-
bundle-like model,” Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat.
C. Generation Nodes Interdiscip. Top., vol. 91, no. 6, p. 062811, 2015.
[14] P. Crucitti, V. Latora, and M. Marchiori, “A topological analysis of
A generation node k is a fixed voltage source. The current the Italian electric power grid,” Phys. A, Statist. Mech. Appl., vol. 338,
emerging from this node depends on its voltage, the power nos. 1–2, pp. 92–97, 2004.
[15] M. J. Eppstein and P. D. H. Hines, “A ‘random chemistry’ algorithm for
consumption of other nodes and the network topology. The identifying collections of multiple contingencies that initiate cascading
nodal equation is failure,” IEEE Trans. Power Syst., vol. 27, no. 3, pp. 1698–1705,
 Aug. 2012.
0 · · · yk · · · 0 ∗ V = v k (32) [16] P. Rezaei, P. D. H. Hines, and M. J. Eppstein, “Estimating cascading
failure risk with random chemistry,” IEEE Trans. Power Syst., vol. 30,
where yk = 1, and v k is the voltage of node k. no. 5, pp. 2726–2735, Sep. 2015.
[17] H. M. Merrill and A. J. Wood, “Risk and uncertainty in power system
planning,” Int. J. Elect. Power Energy Syst., vol. 13, no. 2, pp. 81–90,
D. Transformer Nodes Feb. 1991.
Transformer nodes connect the high-voltage grids with [18] R. Yao et al. (Mar. 2016). “Risk assessment of multi-timescale cas-
cading outages based on Markovian tree search.” [Online]. Available:
mid-voltage or low-voltage grids. Usually two equations are https://arxiv.org/abs/1603.03935
required to describe the relations of electrical parameters [19] Y. Sun, P. Wang, L. Cheng, and H. Liu, “Operational reliability assess-
around a transformer node h. In this study, we perform our ment of power systems considering condition-dependent failure rate,”
IET Generat., Transmiss. Distrib., vol. 4, no. 1, pp. 60–72, Jan. 2010.
analysis in per unit (p.u.). Thus, one equation is sufficient: [20] D. T. Gillespie, “A general method for numerically simulating the

−Yh1 · · · Yhh · · · −Yhn ∗ V = 0 (33) stochastic time evolution of coupled chemical reactions,” J. Comput.
Phys., vol. 22, no. 4, pp. 403–434, Dec. 1976.
Combining equations (30) to (33), we get the following power [21] D. T. Gillespie, “Exact stochastic simulation of coupled chemical
reactions,” J. Phys. Chem., vol. 81, no. 25, pp. 2340–2361, 1977.
system equation: [22] K. Bae and J. S. Thorp, “A stochastic study of hidden failures in
power system protection,” Decision Support Syst., vol. 24, nos. 3–4,
A∗V = B (34) pp. 259–268, Jan. 1999.
[23] F. Yang, A. P. S. Meliopoulos, G. J. Cokkinides, and Q. B. Dam, “Effects
where A is the collection
 of all the row vectors in the equations of protection system hidden failures on bulk power system reliability,”
T
of all nodes, and B = · · · Ii 0 v k 0 · · · . in Proc. 38th IEEE North Amer. Power Symp., 2006, pp. 517–523.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

ZHANG et al.: MODELING THE DYNAMICS OF CASCADING FAILURES IN POWER SYSTEMS 13

[24] J. Chen, J. S. Thorp, and M. Parashar, “Analysis of electric power system Chi K. Tse (M’90–SM’97–F’06) received the
disturbance data,” in Proc. 34th Annu. IEEE Hawaii Int. Conf. Syst. Sci., B.Eng. (Hons.) degree in electrical engineering and
Jan. 2001, pp. 738–744. the Ph.D. degree from the University of Melbourne,
[25] A. R. Bergen and V. Vittal, Power Systems Analysis. Upper Saddle River, Australia, in 1987 and 1991, respectively.
NJ, USA: Prentice-Hall, 1999. He served as Head of the Department of Elec-
[26] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, tronic and Information Engineering, The Hong Kong
“MATPOWER: Steady-state operations, planning, and analysis tools Polytechnic University, Hong Kong, from 2005 to
for power systems research and education,” IEEE Trans. Power Syst., 2012, where he is currently Chair Professor. He has
vol. 26, no. 1, pp. 12–19, Feb. 2011. authored/co-authored ten books, 20 book chapters,
[27] F. Milano, L. Vanfretti, and J. C. Morataya, “An open source power and over 500 papers in research journals and con-
system virtual laboratory: The PSAT case and experience,” IEEE Trans. ference proceedings, and holds five U.S. patents.
Educ., vol. 51, no. 1, pp. 17–23, Feb. 2008. He received a number of research and industry awards, including the Prize
[28] C. Zhan, C. K. Tse, and M. Small, “A general stochastic model Paper Awards by the IEEE T RANSACTIONS O N P OWER E LECTRONICS
for studying time evolution of transition networks,” Phys. A, Statist. in 2001 and 2015, the Best paper Award by the International Journal of Circuit
Mech. Appl., vol. 464, pp. 198–210, Dec. 2016. Theory and Applications in 2003, the RISP Journal of Signal Processing Best
[29] R. Pfitzner, K. Turitsyn, and M. Chertkov. (Apr. 2011). “Controlled trip- Paper Award in 2014, two Gold Medals at the International Inventions Exhibi-
ping of overheated lines mitigates power outages.” [Online]. Available: tion in Geneva in 2009 and 2013, a Silver Medal at the International Invention
https://arxiv.org/abs/1104.4558 Innovation Competition in Canada in 2016, and a number of recognitions
[30] UIUC 150-Bus System, accessed on Jul. 2016. [Online]. Available: by the academic and research communities, including honorary professorship
http://icseg.iti.illinois.edu/synthetic-power-cases/uiuc-150-bus-system/ by several Chinese and Australian universities, the Chang Jiang Scholar
[31] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ Chair Professorship, the IEEE Distinguished Lectureship, the Distinguished
networks,” Nature, vol. 393, no. 6684, pp. 440–442, 1998. Research Fellowship by the University of Calgary, the Gledden Fellowship and
[32] G. A. Pagani and M. Aiello, “The power grid as a complex network: International Distinguished Professorship-at-Large by the University of West-
A survey,” Phys. A, Statist. Mech. Appl., vol. 392, no. 11, pp. 2688–2700, ern Australia. He received the President’s Award for Outstanding Research
2013. Performance twice, the Faculty Research Grant Achievement Award twice,
the Faculty Best Researcher Award, and several teaching awards with The
Xi Zhang received the B.Eng. degree in electri- Hong Kong Polytechnic University.
cal engineering from Beijing Jiaotong University, Dr. Tse serves as panel member of the Hong Kong Research Grants Council
and NSFC, and member of several professional and government committees.
Beijing, China, in 2013. He is currently pursuing
the Ph.D. degree with the Department of Electronic He is on the Editorial Boards of a few other journals. He has served and
and Information Engineering, The Hong Kong Poly- serves as the Editor-in-Chief of the IEEE T RANSACTIONS ON C IRCUITS
AND S YSTEMS II (2016–2017), the IEEE Circuits and Systems Magazine
technic University, Hong Kong.
His research interests include nonlinear analysis of (2012–2015), the Editor-in-Chief of the IEEE C IRCUITS AND S YSTEMS
power systems, modeling of electrical networks, and S OCIETY N EWSLETTER (since 2007), an Associate Editor of three IEEE
applications of complex networks in the assessment Journal/Transactions, and the Editor of the International Journal of Circuit
Theory and Applications.
of robustness of power systems.

Choujun Zhan received the B.S. degree in

automatic control engineering from Sun Yat-sen
University, Guangzhou, China, in 2007, and the
Ph.D. degree in electronic engineering from City
University of Hong Kong in 2012. He was a
Post-Doctoral Fellow with The Hong Kong Poly-
technic University. Since 2016, he has been an Asso-
ciate Professor with the Department of Electronic
Communication and Software Engineering, Nanfang
College of Sun Yat-Sen University, Guangzhou,
China. His research interests include complex net-
works, collective human behavior and systems biology.