2. THE NATURE OF CASCADING FAILURE BLACKOUTS

Bulk electrical power transmission systems are complex networks of large numbers of components that interact in diverse ways. For example, most of America and Canada east of the Rocky Mountains is supplied by a single network running at a shared supply frequency. This network includes thousands of generators, tens of thousands of transmission lines and network nodes, and about 100 control centers that monitor and control the network flows. The flow of power and some dynamical effects propagate on a continental scale. All of the electrical components have limits on their currents and voltages. If these limits are exceeded, automatic protection devices or the system operators disconnect the component from the system. We regard the disconnected component as failed because it is not available to transmit power (in practice, it will be reconnected later). Components can also fail in the sense of misoperation or damage due to aging, fire, weather, poor maintenance, or incorrect design or operating settings. In any case, the failure causes a transient and causes the power flow in the component to be redistributed to other components according to circuit laws and subsequently redistributed according to automatic and manual control actions. The transients and readjustments of the system can be local in effect or can involve components far away, so that a component disconnection or failure can effectively increase the loading of many other components throughout the network. In particular, the propagation of failures is not limited to adjacent network components. The interactions involved are diverse and include deviations in power flows, frequency, and voltage, as well as operation or misoperation of protection devices, controls, operator procedures, and monitoring and alarm systems. However, all of the interactions between component failures tend to be stronger when components are highly loaded. For example, if a more highly loaded transmission line fails, it produces a larger transient, there is a larger amount of power to redistribute to other components, and failures in nearby protection devices are more likely. Moreover, if the overall system is more highly loaded, components have smaller margins so they can tolerate smaller increases in load before failure, the system nonlinearities and dynamical couplings increase, and the system operators have fewer options and more stress.

A typical large blackout has an initial disturbance or trigger events, followed by a sequence of cascading events. Each event further weakens and stresses the system and makes subsequent events more likely. Examples of an initial disturbance are short circuits of transmission lines through untrimmed trees, protection device misoperation, and bad weather. The blackout events and interactions are often rare, unusual, or unanticipated because the likely and anticipated failures are already routinely accounted for in power system design and operation. The complexity is such that it can take months after a large blackout to sift through the records, establish the events occurring, and reproduce with computer simulations and hindsight a causal sequence of events.

The historically high reliability of North American power transmission systems is largely due to estimating the transmission system capability and designing and operating the system with margins with respect to a chosen subset of likely and serious contingencies. The analysis is usually either a deterministic analysis of estimated worst cases or a Monte Carlo simulation of moderately detailed probabilistic models that capture steady-state interactions [2]. Combinations of likely contingencies and some dependencies between events such as common mode or common cause are sometimes considered. The analyses address the first few likely failures rather than the propagation of many rare or unanticipated failures in a cascade.

We briefly review some other approaches to cascading failure in power system blackouts. Carreras, Lynch, Dobson, and Newman [4] represented cascading transmission line overloads and outages in a power system model using the DC load flow approximation and standard linear programming optimization of the generation dispatch. The model shows critical point behavior as load is increased and can show power tails similar to those observed in blackout data. Chen and Thorp [9] modeled power system blackouts using the DC load flow approximation and standard linear programming optimization of the generation dispatch and represented in detail hidden failures of the protection system. The expected blackout size is obtained using importance sampling and it shows some indications of a critical point as loading is increased. Rios, Kirschen, Jawayeera, Nedic, and Allan [30] evaluated expected blackout cost using Monte Carlo simulation of a power system model that represents the effects of cascading line overloads, hidden failures of the protection system, power system dynamic instabilities, and the operator responses to these phenomena. Ni, McCalley, Vittal, and Tayyib [26] evaluate expected contingency severities based on real-time predictions of the power system state to quantify the risk of operational conditions. The computations account for current and voltage limits, cascading line overloads, and voltage instability. Roy, Asavathiratham, Lesieutre, and Verghese [31] constructed randomly generated tree networks that abstractly represent influences between idealized components. Components can be failed or operational according to a Markov model that represents both internal component failure and repair processes and influences between components that cause failure propagation. The effects of the network degree and the intercomponent influences on the failure size and duration were studied. Pepyne, Panayiotou, Cassandras, and Ho [29] also used a Markov model for discrete state power system nodal components, but they propagated failures along the transmission lines of a power systems network with a fixed probability. They studied the effect of the propagation probability and maintenance policies that reduce the probability of hidden failures. The challenging problem of determining cascading failure due to dynamic transients in hybrid nonlinear differential equation models was addressed by DeMarco [15] using Lyapunov methods applied to a smoothed model and by Parrilo, Lall, Paganini, Verghese, Lesieutre, and Marsden [28] using Karhunen–Loeve and Galerkin model reduction. Watts [34] described a general model of cascading failure in which failures propagate through the edges of a random network. Network nodes have a random threshold and fail when this threshold is exceeded by a sufficient fraction of failed nodes one edge away. Phase transitions causing large cascades can occur when the network becomes critically connected by having sufficiently average degree or when a highly connected network has sufficiently low average degree so that the effect of a single failure is not swamped by a high connectivity to unfailed nodes. Lindley and Singpurwalla [24] described some foundations for causal and cascading failure in infrastructures and model cascading failure as an increase in a component failure rate within a time interval after another component fails. Initial versions of the cascading failure model of this article appear in Dobson, Chen, Thorp, Carreras, and Newman [18] and Dobson, Carreras, and Newman [16].

3. DESCRIPTION OF MODEL

The model has n identical components with random initial loads. For each component, the minimum initial load is L^min and the maximum initial load is L^max. For j = 1,2,…,n, component j has initial load L_j that is a random variable uniformly distributed in [L^min,L^max]. L₁,L₂,…,L_n are independent.

Components fail when their load exceeds L^fail. When a component fails, a fixed and positive amount of load P is transferred to each of the components.

To start the cascade, an initial disturbance loads each component by an additional amount D. Some components may then fail depending on their initial loads L_j, and the failure of each of these components will distribute an additional load P that can cause further failures in a cascade. The components become progressively more loaded as the cascade proceeds.

In particular, the model produces failures in stages i = 0,1,2,… according to the following algorithm, where M_i is the number of failures in stage i.

CASCADE Algorithm

The CASCADE algorithm has the property that if there are no failures in stage j so that M_j = 0, then 0 = M_j = M_j+1 =··· so that there are no subsequent failures (in step 2, M_j can be zero either because all the components have already failed or because the loads of the unfailed components are less than L^fail). Since there are n components, it follows that M_n = 0 and that the outcome with the maximum number of stages with nonzero failures is 1 = M₀ = M₁ =···= M_n−1. We are most interested in the total number of failures S = M₀ + M₁ +···+ M_n−1.

When the model in an application is being interpreted, the load increment P need not correspond only to transfer of a physical load such as the power flow through a component. Many ways by which a component failure makes the failure of other components more likely can be thought of as increasing an abstract “load” on the other components until failure occurs when a threshold is reached.

It is useful to normalize the loads and model parameters so that the initial loads lie in [0,1] and L^fail = 1 while preserving the sequence of component failures and M₀,M₁,…. First, note that the sequence of component failures and M₀,M₁,… are unchanged by adding the same constant to the initial disturbance D and the failure load L^fail. In particular, choosing the constant to be L^max − L^fail, the initial disturbance D is modified to D + (L^max − L^fail) and the failure load L^fail is modified to L^fail + (L^max − L^fail) = L^max. Then all of the loads are shifted and scaled to yield normalized parameters. The normalized initial load on component j is [ell ]_j = (L_j − L^min)/(L^max − L^min) so that [ell ]_j is a random variable uniformly distributed on [0,1]. The normalized minimum initial load is zero, and the normalized maximum initial load and the normalized failure load are both one. The normalized modified initial disturbance and the normalized load increase when a component fails are

An alternative way to describe the model follows. It is convenient to use the normalized parameters in Eq. (1). Let N(t) be the number of components with loads in (1 − t,1] . If the n initial component loadings are regarded as n points in

, then N(t) is the number of points greater than 1 − t. Then 0 ≤ N(t) ≤ n, the sample paths of N are nondecreasing, and N(t) = 0 for t ≤ 0 and N(t) = n for t ≥ 1.

Let the number of components failed at or before stage j be S_j = M₀ + M₁ +···+ M_j. Then, assuming S₋₁ = 0, the CASCADE algorithm generates S₀,S₁,… according to

Then 0 ≤ S_j ≤ n, S_j is nondecreasing, and S_k = S_k+1 implies that S_j = S_j+1 for j ≥ k. The minimum such k is the maximum stage number in which failures occur and S₋₁ < S₀ < S₁ <···< S_k = S_k+1 =··· and the total number of failures S = S_k; that is,

Moreover, for j < k and r = 0,1,…,M_j+1 − 1,

Therefore, N(d + sp) > s for s = 0,1,…,S − 1, and this inequality and Eq. (3) allow the total number of failures to be characterized as

If, at stage j, d + S_j p > 1, we say that the model saturates. Saturation implies S_j+1 = n. Saturation never occurs if d and p are small enough that d + np < 1.

The model can be formulated as a queue with a single server. Exactly n customers arrive during a given hour independently and uniformly. The server is available to serve these customers at time d after the start of the hour because of completing some other task. The customer service time is p. Then, S is the number of customers that arrive during the first busy period. The queue saturates when the first busy period runs past the end of the hour. Charalambides [10] and Takács [32] analyzed this queue in the nonsaturating case described in Section 4.3.

The model can also be recast in the form of an approximate and idealized fiber bundle model. There are n identical, parallel fibers in the bundle. The L_j of the unnormalized model now indicates breaking strength: Fiber j has random breaking strength L^fail − L_j that is uniformly distributed in [L^fail − L^max,L^fail − L^min]. Each fiber has zero load initially. Then, an initial force is applied to the bundle that increases the load of each fiber to D and this starts a burst avalanche of fiber breaks of size S. When a fiber breaks, it distributes a constant amount of load P to all the other fibers. In contrast, and with better physical justification, idealized fiber bundle models with global redistribution as described by Kloster, Hansen, and Hemmer [22] redistribute the current fiber load equally to the remaining fibers.

4. DISTRIBUTION OF NUMBER OF FAILURES

The main result is that the distribution of the total number of component failures S is

where p ≥ 0 and the saturation function is

It is convenient to assume that 0⁰ ≡ 1 and 0/0 ≡ 1 when these expressions arise in any formula in this article.

If d ≥ 0 and d + np ≤ 1, then there is no saturation (φ(x) = x) and Eq. (7) reduces to the quasibinomial distribution

The quasibinomial distribution was introduced by Consul [11] to model an urn problem in which a player makes strategic decisions. Burtin [3] derived the distribution of the number of initially uninfected nodes that become infected in an inverse epidemic process in a random mapping. This distribution is quasibinomial, with d the fraction of initially infected nodes and p the uniform random mapping probability. Islam et al. [19] interpreted d and p as primary and secondary infection probabilities and applied the quasibinomial distribution to data on the final size of influenza epidemics. Jaworski [20] generalized the derivation to a random mapping with a general fixed-point probability.

The cascading failure model gives a new application and interpretation of the quasibinomial distribution. Moreover, the saturation in Eq. (7) extends the range of parameters of the quasibinomial distribution to allow d + np > 1. Section 5 shows that this extended parameter range can describe regimes with a high probability of all components failing.

The next two subsections derive Eq. (7) from the CASCADE algorithm in two ways: by means of a recursion and by means of the quasimultinomial joint distribution of M₀,M₁,…,M_n−1.

4.1. Recursion

It is convenient to show the dependence of the distribution of number of failures on the normalized parameters by writing P[S = r] = f (r,d,p,n).

In the case of n = 0 components,

According to the CASCADE algorithm, when the initial disturbance d ≤ 0, no components fail, and when d ≥ 1, all n components fail. Then

We assume n > 0 and 0 < d < 1 for the rest of the subsection.

The initial disturbance d causes stage 0 failure of the components that have initial load [ell ] in (1 − d,1] . Therefore, the probability of any component failing in stage 0 is d and

Suppose that M₀ = k and consider the n − k components that did not fail in stage 0. Since none of the n − k components failed in stage 0, their initial loads [ell ] must lie in [0,1 − d] and the distribution of their initial loads conditioned on not failing in stage 0 is uniform in [0,1 − d] . In stage 1, each of the n − k components has had a load increase d from the initial disturbance and an additional load increase kp from the stage 0 failure of k components. Therefore, the equivalent total initial disturbance for each of the n − k components is D = kp + d.

To summarize, assuming M₀ = k, the failure of the n − k components in stage 1 is governed by the model with initial disturbance D = kp + d, load transfer P = p, L^min = 0, L^max = 1 − d, L^fail = 1, and n − k components. Normalizing the parameters using Eq. (1) yields that the failure of the n − k components is governed by the model with normalized initial disturbance kp/(1 − d) and normalized load transfer p/(1 − d); that is,

Combining Eqs. (12) and (13) yields the recursion

Equations (10), (11), and (14) define f (r,d,p,n) = P[S = r] for all n ≥ 0 and p ≥ 0. Equations (10) and (11) agree with Eq. (7). Moreover, it is routine to prove in the Appendix that Eq. (7) satisfies recursion (14). Therefore, Eq. (7) is the distribution of S in the CASCADE algorithm. Thus, the recursion offers a simple way to derive the saturating quasibinomial distribution that avoids complicated algebra or combinatorics. It is also straightforward to use Eqs. (10) and (14) to confirm by induction on n that Eq. (7) is a probability distribution.

4.2. A Quasimultinomial Distribution

This subsection shows that the joint distribution of M₀,M₁,…,M_n−1 is quasimultinomial and hence derives Eq. (7). It is convenient throughout to assume d ≥ 0, restrict m₀, m₁, … to nonnegative integers, and write s_i = m₀ + m₁ +···+ m_i for i = 0,1,… and s₋₁ = 0.

Let α₀ = φ(d), β₀ = 1, and, for i = 1,2,…,

The identity

can be verified using 1 − φ(x) = φ(1 − x) and d ≥ 0 and considering all of the cases.

In step 2 of stage 0 in the CASCADE algorithm, the probability that the load increment of d causes one of the components to fail is α₀ = φ(d) and the probability of m₀ failures in the n components is

Consider the end of step 2 of stage i ≥ 1 in the CASCADE algorithm. The failures that have occurred are M₀ = m₀,M₁ = m₁,…,M_i = m_i and there are n − s_i unfailed components, but the component loads have not yet been incremented by m_i p in step 3.

Suppose that d + s_i−1 p < 1. Then, conditioned on the n − s_i components not yet having failed, the loads of the n − s_i unfailed components are uniformly distributed in [d + s_i−1 p,1] . In step 3, the probability that the load increment of m_i p causes one of the unfailed components to fail is α_i+1 and the probability of m_i+1 failures in the n − s_i unfailed components is

Suppose that d + s_i−1 p ≥ 1. Then, all of the components must have failed on a previous step and P[M_i+1 = m_i+1|M_i = m_i,…,M₀ = m₀] = 1 for m_i+1 = 0 and is zero otherwise. In this case, α_i+1 = 0 and Eq. (18) is verified.

We claim that for s_i ≤ n,

Equation (19) is proved by induction on i. For i = 0, Eq. (19) reduces to Eq. (17). The inductive step is verified by multiplying Eqs. (18) and (19) and using Eq. (16) to obtain P[M_i+1 = m_i+1,…,M₀ = m₀] in the form of Eq. (19).

An expression equivalent to Eq. (19) obtained using Eq. (16) is

The CASCADE algorithm has the property that if there are no failures in stage j so that M_j = 0, then 0 = M_j = M_j+1 =··· and there are no subsequent failures. This property is verified by Eq. (20) because m_j = 0 implies β_j+1 = β_j+2 so that the factor (β_j+1 − β_j+2)^m_j+1 = 0^m_j+1, which vanishes unless m_j+1 = 0. Iterating this argument gives 0 = M_j = M_j+1 =···. Since the maximum number of failures is n, the longest sequence of failures has n stages with M₀ = M₁ =···= M_n−1 = 1. It follows that 0 = M_n = M_n+1 =··· and that the nontrivial part of the joint distribution is determined by M₀,M₁,….,M_n−1. It also follows that M_n−1 = 0 if there are less than n stages with failures.

Equation (20) can now be rewritten for i = n − 1. Let I be the largest integer not exceeding n such that 1 − d − s_I−2 p > 0. Then, Eq. (20) becomes, for s_n−1 ≤ n,

where A(m,n) = 1 and A(m,I) = 0^m_I+1···0^m_n−10^n−s_n−1 for I < n. It follows from the definition of A(m,I) that Eq. (21) vanishes for I < n unless 0 = M_I+1 =···= M_n−1 and S = M₀ +···+ M_I = n. (Although Eq. (21) was derived assuming d ≥ 0, it also holds for d < 0. In particular, for d < 0, Eq. (21) implies P[M_n−1 = 0,…,M₀ = 0] = 1.)

Equation (21) generalizes the quasibinomial distribution and is a form of quasimultinomial distribution. It is a different generalization of the quasibinomial distribution than the quasitrinomial distribution considered by Berg and Mutafchiev [1] to describe numbers of nodes in central components of random mappings.

Suppose that S = M₀ +···+ M_n−1 = r < n. Then, M_n−1 = 0 and M₀ +···+ M_n−2 = r − M_n−1 = r, and Eq. (21) vanishes unless I = n. Summing Eq. (21) over nonnegative integers m₀,…,m_n−1 that sum to r yields

which reduces to Eq. (7) using a lemma by Katz [21]. (The context of Katz's lemma assumes φ(d)/p is a positive integer, but the generalization is immediate.)

4.3. Applying a Generalized Ballot Theorem

Charalambides [10] explained how the quasibinomial distribution appears as a consequence of generalized ballot theorems in the theory of fluctuations of stochastic processes [32]. We summarize this approach and comment that it derives only the nonsaturating cases of Eq. (7).

We assume 0 < d < 1. Consider p multiplied by the number of components N(t) with loads in (1 − t,1] . For 0 ≤ t ≤ 1, pN(t) is a stochastic process with interchangeable increments whose sample functions are nondecreasing step functions with pN(0) = 0. According to Eq. (6), the first passage time of t − pN(t) through d is min{t|pN(t) = t − d} = min{d + sp|N(d + sp) = s} = d + Sp. Then, according to Takács [32, Sect.17, Thm.4],

for 0 < d ≤ t ≤ 1; that is,

Setting t = d + rp in Eq. (23) for r = 0,1,…,min{n,(1 − d)/p}, differencing the resulting equations, and using the binomial distribution of N(t) for 0 ≤ t ≤ 1 yields the nonsaturating cases of Eq. (7). However, the approach does not extend to the saturating cases because pN(t) does not have interchangeable increments when t > 1.

4.4. Approximate Power Tail Exponent at a Critical Case

We describe standard approximations of the quasibinomial distribution that yield a power tail exponent at the critical case. For parameters satisfying np + d ≤ 1 (no saturation), the distribution of S is quasibinomial and can be approximated by letting n → ∞, p → 0, and d → 0 in such a way that λ = np and θ = nd are fixed to give the generalized (or Lagrangian) Poisson distribution [12,13,14]

which is the distribution of the number of offspring in a Galton–Watson–Bienaymé branching process, with the first generation produced by a Poisson distribution with parameter θ and subsequent generations produced by a Poisson distribution with parameter λ. The critical case for the branching process is np = λ = 1 and Otter [27] proved that at criticality, the distribution of the number of offspring has a power tail with exponent −1.5. Further implications for cascading failure of the branching process approximation are considered in Dobson, Carreras, and Newman [17].

5. EFFECT OF LOADING

How much can an electric power transmission system be loaded before there is undue risk of cascading failure? This section discusses qualitative effects of loading on the distribution of blackout size and then applies the model to describe the effect of loading and illustrate its use.

5.2. Effect of Loading in the Model

This subsection describes one way to represent a load increase in the model and how this leads to a parameterization of the normalized model. Then the effect of the load increase on the distribution of the number of components failed is described.

For purposes of illustration, the system has n = 1000 components. Suppose that the system is operated so that the initial component loadings vary from L^min to L^max = L^fail = 1. Then the average initial component loading L = (L^min + 1)/2 may be increased by increasing L^min. The initial disturbance D = 0.0004 is assumed to be the same as the load transfer amount P = 0.0004. These modeling choices for component load lead, via the normalization of Eq. (1), to the parameterization p = d = 0.0004/(2 − 2L), 0.5 ≤ L < 1. The increase in the normalized power transfer p with increased L can be thought of as strengthening the component interactions that cause cascading failure.

The probability distribution of the number S of components failed as L increases from 0.6 is shown in Figure 1. The distribution for the nonsaturating case L = 0.6 has a tail that is approximately exponential. The tail becomes heavier as L increases, and the distribution for the critical case L = 0.8, np = 1 has an approximate power-law region over a range of S. The power-law region has an exponent of approximately −1.4 and this compares to the exponent of −1.5 obtained by the analytic approximation in Section 4.4. The distribution for the saturated case L = 0.9 has an approximately exponential tail for small r, zero probability of intermediate r, and a probability of 0.80 of all 1000 components failing. If an intermediate number of components fail in a saturated case, then the cascade always proceeds to all 1000 components failing.

Log-log plot of distribution of number of components failed S for three values of average initial load L. Note the power-law region for the critical loading L = 0.8. L = 0.9 has an isolated point at (1000,0.80), indicating probability 0.80 of all 1000 components failed. The probability of no failures is 0.61 for L = 0.6, 0.37 for L = 0.8, and 0.14 for L = 0.9.

The increase in the mean number of failures ES as the average initial component loading L is increased is shown in Figure 2. The sharp change in gradient at the critical loading L = 0.8 corresponds to the saturation of Eq. (7) and the consequent increasing probability of all components failing. Indeed, at L = 0.8, the change in gradient in Figure 2 together with the power-law region in the distribution of S in Figure 1 suggest a type 2 phase transition in the system. If we interpret the number of components failed as corresponding to blackout size, the power-law region is consistent with North American blackout data and blackout simulation results [4,8,18]. In particular, North American blackout data suggest an empirical distribution of blackout size with a power tail with exponent between −1 and −2 [6,7,8]. This power tail indicates a significant risk of large blackouts that is not present when the distribution of blackout sizes has an exponential tail [5].

Mean number of components failed ES as a function of average initial component loading L. Note the change in gradient at the critical loading L = 0.8. There are n = 1000 components and ES becomes 1000 at the highest loadings.

The model results show how system loading can influence the risk of cascading failure. At low loading, there is an approximately exponential tail in the distribution of number of components failed and a low risk of large cascading failure. There is a critical loading at which there is a power-law region in the distribution of number of components failed and a sharp increase in the gradient of the mean number of components failed. As loading is increased past the critical loading, the distribution of number of components failed saturates, there is an increasingly significant probability of all components failing, and there is a significant risk of large cascading failure.