Dynamic Modeling of Cascading Failure in Power Systems
Abstract
The modeling of cascading failure in power systems is difficult because of the many different mechanisms involved; no single model captures all of these mechanisms. Understanding the relative importance of these different mechanisms is important for choosing which mechanisms need to be modeled for particular applications. This work presents a dynamic simulation model of both power networks and protection systems, which can simulate a wider variety of cascading outage mechanisms relative to existing quasisteady state (QSS) models. This paper describes the model and demonstrates how different mechanisms interact. In order to test the model we simulated a batch of randomly selected contingencies for several different static load configurations, and found that the distributions of blackout sizes and event lengths from the simulator correlate well with historical trends. The results also show that load models have significant impacts on the cascading risks. Finally, the dynamic model was compared against a simple dcpowerflow based QSS model; we find that the two models tend to agree for the early stages of cascading, but produce substantially different results for later stages.
IEEEexample:BSTcontrol
I Introduction
The vital significance of studying cascading outages has been recognized in both the power industry and academia [1, 2, 3]. However, since electrical power networks are very large and complex systems [4], understanding the many mechanisms by which cascading outages propagate is challenging. This paper presents the design of and results from a new nonlinear dynamic model of cascading failure in power systems (the Cascading Outage Simulator with Multiprocess Integration Capabilities or COSMIC), which can be used to study a wide variety of different mechanisms of cascading outages.
A variety of cascading failure modeling approaches have been reported in the research literature, many of which are reviewed in [2, 1, 3]. Several have used quasisteady state (QSS) dc power flow models [5, 6, 7], which are numerically robust and can describe cascading overloads. However, they do not capture nonlinear mechanisms like voltage collapse or dynamic instability. QSS ac power flow models have been used to model cascading failures in [8, 9, 10, 11], but these models still require difficult assumptions to model machine dynamics and to deal with nonconvergent power flows. Some have proposed models that combine the dc approximations and dynamic models [12], allowing for more accurate modeling of undervoltage and underfrequency load shedding. These methods increase the modeling fidelity over pure dc models but still neglect voltage collapse. Others developed statistical models that use data from simulations [13, 14] or historical cascades [15] to represent the general features of cascading. Statistical models are useful, but cannot replace detailed simulations to understand particular cascading mechanisms in depth. There are also topological models [16, 17, 18, 19] which have been applied to the identification of vulnerable/critical elements, however, without detailed power grid information, the results that they yield differ greatly and could result in misleading conclusions about the grid vulnerability [20].
Some dynamic models [21, 22] and numerical techniques [23, 24] study the mid/longterm dynamics of power system behavior, and show that mid/long term stability is an important part of cascading outage mechanisms. However, concurrent modeling of power system dynamics and discrete protection events—such as line tripping by overcurrent, distance and temperature relays, undervoltage and underfrequency load shedding—is challenging and not considered in most existing models. In [25] the authors describe an initial approach using a system of differentialalgebraic equations with an additional set of discrete equations to dynamically model cascading failures.
The paper describes the details of and results from a new nonlinear dynamic model of cascading failure in power systems, which we call “Cascading Outage Simulator with Multiprocess Integration Capabilities” (COSMIC). In COSMIC, dynamic components, such as rotating machines, exciters, and governors, are modeled using differential equations. The associated power flows are represented using nonlinear power flow equations. Load voltage responses are explicitly represented, and discrete changes (e.g., components failures, load shedding) are described by a set of equations that indicate the proximity to thresholds that trigger discrete changes. Given dynamic data for a power system and a set of exogenous disturbances that may trigger a cascade, COSMIC uses a recursive process to compute the impact of the triggering event(s) by solving the differentialalgebraic equations (DAEs) while monitoring for discrete events, including events that subdivide the network into islands.
There are only a few existing commercial and researchgrade simulation tools that specifically address cascading failure events and their consequences [2]. Static modeling is still dominant in these tools, although at least two commercial packages (Eurostag ASSESS [26] and POMPCM [27]) have introduced dynamic simulation to their cascading failure analysis. The model that we propose here supplements these commercial tools by providing an open platform for research and development that allows one to explicitly test the impact of the many assumptions that are necessary for dynamic cascading failure modeling. For example, users can modify the existing system components, add new ones, and integrate advanced remedial control actions. Additionally, the dynamic/adaptive time step and recursive islanded time horizons (see Section IIC) implemented in this simulator allow for faster computations during, or near, steadystate regimes, and fine resolution during transient phases. Moreover, this tool can be easily integrated with High Performance Computing (HPC) clusters to run many simulations simultaneously at a much lower cost, relative to commercial tools.
The remainder of this paper proceeds as follows: Section II introduces the components of the model mathematically and describes how different modules interact. In Section III, we present results from several experimental validation studies. Section IV presents our conclusions from this study. Finally, Section V provides further details on the model and its settings.
Ii Hybrid system modeling in COSMIC
Iia Hybrid differentialalgebraic formulation
Dynamic power networks are typically modeled as sets of DAEs [28]. If one also considers the dynamics resulting from discrete changes such as those caused by protective relays, an additional set of discrete equations is added, which results in a hybrid DAE system [29].
Let us assume that the state of the power system at time can be defined by three vectors: , , and , where:

is a vector of continuous state variables that change with time according to a set of differential equations
(1) 
is a vector of continuous state variables that have pure algebraic relationships to other variables in the system:
(2) 
is a vector of state variables that can only take integer states ()
(3)
When constraint fails, an associated counter function (see IIB) activates. Each changes state if reaches its limit.
The set of differential equations (1) represent the machine dynamics (and/or load dynamics if dynamic load models are included). In COSMIC the differential equations include a third order machine model and somewhat simplified governor and exciter models in order to improve computational efficiency without compromising the fundamental functions of those components (see Section VA for more details). In particular, the governor is rate and rail limited to model the practical constraints of generator power control systems. The governor model incorporates both droop control and integral control, which is important to mid/long term stability modeling, especially in isolated systems [28].
The algebraic constraints (2) encapsulate the standard ac power flow equations. In this study we implemented both polar and rectangular power flow formulations. Load models are an important part of the algebraic equations, which are particularly critical components of cascading failure simulation because i) they need to represent the aggregated dynamics of many complicated devices and ii) they can dramatically change system dynamics. The baseline load model in COSMIC is a static model, which can be configured as constant power (), constant current (), constant impedance (), exponential (), or any combination thereof () [30].
As Fig. 1 illustrates, load models can have a dramatic impact on algebraic convergence. Constant power loads are particularly difficult to model for the offnominal condition. Numerical failures are much less common with constant I or Z loads, but are not accurate representations of many loads. This motivated us to include the exponential component in COSMIC.
During cascading failures, power systems undergo many discrete changes that are caused by exogenous events (e.g., manual operations, weather) and endogenous events (e.g., automatic protective relay actions). The discrete event(s) will consequently change algebraic equations and the systems dynamic response, which may result in cascading failures, system islanding, and large blackouts. In COSMIC, the endogenous responses of a power network to stresses are represented by (3). These discrete responses are described in detail in IIB and IIC.
IiB Relay modeling
Major disturbances cause system oscillations as the system seeks a new equilibrium. These oscillations may naturally die out due to the interactions of system inertia, damping, and exciter and governor controls. In order to ensure that relays do not trip due to brief transient state changes, timedelays are added to each protective relay in COSMIC.
We implemented in this model two types of timedelayed triggering algorithms: fixedtime delay and timeinverse delay. These two delay algorithms are modeled by a counter function, , which is triggered by (3). The fixedtime delay triggering activates its counter/timer as soon as the monitored signal exceeds its threshold. If the signal remains beyond the threshold, this timer will continue to count down from a preset value until it runs out then the associated relay takes actions. Similarly, the timer will recover if the signal is within the threshold and will max out at the preset value. For the timeinverse delay algorithm, instead of counting for the increment of time beyond (or within) the threshold, we evaluate the area over (or under) the threshold by integration (based on Euler’s rule).
Five types of protective relays are modeled in COSMIC: overcurrent (OC) relays, distance (DIST) relays, and temperature (TEMP) relays for transmission line protection; as well as undervoltage load shedding (UVLS) and underfrequency load shedding (UFLS) relays for stress mitigation. OC relays monitor the instantaneous current flow along each branch. DIST relays represent a Zone 1 relay that monitors the apparent admittance of the transmission line. TEMP relays monitor the line temperature, which is obtained from a first order differential equation
(4) 
where is the temperature difference relative to the ambient temperature (20 C) for line , is the current flow of line , and are the are heating and time constants for line [12]. and are chosen so that line ’s temperature reaches 75 C (ACSR conductors) if current flow hits the rateA limit, and its TEMP relay triggers in 60 seconds when current flow jumps from rateA to rateC. The threshold for each TEMP relay is obtained from the rateB limit. While it would have been possible to incorporate the temperature relays into the trapezoidal integration used for the other differential equations, these variables are much slower than the other differential variables in . Instead, we computed outside of the primary integration system using Euler’s rule.
When voltage magnitude or frequency signals at load Bus are lower than the specified thresholds, the UVLS relay or UFLS relay will shed a (default setting) of the initial to avoid the onset of voltage instability and reduce system stress. In order to monitor frequency at each load bus, Dijkstra’s algorithm [31] and electrical distances [32] were used to find the generator (and thus frequency from ) that is most proximate to each load bus. Both the UVLS and UFLS relays used a fixedtime delay of seconds.
IiC Solving the hybrid DAE
Because of its numerical stability advantages, COSMIC uses the trapezoidal rule [33] to simultaneously integrate and solve the differential and algebraic equations (see Section VB for more details).
Whereas many of the common tools in the literature [34, 35] use a fixed time stepsize, COSMIC implements a variable time stepsize in order to tradeoff between the diverse timescales of the dynamics that we implement. We select small step sizes during transition periods that have high deviation or oscillations in order to keep the numerical within tolerance; the step sizes increase as the oscillations dampen toward steadystate values.
When a discrete event occurs at , the complete representation of the system is provided in the following equations:
(5)  
(6)  
(7)  
(8) 
where, is the previous time point, and is the counter function mentioned previously. Because of the adaptive stepsize, COSMIC retains from , in which is found by linear interpolation of two time steps.
Every time a discrete event happens ( and ), COSMIC stops solving the DAE for the previous network configuration, processes the discrete event(s), then resumes the DAE solver using the updated initial condition. Specifically, it uses updated algebraic variables while holding the differential variables. In some cases the discrete event causes the system to split into multiple islands. COSMIC checks for this condition by inspection of the network admittance matrix () (see Algorithm 1). COSMIC deals with the separation of a power network into subnetworks (unintentional islanding) using a recursive process that is described in Algorithm 1. If islanding results from a discrete event, the present hybrid DAE separates into two sets of DAEs. COSMIC treats the two subnetworks the same way as the original one, integrates and solves these two DAE systems in parallel, and synchronizes two result sets in the end.
IiD Validation
To validate COSMIC, we compared the dynamic response in COSMIC against commercial software—PowerWorld [35]—using the classic 9bus test case [36]. From a random contingency simulation, the Mean Absolute Error (MAE) between the results produced by COSMIC and PowerWorld was within 0.11%. Since COSMIC adopted simplified exciter and governor models that are not included in many commercial packages, several of the time constants were set to zero or very close to zero in order to obtain agreement between the two models.
Iii Experiments and Results
In this section we present the results from several computational experiments, which illustrate this model and the types of insight that one can gain with dynamic cascading failure simulation. Three test systems are used: the 9bus system [36], the 39bus system, and the 2383bus system [37], which is an equivalenced system based on the year 2000 winter snapshot for the Polish network. Section IIIA compares the computational efficiency of the polar and rectangular formulations of the model. Section IIIB uses the 9bus test case to illustrate the different relay functions and their time delay algorithms. Section IIIC demonstrates how COSMIC processes cascading events such as line branch outages, load shedding, and islanding. Section IIID studies the impact of different load modeling assumptions on cascading failure sizes. Finally, Section IIIE compares results from COSMIC with results from a dcpower flow based model of cascading failure in order to understand similarities and differences between these two different modeling approaches.
Iiia Polar formulation vs. rectangular formulation in computational efficiency
COSMIC includes both polar and rectangular power flow formulations. In order to compare the computational efficiency of the two formulations, we conducted a number of and experiments using two different cases, the 39bus and 2383bus systems. The amount of time that a simulation requires and the number of linear solves () are two measures commonly used to evaluate the computational efficiency of a model. Compared to the first metric, the number of linear solves describes computational speed independently of specific computing hardware, and it is adopted in this study.
Tables II and II compare the two formulations with respect to the number of linear solves. These comparisons are grouped by the amount of demand lost in each simulation, given that longer cascades tend to require more linear solves and have different numerical properties. For the 39bus case, 45 experiments and 222 randomly selected experiments were conducted, and each of them finished at 50 seconds and lost the same amount of power demand. As shown in Table II, the performances of the two methods were similar in terms of number of linear solves; however, the rectangular formulation required fewer linear solves and shows various improvements (e.g., positive decrease rectangular vs. polar) for different demand losses.
For the 2383bus test case, we simulated 2494 and 556 contingencies; Table II shows the results. There was no significant improvement for the rectangular formulation over the polar formulation, and the number of linear solves that resulted from both forms were almost identical.
One can also notice that solving the 2383bus case required fewer linear solves than for the 39bus case. This suggests that some branch outages have a higher impact on a smaller network and cause more dynamic oscillations than on a larger network such as the 2383bus case.
01  1200  200500  5001000  10003000  30006000  
MW  MW  MW  MW  MW  MW  
Polar  3084  6839  10523  5994  41067  10278 
Rec.  3035  6724  10324  5846  4052  9501 
Tests  176  38  14  32  4  3 
% Dec.  1.59%  1.69%  1.89%  2.46%  1.33%  7.6% 
1: The density (nonzero rates) of the Jacobian matrices for polar and rectangular forms are 0.0376% and 0.0383% respectively; 2: the number of tests; 3: The percentage decrease in the number of linear solves, rectangular vs. polar.
01 MW
12000 MW
Polar
867.77
692.81
Rec.
867.84
693.10
Tests
2494
556
% Dec.
0.009%
0.04%
4: The density (nonzero rates) of the Jacobian matrices for polar and rectangular formulations are 0.000856% and 0.000872% respectively.
IiiB Relay event illustration
To depict the functionality of how protective relays integrate with COSMIC’s time delay features, we implemented the following example using the 9bus system. The initial event was a singleline outage from Bus 6 to Bus 9 at seconds. The countdown timer of the DIST relay for branch 57 was activated with a seconds. As shown in Fig. 2, the system underwent a transient swing following the oneline outage. Right after seconds, ran out () and branch 57 was tripped by its DIST relay, which resulted in two isolated islands. Meanwhile, the thresholds for UVLS relays were set to pu. Note that the magenta voltage trace violated this voltage limit at about seconds (). Because this trace continued under the limit after that, its UVLS timer counted down from seconds till seconds (), where its UVLS relay took action and shed of the initial load at this bus. The adjacent yellow trace illustrates that the UVLS relay timer was activated as well but with a small lag. This UVLS relay never got triggered because before its emptied out the load shedding at put this yellow trace back upon threshold, and its was restored to seconds.
IiiC Cascading outage examples using the 39bus and the 2383bus power systems
The following experiment demonstrates a cascading outage example using the IEEE 39bus case (see Table III for a summary of the sequential events). The system suffered a strong dynamic oscillation after the initial two exogenous events (branches 225 and 56). After approximately 55 seconds the first OC relay at branch 45 triggered. Because the monitored current kept upcrossing and downcrossing its limit, it delayed the relay triggering based on the time delay algorithms for the protection devices. Load shedding at two buses (Bus 7 and Bus 8) occurred around seconds, then another two branches (1013 and 1314) shut down after OC relay trips at seconds. These events separated the system into two islands. At seconds, two branches (34 and 1718) were taken off the grid and this resulted in another island. The system eventually ended up with three isolated networks. However, one of them was not algebraically solvable due to a dramatic power imbalance, and it was declared as a blackout area.
No.  Time (sec)  Events 

1  3  Initial events: branches 225 and 56 trip 
2  54.78  Branch 67 is tripped by OC relay 
3  55.06  58.45 MW load shedding at Bus 7 
4  55.07  130.50 MW load shedding at Bus 8 
5  55.28  Branch 414 is tripped by OC relay 
6  55.28  Branch 1013 is tripped by OC relay 
7  55.28  Branch 1314 is tripped by OC relay 
—  55.28  1st islanding event 
8  55.78  Branch 34 is tripped by OC relay 
9  55.78  Branch 1718 is tripped by OC relay 
—  55.78  2nd islanding event 
Fig. 3 illustrates a sequence of branch outage events using the 2383bus network. This cascading was initiated by two exogenous branch outages (Branches 3132 and 388518, which are marked as number 0 in Fig. 3) and resulted in a total of 92 discrete events. Out of these, 24 were branch outage events, and they are labeled in order in Fig. 3. These events consequently caused a small island (light blue colored dots). The dots with additional red square highlighting indicate buses where load shedding occurred. From this figure we can see that cascading sequences do not follow an easily predictable pattern, and the affected buses with low voltage or frequency may be far from the initiating events.
The top panel in Fig. 4 shows the timeline of all branch outage events for the above cascading scenario, and the lower panel zooms in the loadshedding events. In the early phase of this cascading outages, the occurrence of the components failed relatively slowly, however, it speeded up as the number of failures increased. Eventually the system condition was substantially compromised, which caused fast collapse and the majority of the branch outages as well as the load shedding events (see lower panel in Fig. 4).
IiiD contingency analysis using the 2383bus case
Power systems are operated to ensure the security criterion so that any single component failure will not cause subsequent contingencies [38]. The modified 2383bus system that we are studying in this paper satisfies this criterion for transmission line outages. Thus, we assume here that branch outages capture a wide variety of exogenous contingencies that initiate cascades, for example a transformer tripping due to a generator failure.
The experiment implemented here included four groups of 1200 randomly selected contingencies for the 2383bus system. We measured the size of the resulting cascades using the number of relay events and the amount of demand lost. Each group had a different static load configuration. The load configuration for the first group was 100% constant Z load; the second group used 100% E load; the third had 100% constant P load; and the fourth one included 25% of each portion in the ZIPE model. We set TEMP, DIST, UVLS, and UFLS relays active in this experiment and deactivated OC relay due to its overlapping/similar function with TEMP relay.
Fig. 5 shows the Complementary Cumulative Distribution Function (CCDF) of demand losses for these four groups of simulations. The CCDF plots of demand losses exhibit a heavytailed blackout size distribution, which are typically found in both historical blackout data and cascading failure models [39]. The magenta trace indicates constant Z load, and shows the best performance –in terms of the average power loss and the probability of large blackout– within this set of 1200 random contingencies (listed in Table IV). In contrast, the blue trace (constant E load) reveals the highest risk of large size blackouts ( MW). The constant P load has a similar trend as the constant E load, due to their similar stiff characteristics; however, the constant E load with this particular exponent, 0.08, demonstrates a negative effect on the loss of load. The one with 25% Z, 25% I, 25% P and 25% of E performs in the middle of constant P load and constant Z load.
As can be seen in Table IV, the probabilities of large demand losses varies from 2.5% to 3.5% for those four load configurations. These results show that load models play an important role in dynamic simulation and may increase the frequency of nonconvergence if they are not properly modeled.
Load Model  Avg. Loss  Avg. BO  Prob. of Largest Loss 

644.19 MW  0.2092  0.025  
889.02 MW  0.1800  0.035  
848.21 MW  0.1775  0.033  
807.11 MW  0.2042  0.032 
1: BO — branch outages.
Fig. 7 shows the CCDF plots of total dependent event length, including all event types after the initial contingencies, such as branch outages caused by TEMP and DIST relays, and load shedding events by UVLS and UFLS relays. Fig. 8 shows the CCDF of the branch outage lengths only. We can see from these two figures that the distributions of constant P and constant E loads have a comparable pattern, and they are in general less likely to have the same amount of branch outages, relative to the other two configurations. In Fig. 6 we compare the amount of load shedding in the 88 simulated cascades to the number of discrete events (including the two exogenous line outages), the total number of events, and the cascading time. Each of these 88 samples includes at least one dependent event (31 of the 1200 simulations did not converge and were declared as full blackouts, and 1081 of them were resilient cases without any dependent events). The results show that there is a positive, but weak, correlation between the number of events and the blackout size in MW. We also find that there is a positive, but weak, correlation between cascading failure size in MW and the length of the cascade in seconds.
IiiE Comparison with a dc cascading outage simulator
A number of authors have implemented quasisteady state (QSS) models using the dc power flow equations to investigate cascading outages [5, 6, 7]. The dc model includes numerous simplifications that are substantially different from the “real” system. However, the dc model is numerically stable, making it possible to produce results that can be statistically similar to data from real power systems [40]. On the other hand, dynamic models, such as the one presented here, include many mechanisms of cascading that cannot be represented in the dc model. In order to understand the implications of these differences, we conducted two experiments to compare COSMIC with the dc QSS model described in [5] with respect to the overall probabilities of demand losses and the extent to which the patterns of cascading from the two models agree. These experiments provide helpful insights into the dominant cascading mechanisms in different phases of cascading sequences, and into the conditions under which one would want to use a dynamic model as opposed to a simpler model.
IiiE1 The probabilities of demand losses
For the first experiment, we computed the CCDF of demand losses in both COSMIC (with the constant impedance load model) and the dc model using the same 1200 branch outage pairs from IIID. From Fig. 9 one can learn that the probability of demand losses in the dc simulator is lower than that of COSMIC for the same amount of demand losses. In particular, the largest demand loss in dc simulator is much smaller than in COSMIC (2639 MW vs. 24602 MW, with probabilities 0.08% vs. 2.5%). This large difference between them is not surprising because the dc model is much more stable and does not run into problems of numerical nonconvergence. Also, the protection algorithms differ somewhat between the two models. In addition, some of the contingencies do produce large blackouts in the dc simulator, which causes the fat tail that can be seen in Fig. 9.
Numerical failures in solving the DAE system greatly contributed to the larger blackout sizes observed in COSMIC, because COSMIC assumes that the network or subnetwork in which the numerical failure occurred experienced a complete blackout. This illustrates a tradeoff that comes with using detailed nonlinear dynamic models: while the component models are more accurate, the many assumptions that are needed substantially impact the outcomes, potentially in ways that are not fully accurate.
It is possible that this result may be affected to some extent by the size of the sample set. The 1200 randomly selected contingency pairs represent 0.0278% of the total branch outage pairs. Extensive investigation of how different sampling approaches might impact the observed statistics remains for future work.
336 critical pairs  
COSMIC  dc simulator  
Avg. Demand Loss  23873 MW  3039.5 MW 
Avg. Branch Outage Length  7.5476  26.7351 
Avg. R  0.1948  
Std. Dev. of R  0.1381  
Max R  0.6  
Avg. R (first 10)  0.3487  
Std. Dev. of R (first 10)  0.2445  
Max R (first 10)  1 
1: only the first 10 branch outages.
IiiE2 Path Agreement Measurement
The second comparison experiment was to study the patterns of cascading between these two models. This comparison provides additional insight into the impact of dynamics on cascade propagation patterns. In order to do so, we compared the sets of transmission lines that failed using the “Path Agreement Measure” introduced in (9) [12]. The relative agreement of cascade paths, , is defined as follows. If models and are both subjected to the same set of exogenous contingencies: . measures the average agreement in the set of dependent events that result from each contingency in each model. If contingency results in the set of dependent branch failures in model and the set of dependent branch failures in , is defined as:
(9) 
The experiment measured between COSMIC (with a load configuration) and the dc simulator for 336 critical branch outage pairs. These branch outage pairs are the ones reported in [5] to produce large cascading blackouts. Table V shows that the average between the two models for the whole set of sequences is , which is relatively low. This indicates that there are substantial differences between cascade paths in the two models. Part of the reason is that the dc model tends to produce longer cascades and consequently increase the denominator in (9). In order to control this, we computed only for the first 10 branch outage events. The average R increased to 0.3487, and some of the cascading paths showed a perfect match (). This suggests that the cascading paths resulting from the two models tend to agree during the early stages of cascading, when nonlinear dynamics are less pronounced, but disagree during later stages.
Iv Conclusions
This paper describes a method for and results from simulating cascading failures in power systems using full nonlinear dynamic models. The new model, COSMIC, represents a power system as a set of hybrid discrete/continuous differential algebraic equations, simultaneously simulating protection systems and machine dynamics. Several experiments illustrated the various components of COSMIC and provided important general insights regarding the modeling of cascading failure in power systems. By simulating 1200 randomly chosen contingencies for a 2383bus test case, we found that COSMIC produces heavytailed blackout size distributions, which are typically found in both historical blackout data and cascading failure models [39]. However, the relative frequency of very large events may be exaggerated in dynamic model due to numerical nonconvergence (about of cases). More importantly, the blackout size results show that load models can substantially impact cascade sizes—cases that used constant impedance loads showed consistently smaller blackouts, relative to constant current, power or exponential models. In addition, the contingency simulation results from COSMIC were compared to corresponding simulations from a dc power flow based quasisteadystate cascading failure simulator, using a new metric. The two models largely agreed for the initial periods of cascading (for about 10 events), then diverged for later stages where dynamic phenomena drive the sequence of events.
Together these results illustrate that detailed dynamic models of cascading failure can be useful in understanding the relative importance of various features of these models. The particular model used in this paper, COSMIC, is likely too slow for many largescale statistical analyses, but comparing detailed models to simpler ones can be helpful in understanding the relative importance of various modeling assumptions that are necessary to understand complicated phenomena such as cascading.
V Appendix
Section VA presents the detailed dynamic component models in COSMIC, to further explain the dynamic equations (in Section II). Section VB explains the formulation of the numerical solver in this model to solve the continuous integration. Section VC shows the relevant settings for the numerical solver and protective relays in COSMIC. This is free open source code with a GNU General Public License and is available for download and contributions in a web repository [41].
Va Differential equations in dynamic power system
The following differential equations are used to represent machine dynamics in COSMIC, based on the polar formulation.
VA1 Equation for rotor speed — Swing Equation
The rotor speed related to bus is represented in the standard second order swing equation [28], which describes the rotational dynamics of a synchronous machine. The normalized rotor speed is fixed during the normal operation, and will accelerate or decelerate when there’s disturbance.
(10) 
where ( is the set of all generator buses), is a machine inertia constant, is a damping constant, is the mechanical power input, and is the generator power output.
VA2 Equation for rotor angle
The rotor angle, , is the integral of the relative rotor speed change with respect to synchronous speed (1.0 in per unit notation). It is given by the equation:
(11) 
VA3 Generator
The salientpole model is adapted in COSMIC. The active and reactive power outputs are given by the nonlinear equations [42]:
(12) 
(13) 
where , and are the direct axis generator synchronous and transient reactances, respectively. The transient open circuit voltage magnitude [42], , is determined by the differential equation:
(14) 
where is the direct axis transient time constant, and is the machine exciter output. This equation, along with Eqs. 1011, describe the basic physical properties of the generation machine, which results in a third order differential equation system.
VA4 Exciter
The machine exciter control in COSMIC utilizes a generic model, which is a second order differential system:
(15) 
(16) 
where is the desired reference voltage, is the actual terminal voltage, , and are the exciter time constants, and sigm is a differentiable sigmoidal function that acts as a limiter between and . Fig. 10 illustrates the simplified exciter configuration.
VA5 Governor
The machine governor controls the mechanical forcing given a deviation in the machine speed. COSMIC adapts a governor model that includes both proportional (droop) control and an approximate representation of secondary (integral) control in order to ensure that the system can restore frequency deviations to zero. It can be described by Eqs. (17)(18):
(17) 
(18) 
where , and are the droop, PI time constants and servomotor time constants, respectively. is a intermediate differential variable in the PI controller. Fig. 11 illustrates the interactions among the governor variables.
VB Trapezoidal method to solve DAE system
Let be the current time step, and be the next time step. Assuming that we have a consistent set of variables ,and for the current time , and can be calculated from them. The trapezoidal solution method solves the following nonlinear system to obtain and :
(19) 
(20) 
VC The relevant settings in COSMIC for this study
Table VI and Table VII list the parameters that are used in the numerical solution, as well as the protective relay settings in COSMIC, respectively.
Item  Setting 1  Setting 2 

Convergence Tolerance  —  
Maximum Number of Iterations  20  — 
Settings to Increase  0.01 [mismatch]  
Settings to Decrease  0.05 [mismatch]  
Maximum  1 sec.  — 
Minimum  0.005 sec.  — 
1: If , increase to , except when a discrete event occurs. 2: If , reduce to , except when a discrete event occurs.
Activation Threshold  Delay  Other  

OC relays  rateB /  Timeinverse delay  sec. 
DIST relays  90% of the line impedance  0.5 sec.  — 
TEMP relays  temperature converted from rateB  0 sec.  — 
UVLS relays  0.87 pu  0.5 sec.  25% 
UFLS relays  0.95 pu  0.5 sec.  25% 
3: The triggering threshold for the OC relays. 4: Load shedding factor, indicating the percentage of shed power over total load for each load shedding action.
References
 [1] R. Baldick et al., “Initial review of methods for cascading failure analysis in electric power transmission systems: IEEE PES CAMS task force on understanding, prediction, mitigation and restoration of cascading failures,” in Power and Energy Society General Meeting, 2008 IEEE, 2008, pp. 1–8.
 [2] M. Papic et al., “Survey of tools for risk assessment of cascading outages,” in Power and Energy Society General Meeting, 2011 IEEE, 2011, pp. 1–9.
 [3] M. Vaiman et al., “Risk assessment of cascading outages: Methodologies and challenges,” IEEE Trans. Power Syst., vol. 27, no. 2, pp. 631–641, 2012.
 [4] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, “Complex systems analysis of series of blackouts: Cascading failure, critical points, and selforganization,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 17, no. 2, p. 026103, 2007.
 [5] M. Eppstein and P. Hines, “A ‘random chemistry’ algorithm for identifying collections of multiple contingencies that initiate cascading failure,” IEEE Trans. Power Syst., vol. 27, no. 3, pp. 1698–1705, 2012.
 [6] B. A. Carreras, V. E. Lynch, I. Dobson, and D. E. Newman, “Critical points and transitions in an electric power transmission model for cascading failure blackouts,” Chaos: An interdisciplinary journal of nonlinear science, vol. 12, no. 4, pp. 985–994, 2002.
 [7] J. Yan, Y. Tang, H. He, and Y. Sun, “Cascading failure analysis with DC power flow model and transient stability analysis,” IEEE Trans. Power Syst., vol. PP, no. 99, pp. 1–13, 2014.
 [8] D. P. Nedic, I. Dobson, D. S. Kirschen, B. A. Carreras, and V. E. Lynch, “Criticality in a cascading failure blackout model,” International Journal of Electrical Power & Energy Systems, vol. 28, no. 9, pp. 627 – 633, 2006.
 [9] Q. Chen and L. Mili, “Composite power system vulnerability evaluation to cascading failures using importance sampling and antithetic variates,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2321–2330, Aug 2013.
 [10] S. Mei, Y. Ni, G. Wang, and S. Wu, “A study of selforganized criticality of power system under cascading failures based on ACOPF with voltage stability margin,” IEEE Trans. Power Syst., vol. 23, no. 4, pp. 1719–1726, Nov 2008.
 [11] D. Bienstock, “Adaptive online control of cascading blackouts,” in Power and Energy Society General Meeting, 2011. IEEE, 2011, pp. 1–8.
 [12] R. Fitzmaurice, E. CotillaSanchez, and P. Hines, “Evaluating the impact of modeling assumptions for cascading failure simulation,” in Power and Energy Society General Meeting, 2012 IEEE, July 2012, pp. 1–8.
 [13] M. RahnamayNaeini, Z. Wang, N. Ghani, A. Mammoli, and M. Hayat, “Stochastic analysis of cascadingfailure dynamics in power grids,” IEEE Trans. Power Syst., vol. 29, no. 4, pp. 1767–1779, July 2014.
 [14] P. D. Hines, I. Dobson, E. CotillaSanchez, and M. Eppstein, “Dual graph and random chemistry methods for cascading failure analysis,” in 2013 46th Hawaii International Conference on System Sciences (HICSS), Jan 2013, pp. 2141–2150.
 [15] I. Dobson, “Estimating the propagation and extent of cascading line outages from utility data with a branching process,” IEEE Trans. Power Syst., vol. 27, no. 4, pp. 2146–2155, Nov 2012.
 [16] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.U. Hwang, “Complex networks: Structure and dynamics,” Physics reports, vol. 424, no. 4, pp. 175–308, 2006.
 [17] P. Hines, E. CotillaSanchez, and S. Blumsack, “Topological models and critical slowing down: Two approaches to power system blackout risk analysis,” in 2011 44th Hawaii International Conference on System Sciences (HICSS). IEEE, 2011, pp. 1–10.
 [18] J.W. Wang and L.L. Rong, “Cascadebased attack vulnerability on the US power grid,” Safety Science, vol. 47, no. 10, pp. 1332–1336, 2009.
 [19] R. Albert, I. Albert, and G. L. Nakarado, “Structural vulnerability of the North American power grid,” Physical review E, vol. 69, no. 2, p. 025103, 2004.
 [20] P. Hines, E. CotillaSanchez, and S. Blumsack, “Do topological models provide good information about electricity infrastructure vulnerability?” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 20, no. 3, p. 033122, 2010.
 [21] H. T. Ma, M. L. Crow, B. H. Chowdhury, and A. Lininger, “Cascading line outage prevention with multiple UPFCs,” in 2007 39th North American Power Symposium (NAPS’07). IEEE, 2007, pp. 273–278.
 [22] D. Fabozzi and T. Van Cutsem, “Simplified timedomain simulation of detailed longterm dynamic models,” in Power Energy Society General Meeting, 2009. PES’09. IEEE, 2009, pp. 1–8.
 [23] S. K. Khaitan, C. Fu, and J. McCalley, “Fast parallelized algorithms for online extendedterm dynamic cascading analysis,” in Power Systems Conference and Exposition, 2009. PSCE’09. IEEE, 2009, pp. 1–7.
 [24] S. Abhyankar and A. J. Flueck, “Realtime power system dynamics simulation using a parallel BlockJacobi preconditioned NewtonGMRES scheme,” in 2012 SC Companion: High Performance Computing, Networking, Storage and Analysis (SCC). IEEE, 2012, pp. 299–305.
 [25] C. Parmer, E. CotillaSanchez, H. K. Thornquist, and P. D. Hines, “Developing a dynamic model of cascading failure for high performance computing using Trilinos,” in Proc. of the 1st Int’l Workshop on High Performance Computing, Networking and Analytics for the Power Grid, ser. HiPCNAPG ’11. New York, NY, USA: ACM, 2011, pp. 25–34.
 [26] J. Antoine and M. Stubbe, “Eurostag, software for the simulation of power system dynamics. its application to the study of a voltage collapse scenario,” in IEEE Colloquium on Interactive Graphic Power System Analysis Programs, Mar 1992, pp. 5/1–5/4.
 [27] N. Bhatt et al., “Assessing vulnerability to cascading outages,” in Power Systems Conference and Exposition, 2009. PSCE ’09. IEEE/PES, March 2009, pp. 1–9.
 [28] P. Kundur, Power System Stability and Control. Tata McGrawHill Education, 1994.
 [29] A. J. Van Der Schaft, J. M. Schumacher, A. J. van der Schaft, and A. J. van der Schaft, An Introduction to Hybrid Dynamical Systems. Springer London, 2000, vol. 251.
 [30] J. Song, E. CotillaSanchez, and T. Brekken, “Load modeling methodologies for cascading outage simulation considering power system stability,” in 2013 1st IEEE Conference on Technologies for Sustainability (SusTech), 2013, pp. 78–85.
 [31] Dijkstra’s Algorithm. [Online]. Available: http://en.wikipedia.org/wiki/Dijkstra’s_algorithm
 [32] E. CotillaSanchez, P. Hines, C. Barrows, and S. Blumsack, “Comparing the topological and electrical structure of the North American electric power infrastructure,” IEEE Systems Journal, vol. 6, no. 4, pp. 616–626, 2012.
 [33] F. Milano, Power System Modeling and Scripting. Springer, 2010.
 [34] SIEMENS, PTI, “PSS/E,” E 30.2 Program Operational Manual, 2011.
 [35] PowerWorld. [Online]. Available: http://www.powerworld.com
 [36] P. M. Anderson and A. A. Fouad, Power System Control and Stability. John Wiley & Sons, 2008.
 [37] R. Zimmerman, C. MurilloSánchez, and R. Thomas, “MATPOWER: Steadystate operations, planning, and analysis tools for power systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 12–19, Feb 2011.
 [38] H. Ren, I. Dobson, and B. Carreras, “Longterm effect of the n1 criterion on cascading line outages in an evolving power transmission grid,” IEEE Trans. Power Syst., vol. 23, no. 3, pp. 1217–1225, 2008.
 [39] 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.
 [40] B. A. Carreras, D. E. Newman, I. Dobson, and N. S. Degala, “Validating OPA with WECC data,” in 2013 46th Hawaii International Conference on System Sciences (HICSS), Jan 2013, pp. 2197–2204.
 [41] [Online]. Available: web.engr.oregonstate.edu/â¼cotillaj/ecs/Cosmic.html
 [42] A. R. Bergen, Power Systems Analysis, 2nd ed. Pearson Education India, 2009.
Jiajia Song (S‘12) is currently pursuing his doctoral degree in Electrical Engineering at Oregon State University. His research interests focus on dynamic power system and protection modeling, cascading outages analysis, and phasor measurement unit applications. 
Eduardo CotillaSanchez (S‘08, M‘12) received the M.S. and Ph.D. degrees in Electrical Engineering from the University of Vermont in 2009 and 2012, respectively. He is currently an Assistant Professor in the School of Electrical Engineering & Computer Science at Oregon State University. His primary field of research is the vulnerability of electrical infrastructure, in particular, the study of cascading outages. CotillaSanchez is the Secretary of the IEEE Working Group on Understanding, Prediction, Mitigation and Restoration of Cascading Failures. 
Goodarz Ghanavati (S‘11) received the B.S. and M.S. degrees in electrical engineering from Amirkabir University of Technology, Tehran, Iran, in 2005 and 2008, respectively. Currently, he is pursuing the Ph.D. degree in Electrical Engineering at the University of Vermont, Burlington, VT, USA. He worked as a Design Engineer at Monenco Iran Co. His research interests include power system dynamics, stochastic modeling of power systems, phasor measurement unit applications and smart grid. 
Paul Hines (S‘96, M‘07) received the Ph.D. in Engineering and Public Policy from Carnegie Mellon University in 2007 and M.S. (2001) and B.S. (1997) degrees in Electrical Engineering from the University of Washington and Seattle Pacific University, respectively. He is currently an Associate Professor in the School of Engineering at the University of Vermont, and a member of the adjunct research faculty at the Carnegie Mellon Electricity Industry center. Currently Dr. Hines serves as the vicechair of the IEEE Working Group on Understanding, Prediction, Mitigation and Restoration of Cascading Failures. 