Finding $K$ Contingency List in Power Networks using a New Model of Dependency

TL;DR

This paper introduces the MIIR model to better understand power network dependencies and proposes a method to identify the most critical $K$ components whose failure could cause widespread outages, validated through simulations.

Contribution

The paper presents the MIIR model for power network dependencies and formulates the $K$ Contingency List problem, offering both optimal and heuristic solutions for complex failure analysis.

Findings

The heuristic provides near-optimal solutions significantly faster than the MIP.

MIIR effectively models complex power network dependencies.

Validation with 2011 Southwest blackout demonstrates model relevance.

Abstract

Smart grid systems are composed of power and communication network components. The components in either network exhibit complex dependencies on components in its own as well as the other network to drive their functionality. Existing, models fail to capture these complex dependencies. In this paper, we restrict to the dependencies in the power network and propose the Multi-scale Implicative Interdependency Relation (MIIR) model that address the existing limitations. A formal description of the model along with its working dynamics and a brief validation with respect to the 2011 Southwest blackout are provided. Utilizing the MIIR model, the $K$ Contingency List problem is proposed. For a given time instant, the problem solves for a set of $K$ entities in a power network which when failed at that time instant would cause the maximum number of entities to fail eventually. Owing to the problem being NP-complete we devised a Mixed Integer Program (MIP) to obtain the optimal solution and a polynomial time sub-optimal heuristic. The efficacy of the heuristic with respect to the MIP is compared by using different bus system data. In general, the heuristic is shown to provide near optimal solution at a much faster time than the MIP.

Tables

Table 1. TABLE I : IDRs of the Southwest Power System

Dependency Equations
$SE\leftarrow T1\cdot WECC+T4\cdot DE$
$DE\leftarrow T3\cdot PV+T2\cdot WECC$
$SONGS\leftarrow T5\cdot SE$
$IID\leftarrow T8\cdot DE+T15\cdot WALC+T9\cdot IV$
$MI\leftarrow T12\cdot IV$
$IV\leftarrow T13\cdot NG$
$CFE\leftarrow T14\cdot IV+T13\cdot MI$
$WALC\leftarrow T16\cdot NG$
$NG\leftarrow T11\cdot HA$
$HA\leftarrow T10\cdot PV$

Table 2. TABLE II : Quality of solution comparison of Mixed Integer Program and Heuristic for different Power Network bus systems and varying K 𝐾 K

		$K$ values and number of entities dead in MIP and heuristic
DataSet	Number of	1		2		3		4		5
	Entities	MIP	Heu	MIP	Heu	MIP	Heu	MIP	Heu	MIP	Heu
9 bus	$24$	$15$	$15$	$19$	$19$	$24$	$24$	$24$	$24$	$24$	$24$
14 bus	$44$	$29$	$29$	$32$	$31$	$36$	$35$	$42$	$41$	$44$	$43$
24 bus	$80$	$48$	$48$	$51$	$51$	$54$	$54$	$58$	$58$	$61$	$61$
30 bus	$83$	$54$	$54$	$58$	$58$	$64$	$63$	$71$	$70$	$74$	$73$
39 bus	$105$	$63$	$63$	$70$	$70$	$74$	$74$	$79$	$78$	$84$	$82$
57 bus	$149$	$96$	$96$	$113$	$108$	$128$	$119$	$136$	$128$	$140$	$136$
118 bus	$405$	$240$	$240$	$245$	$245$	$250$	$250$	$255$	$254$	$260$	$260$
145 bus	$666$	$499$	$499$	$511$	$511$	$517$	$516$	$522$	$517$	$527$	$523$
300 bus	$847$	$493$	$491$	$506$	$506$	$513$	$513$	$524$	$524$	$531$	$530$
2383wp bus	$5923$	$3249$	$3249$	$3296$	$3278$	$3326$	$3308$	$3354$	$3332$	$3382$	$3353$

Table 3. TABLE III : Run time comparison of Mixed Integer Program and Heuristic for different Power Network bus systems and varying K 𝐾 K

		$K$ values and running time in sec for MIP and in ms for Heuristic
DataSet	IDR Generation	1		2		3		4		5
	Time (ms)	MIP	Heu	MIP	Heu	MIP	Heu	MIP	Heu	MIP	Heu
9 bus	$0.31$	$0.026$	$2$	$0.13$	$1$	$0.02$	$1$	$0.02$	$1$	$0.01$	$1$
14 bus	$0.34$	$0.09$	$2$	$0.14$	$1$	$0.08$	$1$	$0.05$	$1$	$0.04$	$1$
24 bus	$0.60$	$0.27$	$4$	$0.16$	$4$	$0.19$	$6$	$0.18$	$3$	$0.27$	$3$
30 bus	$0.55$	$0.24$	$5$	$0.20$	$5$	$0.17$	$6$	$0.16$	$3$	$0.18$	$4$
39 bus	$0.62$	$0.30$	$8$	$0.17$	$14$	$0.28$	$6$	$0.28$	$7$	$0.25$	$5$
57 bus	$0.76$	$1.44$	$17$	$0.65$	$13$	$0.60$	$15$	$0.74$	$15$	$0.60$	$17$
118 bus	$1.49$	$1.78$	$33$	$1.09$	$52$	$1.16$	$64$	$0.73$	$33$	$1.02$	$41$
145 bus	$2.22$	$10.48$	$60$	$9.45$	$83$	$9.49$	$91$	$10.26$	$147$	$12.72$	$76$
300 bus	$2.09$	$11.09$	$146$	$4.23$	$249$	$4.68$	$364$	$3.63$	$253$	$4.34$	$293$
2383wp bus	$7.66$	$315$	$5$	$220$	$7$	$225$	$11$	$218$	$14$	$262$	$18$

Equations

$$\sum _{e_m\in O_{e_g}}{e}_{m,c,\tau }=\sum _{e_n\in I_{e_g}}{e}_{n,c,\tau }+e_{g,c,\tau },\forall e_g\in G$$

$$\sum _{e_m\in O_{e_l}}{e}_{m,c,\tau }=\sum _{e_n\in I_{e_l}}{e}_{n,c,\tau }-e_{l,c,\tau },\forall e_l\in L$$

$$\sum _{e_m\in O_{e_k}}{e}_{m,c,\tau }=\sum _{e_n\in I_{e_k}}{e}_{n,c,\tau }+e_{k,c,\tau },\forall e_k\in N$$

$$max\stackrel{|E|-1}{\sum _{i=1}}{x}_{i|E^2|}$$

Full text

I Introduction

The smart grid system constitutes both Power and Communication entities to execute different tasks and can be modeled as an interdependent Power-Communication network. Consider entities in either network. The electricity generation and power flows are partially controlled by the Supervisory Control and Data Acquisition System (SCADA) through signals received from Remote Terminal Units (RTUs). Meanwhile, every entity involved in sensing, sending and controlling the power grid (i.e. the Phasor Measurement Units (PMUs)) are dependent on power network entities to drive their functions. Owing to these dependencies, failure of some entities in either network may eventually result in a cascading failure causing widespread power blackouts.

In this article, we restrict to the analysis of cascading failure in smart grids considering power network dependencies in isolation. Consider that the control centers are not prone to cyber-attacks/physical failures and have additional power backup (independent of the underlying power grid) to make them operational under all circumstances. Even with such an assumption, studying the power network dependencies would provide potential information regarding the vulnerability and reliability of the smart grid system. It is to be noted that analysis under this assumption does bring into effect the failure of PMUs. PMUs are directly dependent for power on the associated power network entity. Hence failure of a power network entity would make the associated PMU non-operational. A brief description of the underlying analysis is provided. The power network entities can be broadly categorized into – Generation Bus, Load Bus, Neutral Bus and Transmission Lines. Consider a scenario in which there is an initial trip of a transmission line. This would cause some power to be re-routed through other lines to satisfy the load demands. Occasionally, this may result in some transmission lines to have power flow beyond their rated capacity limits causing them to trip. A series/cascade of such line trips might eventually cause a blackout. Regional power blackouts caused by cascading failure of components in power network have been seen previously in New York (2003)  [1], San Diego (2011)  [2] and India (2012)  [3]. From a system operator’s point of view, understanding the vulnerable components in the power network is critical to making well-judged decisions for prevention of power failure. An abstract model that captures these dependencies and an algorithm that uses this model to identify the critical components fast is beneficial from the operator’s point of view.

For the last few years, there has been a considerable surge of research being done in studying the complex dependencies in different critical infrastructures that cause cascading failures. Existing literature relies mostly on graph based models to represent the dependencies [4], [5], [6], [7], [8], [9], [10]. These models, however, fail to capture the complex dependencies that might exist which limits their applicability to real world systems. Authors in [11] bring out the need to address the complex dependencies in an inter-dependent or intra-dependent critical infrastructure(s) which can be explained through the following example. Let $e_{x},e_{w},e_{y},e_{z}$ (which can be a generator, substation, transmission line etc.) be power network entities. Consider the dependency where the entity $e_{x}$ is operational if (i) entities $e_{w}$ and (logical AND) $e_{y}$ are operational, or (logical OR) (ii) entity $e_{z}$ is operational. Models in [4], [5], [6], [7], [8], [9], [10], [12] fails to capture this kind of dependency. Motivated by these findings and limitations of the existing models, the authors in [11] proposed a Boolean logic based dependency model termed as Implicative Interdependency Model (IIM). For the example stated above, the dependency of $e_{x}$ on $e_{w},e_{y},e_{z}$ can be represented as $e_{x}\leftarrow e_{w}\cdot e_{y}+e_{z}$. This equation representing the dependency of an entity is termed as Interdependency Relation (IDR).

The biggest challenge of the IIM formulation is the identification of the dependency relations between the different entities. In the past, this was done based on information obtained from subject matter experts. However, it has now become clear that this is not a very reliable procedure. The second problem is that IIM operates on Boolean logic, implying that the different entities can only have two values — [math] or $1$, representing the state of the entity being operational or non-operational. However, this does not provide information about entities which are operating at near-failure state. For example, consider a generator that is operational but operating at $95\%$ of its peak generation capacity. The IIM model only provides the state of the generator (operational or non-operational). Hence the operator won’t be alerted even though the generator is reaching its peak generation capacity (which might eventually fail due to over generation). Thus, such scenarios limits the applicability of the IIM model.

In order to overcome the limitations of IIM and to extend the application domain of IDRs for long-term planning and short-term operational management, the Multi-scale Implicative Interdependency Relation (MIIR) model is proposed. The MIIR model uses the notion of IDRs with added features to capture the power flow in transmission lines and demand/generation of buses. PMU data is used to generate the dependency equations as well as obtaining power flow and demand/generation values.

Using the MIIR model we study the $K$ contingency list problem in this paper. At a given time $t$ the problem solves for a set of $K$ components in the power network which when made non-operational at time $t$ would cause the maximum number of entities to fail. Additionally, the solution would provide insights into components which are operating at their near capacity limits. Such a solution would have an immediate benefit to a system operator making decisions at that time to prevent large-scale power failures.

The paper is segregated as follows. In Section III, we describe the MIIR model in details and its applicability through a case study of real world power failure. The $K$ contingency list problem and its computational complexity are detailed out in Section IV. Owing to the problem being NP-complete, a mixed integer program is provided that solves the problem optimally in Section V. As it is critical to solving the problem fast, a sub-optimal heuristic with polynomial time complexity is discussed in the same section. A detailed experimental analysis for different bus systems is provided in Section VI.

II Related Work

The Presidential Commission on Critical Infrastructures, appointed by President Clinton in 1996 [13], was one of the earliest studies done on understanding robustness and resiliency of Critical Infrastructure Interdependency. Rinaldi et al.[14], [15] are among the first group of researchers to study the inter-dependency between critical infrastructures and proposed a complex adaptive system to model the inter-dependencies. A survey of Critical Infrastructure Interdependency modeling is presented in Pederson et al. in [16] which was undertaken by the U.S and a group of international researchers. As discussed before, there have been significant research to capture these inter-dependencies using graph-based models [4], [5], [6], [7], [8], [9], [10]. In general, the abstractions provided by the graph-based models fail to capture the complex inter-dependencies that might exist.

An event driven co-simulation framework for interdependent power-communication network was presented by Lin et al. in [17], [18]. The authors in [12] presented a game theoretic model for a multi-layer infrastructure networks using flow equilibrium. Authors in [19] studied the security of interdependent and identical Networked Control System (NCS). They used a discrete-time stochastic linear system to model each plant with a shared communication network controlling the systems. In [20] the importance of considering power and communication infrastructures simultaneously was highlighted. Based on the Polish power grid blackout, a study of the impact of overhead lines tripping order on the severity of failure is presented in [21]. Analyzing failure in the smart grid under targeted initial attack was studied in [22]. The effect of cyber (communication) and power network dependencies in the smart grid was studied in [23] for reliability assessments. Recovery of information of the failed entities in a power grid after a failure event was studied in [24]. Topology based vulnerability assessment in power netowkr has been reported in [25], [26], [27]. As we rely on PMU data of the power network for generating the dependency equations, we investigated the current status of utilizing PMU data in addressing the complex inter-dependency issues. The roles of real-time measurements obtained from IEDs such as PMUs in mitigating failures have been investigated in [28], [29], [30], [31]. However, most of the analysis in these papers does not gives insight about the instantaneous identification of most vulnerable entities in a power network.

III The MIIR Model

In this paper, the MIIR model is developed based on the power network. It is aimed to have a near accurate abstraction of the power flow dynamics and capture cascading failure propagation in the same. It is to be noted with proper modification the model has the potential to be extended for performing a similar analysis in different inter/intra-dependent critical infrastructure system(s).

III-A Model Variables

We consider load buses, generator buses, neutral buses (or zero injection buses) and transmission lines as different types of entities. Let $E=\{e_{1},e_{2},\ldots,e_{n}\}$ denote the set of entities in the power network. Each entity $e_{i}\in E$ has three values associated with it — (i) a lower bound $e_{i,l}$, (ii) an upper bound $e_{i,u}$, and (iii) the instantaneous power value $e_{i,c,t}$ at time $t$ of the entity. For a transmission line type entity $e_{k}$, the value of $e_{k,c,t}$ provides the power flow in that line at time $t$. Corresponding, for a load bus $e_{m}$ and generator bus $e_{n}$ the values $e_{m,c,t}$ and $e_{n,c,t}$ provides the load demand and generating value at time $t$. For the power transmission and distribution systems, PMU data can be used to get the instantaneous power value $e_{i,c,t}$ of the entity $e_{i}$. The values of $e_{i,l}$ and $e_{i,u}$ can be easily obtained from the entity rating data. For a given time $t$ the state of the entity is still Boolean (operational or not operational) and is guided by the following two factors — (a) $e_{i,c,t}$ satisfies the property $e_{i,l}<e_{i,c,t}<e_{i,u}$, (b) the corresponding dependency equation of $e_{i}$ at time $t$ is satisfied. Hence, if $e_{i}$ has an IDR $e_{i}\leftarrow e_{j}\cdot e_{k}+e_{l}$ then for $e_{i}$ to be operational both the properties — (a) ($e_{i}$ and $e_{k}$) or $e_{l}$ is operational at $t-\delta$, and (b) $e_{i,t}\leq e_{i,c,t}\leq e_{i,u}$ has to be satisfied. Here $\delta$ refers to the time within which the effect of failure of an entity is propagated to its dependent entity. So using MIIR model the power network at time $t$ is mathematically represented as $P(E,B,C_{t},F)$ where $E$ is the set of of entities, $B$ is a set of tuples $\{e_{i,l},e_{i,u}\}$ ($\forall e_{i}\in E$) denoting the power value bound on the entity, $C_{t}$ consist of instantaneous power value $e_{i,c,t}$ ($\forall e_{i}\in E$) at time $t$ and $F$ contains the set of dependency equations for the entities in $E$. A similar notation has been used in [26] but our notation brings out a completely different topological aspect of the power network.

III-B Generating and Obtaining the Sets $B,C_{t}$ and $F$ for a Power Network

We illustrate our strategy to generate the dependency equations $F$ and the set $C_{t}$ of a power network $P(E,B,C_{t},F)$ at a given time $t$. The PMU data is substituted with simulated data due to its unavailability for different bus systems. The MATPOWER [32] software is used to generate the simulated data. For a given time $t$ and a standard bus system (containing a set of buses and transmission lines $E$), the software use load demand of the bus, the impedance of the transmission lines etc. to solve the power flow. The software produces the voltage of each bus in the system as the output. The software suite also includes a wide range of bus systems along with power ratings of the components for all such systems. We restrict ourselves to analyze the real power flow. Firstly, for a given solution, we formally state the procedure to obtain the tuple values of the set $B$ and instantaneous power value contained in the set $C$ (all values in MW) for generator buses, load buses, neutral buses and transmission lines —

• •

*Generator Bus: *The real part of the power generated is taken as the value of $e_{a,c,t}$ for a generator bus $e_{a}\in E$. The upper bound $e_{a,u}$ is set to its real generation capacity (supplied in the MATPOWER suite) and the lower bound $e_{a,l}$ is set to [math]. It is to be noted that some generator buses have load demand. Consider $e_{x}$ be a generator bus with load demand $d$ units and real instantaneous power generated $e_{x,c,t}$ units. Without the loss of generality, such a bus is split into a generator bus $e_{x1}$ with [math] load demand (instantaneous power generated $e_{x1,c,t}$) and a load bus $e_{x2}$ with instantaneous load demand $d$ units (instantaneous power generated [math]). A transmission line $e_{x12}$ is constructed that connects $e_{x1}$ to $e_{x2}$ with an instantaneous power flow of $d$ units flowing from $e_{x1}$ to $e_{x2}$.

• •

*Load Bus: *The real part of the load demand is taken as instantaneous demand value $e_{b,c,t}$ of a load bus $e_{b}\in E$. For a load bus $e_{b}$, both its upper and lower bound is set to the instantaneous demand value $e_{b,c,t}$. Essentially, our assumption is that a load bus does not change its demand value irrespective of any failure.

• •

*Neutral Bus: *For a netrual bus $e_{d}\in E$ the values of $e_{d,l},e_{d,u}$ and $e_{d,c,t}$ are set to [math].

• •

*Transmission Lines: *For two buses $e_{1}$ and $e_{2}$ connected by a transmission line $e_{12}$ the power flowing through the transmission line is calculated as $P_{12}=Real(V_{1}*(\frac{V_{1}-V-2}{I_{12}})^{*})$, where $V_{1}$ is the voltage at bus $e_{1}$, $V_{2}$ is the voltage at bus $e_{2}$ ($V_{1}$ and $V_{2}$ returned by the MATPOWER solver) and $I_{12}$ is the impedance of the transmission line $e_{12}$ (obtained from the supplied bus system file of MatPower). $P_{12}$ is the real component of the power flowing in the transmission line $e_{12}$. The lower bound is set to [math] and the upper bound is taken as the rated capacity of the transmission line. The instantaneous power value $e_{12,c,t}$ is set to $|P_{12}|$ (absolute value).

The description of the transmission line type entity provided above shows that power flows from bus $e_{1}$ to $e_{2}$ if $P_{12}$ is positive and vice versa otherwise. As a result, we can interpret the direction of power flow in the line from the solution which we obtained from MatPower. We use this solution to generate the set F which is the set of dependency equations. As an example, we consider the nine bus system which is shown in Figure 1. This figure describes a power network $P(E,B,C_{t},F)$, at time instance $t$ with $E$ being the set of entities containing generator buses from $G_{1}$ to $G_{3}$, load buses $L_{1}$ to $L_{4}$, neutral buses $\{N_{1},N_{2}\}$ and transmission lines $T_{1}$ through $T_{9}$. The figure also provides the instantaneous power values by solving the power network flow based on demand/generation at some time instant $t$. The red blocks denote the instantaneous real power generated by a generator, the green blocks denote instantaneous real load demands and the blue nodes are neutral. The values in the grey blocks denote the flow of power in the transmission lines with the arrows denoting the direction of power flow. There aren’t any IDRs for the transmission lines. The IDRs for bus $b_{1}$ is created as follows: (a) let $b_{2}$, $b_{3}$ be the buses and $b_{12}$ (between $b_{1}$ and $b_{2}$) and $b_{13}$ between ($b_{1}$ and $b_{3}$) be the transmission lines for which power flows from these buses to $b_{1}$, (b) the dependency equation for the bus $b_{1}$ is constructed as disjunction of minterms of size $2$ (this consists of the bus from which the power is flowing and the respective transmission line) with each disjunction corresponding to buses from which power is flowing to it. For this example the dependency equation $b_{1}\leftarrow b_{12}b_{2}+b_{13}b_{3}$ is created. Using this definition, the dependency equations for the buses in Figure 1 are as follows — (a) $L_{1}\leftarrow T_{1}\cdot G_{1}$ , (b) $L_{2}\leftarrow T_{2}\cdot L_{1}+T_{7}\cdot N_{2}$, (c) $L_{3}\leftarrow T_{3}\cdot L_{1}+T_{4}\cdot N_{1}$, (d) $L_{4}\leftarrow T_{6}\cdot N_{1}+T_{8}\cdot N_{2}$, (e) $N_{1}\leftarrow T_{5}\cdot G_{3}$, (f) $N_{2}\leftarrow T_{9}\cdot G_{2}$.

III-C Dynamics of the MIIR model

To understand the dynamics of cascading failure in power network based on the MIIR model we first create the abstract representation $P(E,B,C_{0},F)$ (which is constructed using technique discussed in Section III-B) for a power network at time $t=0$. An event of initial failure is assumed to occur at time $t=0$ with failure cascade propagating in unit time steps. For an entity that is operational at time step $t=\tau$ the following equations are required to be satisfied —

[TABLE]

[TABLE]

[TABLE]

Equations 1-3 dictates the law of conservation of energy for each bus in the system. We assume that there is no delay in time between the power flowing out from a bus and power flowing into it.

In equation 1, the transmission line entities that flows power out and into the generator bus $e_{g}$ (where $G\subset E$ contains all generator buses) are represented by sets $I_{e_{g}}$ and $O_{e_{g}}$ respectively. Equation 2 and 3 uses the same notations of transmission line for load bus $e_{l}$ (where $L\subset E$ contains all load buses) and neutral bus $e_{k}$ (where $N\subset E$ contains all neutral buses) respectively. In Equation 3, the value of $e_{k,c,t}=0$ for all time steps and hence it can be simplified as $\sum_{e_{m}\in I_{e_{k}}}e_{m,c,t}=\sum_{e_{n}\in O_{e_{k}}}e_{n,c,t}$. Additionally, for a generator bus there is no power injected to it through transmission lines. Hence equation 1 can be re-written as $\sum_{e_{m}\in I_{e_{g}}}e_{m,c,t}=e_{g,c,t},\forall e_{g}\in G$. We use this abstract representation of the power network i.e. $P(E,B,C_{0},F)$ at time $t=0$. Using this the cascading failure process of the power network on an event of initial failure of $E^{\prime}\subset E$ at time step $t=0$ is detailed out in Algorithm 1.

In Algorithm 1, at every iteration the flow values are adjusted based on Equations 1-3 at line $8$ and entities are made non-operational based on the two conditions mentioned at lines $7$ and $9$. The cascading process continues if a new entities fail in the previous time step (condition $size\neq|S|$). As evident we assume that there is a delay of $1$ time unit for an entity to become non operational if its dependency equations are unsatisfied. All entities whose bound values are not satisfied are made non operational at that time step after power flow calculation. It is to be noted that the dependency equations are generated from a graph which is directed acyclic. Owing to this property the cascade reaches a steady (no new entities are non-operational) within $\mathcal{O}(|E|)$ times steps. This can be explained as follows. Consider a single initial failure of an entity. If no entity fails the cascading algorithm would continue till at most $|E|-1$ time steps since the maximum distance between two nodes in a directed acyclic graph is $|E|-1$ (considering each edge having a weight $1$). If more than one entity fails then the cascade is expected to stop before $|E|-1$. Hence the number of cascading time steps is strictly upper bounded by $|E|-1$.

Algorithm 1 assumes that there exist a method to compute the flow value equations to get the instantaneous power values of entities. This is equivalent to computing the AC power flow equations again (using MatPower) which would be time intensive and does not make use of the abstraction created by the MIIR model to make any decision. Moreover, using general graph theoretical algorithmic techniques might result in multiple solutions of instantaneous power values when solving a given set of power flow equations thus resulting in ambiguity. To counteract this, in our abstraction, we use the notion of Worst-Case Cascade Propagation (WCCP) in Algorithm 1. Qualitatively, the instantaneous value of power flows and power generator at every time step $t>0$ of the cascade is set to a value that would cause the maximum number of entities to fail at the end of the cascade. Computation of this power flow values using WCCP is proved to be NP-complete in Section IV. We devise a mixed integer porgram to get the optimal solution and a greedy heuristic to get a sub-optimal solution in polynomial time in Section V that addresses WCCP.

III-D Case Study: The 2011 Southwest Blackout

In this subsection the performance of MIIR with WCCP is tested on a real power system event: the 2011 Southwest Blackout. All data used in this analysis are obtained from [2]. An abstraction of the Southwest Power System is provided in Figure 2. The abbreviations used in Figure 2 are — Western Electricity Coordinating Council (WECC), Serrano (SE), Devers (DE), San Onofre Nuclear Generating Station (SONGS), San Diego Gas & Electric (SDG & E), Miguel (MI), Imperial Valley (IV), Imperial Irrigation District (IID), Comisión Federal de Electricidad’s (CFE, corresponding to Baja California Control Area), North Gila (NG), Hassayampa (HA), Palo Verde (PV), and Western Area Power Administration-Lower Colorado (WAPA). The blue, orange and green blocks in Figure 2 represents neutral, load and generator buses respectively. The transmission lines are labeled $T1-T17$ with the arrows indicating the directions of the pre-disturbance power flows. On September 8, 2011, an initial trip of the HA-NG transmission line ($T11$) caused blackout in SDG & E region. The objective here is to see whether MIIR model with worst case cascade propagation is able to capture the power outage.

The dependency equations in Table I without the bounds and instantaneous power values of the entities (buses and transmission lines) corresponds to the set $F$. Consider tripping of the entity $T11$ at $t=0$. Just considering the IDRs itself, the component $NG$ fails at $t=1$, $WALC$ and $IV$ at $t=2$, $CFE$ and $MI$ at $t=3$. The pre-disturbance load demands of SDG&E and IID were approximately 5000 MW and 900 MW, respectively, while the generation bounds on PV and WECC were [0, 4000 MW] and [0, 10000 MW], respectively. After failure of T11, SDG&E and IID would try to meet their bulk load demands through the generator buses PV and WECC via T6 and T8. The bound on T6 is [0, 2200 MW] and T8 is [0, 1800 MW]. Both PV and WECC have enough generation capacity to meet the load demand of SDG&E and IID. At $t=3$ owing to the load demand of SDG & E the transmission line $T6$ would have try to have a power flow of $5000$ MW instantly. Thus $T6$ would trip at $t=3$ causing SDG&E to trip at $t=4$. Owing to this the power flowing through $T1,T4$ and $T5$ would reduce down to [math] at $t=4$. The power flow in $T8$ would increase to $900$ MW at $t=3$ for supplying power to IID. Thus the steady state is reached at $t=4$ and MIIR model accurately predicts the blackout of SDG & E region.

IV Problem Formulation

It is important from a power system operator’s point of view to understand and know the most critical entities in the network at a given time. This would enable the operator to make more reliable decisions on an event of some failure. For larger systems, an automation that provides the operator with this information would be hugely beneficial. Owing to this we develop the $K$ Contingency List (KCoL) problem using MIIR model with WCCP. For a given time $t$ and an integer $K$ the problem provides the operator with a list of $K$ entities which when failed initially causes the maximum number of entities to fail at the steady state of cascade propagation. Qualitatively, for a given integer $K$ the problem finds a set $E^{\prime}$ ($|E^{\prime}|=K$) entities which when failed initially maximizes the total number of entities failed at the end of the cascading process. A formal description of the KCoL problem using WCCP of MIIR model for Power Network is provided —

Input: (a) A power network $P(E,B,C_{t},F)$ where $E=G\cup L\cup N\cup T$. Set of entities $G$, $L$, $N$ and $T$ are disjoint and contains the generator buses, load buses, neutral buses and transmission lines respectively. (b) two positive integers $K$ and $S$.

**Decision Version: ** Does there exist a set of $K$ entities in $E$ whose failure at time $t$ would result in a failure of at least $S$ entities in total at the end of the cascading process.

**Optimization Version: ** Compute the a set $K$ entities in a power network $P(E,B,C_{t},F)$ whose failure at time $t$ would maximize the number of entities failed at the steady state of cascade propagation.

We prove the problem is NP-complete to solve in Theorem 1

Theorem 1**.**

The KVNE problem using MIIRA model is NP-complete.

Proof.

The problem is proved to NP-complete by a reduction from the densest $p-$subhypergraph problem [33]. An instance of the densest $p-$subhypergraph problem consists of a hypergraph $H=(V,E)$ and two parameters $p$ and $M$. The decision version of the problem finds the answer to whether there exists a set of vertices $V^{\prime}\subset V$ and $|V^{\prime}|\leq p$ which completely covers at least $M$ hyper-edges.

From an instance of the densest $p-$subhypergraph we create an instance of the KCoL problem as follows. We start with an empty set of entities $G$, $L$ and $T$ and an empty set $F$ that would comprise of the dependency equations. A load type entity $L_{j}$ is added to set $L$ for each hyper-edge $E_{j}\in E$ with instantaneous load demand $L_{i,c,t}$ set to the number of vertices that comprise this hyper-edge. For each vertex $V_{i}\in V$ we add a generator type entity $G_{i}$ to set $G$. The upper bound on the capacity of the generator $G_{i}$ is set to the sum of all instantaneous load demands $L_{i,c,t}+1$ for which the corresponding hyper-edge $E_{j}$ contains the vertex $V_{i}$. For each hyper-edge $E_{j}$ consisting of vertices $V_{x},V_{y},V_{z}$ (say) three transmission line type entities $T_{x},T_{y}$ and $T_{z}$ are added to set $T$ and a dependency equation $L_{j}\leftarrow T_{x}\cdot G_{x}+T_{y}\cdot G_{y}+T_{z}\cdot G_{z}$ is created and added to set $F$. The upper bound of the transmission line is set to the load demand $+1$ of the entity it connects to (e.g., in this case, the maximum capacity of each transmission line $T_{x},T_{y},T_{z}$ are set to the instantaneous load demand $L_{i,c,t}+1$). The parameter $S$ of KCoL problem is set to $p+M$ and $K$ is set to $p$ (i.e. $p$ entities fail at time $t$). Thus the created instance satisfy the property of the graph from which the dependency relations are computed being Directed Acyclic. In the initial operating condition at time $t$, all transmission lines have a line flow value of $1$ unit with each generator $G_{i}$ producing $P_{i}$ units of power, where $P_{i}$ is the number of load entities it is connected to. Hence all load demands are satisfied. It can be directly followed that an instance of KCoL problem can be created from an instance of densest $p-$subhypergraph problem in polynomial time.

It is to be noted that for the created instance — (1) Each transmission line has the capacity to satisfy the complete load demand of the load type entity it is connecting, (ii) Each generator has the capacity to satisfy the load demand of all the load type entities it is connected to. Hence an initial failure of one or more entities would not cause any transmission line or generator to trip (fail) because of exceeding its maximum capacity. Thus the generators and transmission lines are susceptible only to initial failure whereas the load entities are vulnerable to both initial and induced failures. However, failure of load entities can not cause any induced failure. Induced failure of the load entity can be caused only when each minterm in its dependency equation have at least one failed entity. Thus no entity fail due to change in power flow values.

Now consider there exist a solution to the densest $p-$subhypergraph problem. Hence there exist a set of $p$ vertices $V^{\prime}$ that completely covers $M$ hyper-edges. Failing the generator type entities corresponding to the vertices in $V^{\prime}$ would thus fail at least $M$ load entities at $t+1$ according to the instance construction. Thus a total of at least $p+M$ entities would fail which solves the KCoL problem. On the other way round consider there exist a solution to the KCoL problem. As reasoned earlier, a load entity cannot cause any induced failure. Hence if a load entity is in the solution then it can be substituted with any operational generator entity without loss of correctness. Similarly, if a transmission line type entity is in the solution it can be replaced by a generator type entity it is connected to. Using this substitution a solution thus comprises of entities $G^{\prime}\subset G$. All $M$ (or greater than $M$) entities that fail due to the initial failure of $p$ entities belongs to set $L$. Thus the substituted solution (or original solution if no substitution is required) would consist of generator type entities that cause failure of these $M$ (or greater than $M$) load entities. Hence selecting the vertices corresponding to $G^{\prime}$ would ensure that at least $M$ hyper-edges are completely covered solving the densest $p-$subhypergraph problem. Hence proved.

∎

V Solutions to the Problem

Owing to the KCoL problem being NP-complete, we obtain the optimal solution using an Mixed Integer Program (MIP). However, as we require to compute the contingency list fast, we devise a polynomial time heuristic that provides a sub-optimal solution to the problem.

V-A Optimal Solution using Mixed Integer Program

As a reference frame, we consider that the initial failure occurs at time step $t=0$. It is shown in Section III-C that the number of time steps in the cascade is upper bounded by $|E|-1$. We devise a Mixed Integer Program that solves the KCoL problem optimally for a power network $P(E,B,C_{0},F)$ (the abstraction constructed for $t=0$). Irrespective of whether the steady state is reached before or at time step $|E|-1$, in our MIP we try to maximize the number of entities failed at $t=|E|-1$ when $K$ entities fail at $t=0$. Moreover, it can not be predicted when the cascading failure stops. Hence, the MIP is bound to check for solution to compute the maximum number of entities that can fail till the maximum possible time step, i.e. $|E|-1$. Firstly, the list of variables used in the MIP formulation are discussed—

• •

*Variable List 1: * For each entity $e_{i}\in E$ a variable set $x_{i,t},\forall t,0\leq t\leq|E|-1$ are created. The value of $x_{i,t}$ is [math] if the entity is operational at time step $t$ and $1$ otherwise.

• •

*Variable List 2: * For each entity $e_{i}\in E$ a variable set $y_{i,t},\forall t,0\leq t\leq|E|-1$ is created. From the set $C_{0}$ we can get the initial instantaneous power value $e_{i,c,0}$ of an entity $e_{i}$. The value of $y_{i,0}$ is set to $e_{i,c,0}$. All the instantaneous values are real thus comprising the set of non integer variables in the program.

Using these definitions and the list of variables created, the objective of the MIP is provided in 4 and the constraints of the MIP are formally described.

[TABLE]

Subjected to:

Constraint Set 1: $\sum_{i=1}^{|E|}x_{i,0}=K$. This constraint sets the number of entities failed at time step $t=0$ to $K$.

Constraint Set 2: $x_{i,d}\geq x_{i,t-1},\forall t,1\leq t\leq|E|-1$. This ensures that an entity that is not operational at time step $t=d$ would remain non-operational in all times step $t>d$.

Constraint Set 3: Consider an IDR of form ${e_{i}}\leftarrow e_{a}\cdot e_{b}+e_{c}\cdot e_{d}$. To capture the cascading failure process, a set of constraints is developed and described below —

Step 1: New variables are introduced to represent the minterms. In this example, two new variables $c_{ab}$ and $c_{cd}$ are created to represent the terms $e_{a}\cdot e_{b}$ and $e_{c}\cdot e_{d}$. This is equivalent of adding two new IDRs $c_{ab}\leftarrow e_{a}\cdot e_{b}$ and $c_{cd}\leftarrow e_{c}\cdot e_{d}$ with the transformed IDR being ${e_{i}}\leftarrow c_{ab}+c_{cd}$.

Step 2: A linear constraint is developed for the $c$ type variables to capture the failure propagation. For an IDR $c_{ab}\leftarrow e_{a}\cdot e_{b}$, the constraint is represented as $c_{ab,t}\leq x_{a,t-1}+x_{b,t-1},\forall t,1\leq t\leq|E|-1$. This captures the condition that $c_{ab,t}$ is equal to $1$ only if at least one of the entities $e_{a}$ or $e_{b}$ is non operational.

Step 3: For each transformed IDR a linear constraint is introduced. For an IDR ${e_{i}}\leftarrow c_{ab}+c_{cd}$ the constraint is represented as $N\times x_{i,t}\leq c_{ab,t-1)}+c_{cd,t-1},\forall t,1\leq t\leq|E|-1$. Here $N$ is the number of minterms in the IDR (in this example $N=2$).

Constraint Set 4: For a given load bus entity $e_{l}$, the constraint $y_{l,t}=0,\forall t,0\leq t\leq|E|-1$ is added denoting that the instantaneous power demand of all the load bus remain constant at each time step. Similarly, for a given neutral bus entity $e_{n}$, the constraint $y_{n,t}=0,\forall t,0\leq t\leq|E|-1$ is added.

Constraint Set 5: For a given generator bus entity $e_{p}$ and transmission line entity $e_{q}$, the constraints $x_{p,t}\leq\frac{y_{p,t+1}}{e_{p,u}},\forall t,1\leq t\leq|E|-1$ and $x_{q,t}\leq\frac{y_{q,t+1}}{e_{q,u}},\forall t,1\leq t\leq|E|-1$ are added. As this is a maximization problem, the $x$ type variable of the corresponding generator/transmission line entity would be set to $1$ when it operates beyond its rated upper bound. The constraints $y_{a,t}\geq 0$ and $y_{a,t}\leq e_{a,u}+1$ are added at all time steps for each generator or transmission line type entity $e_{a}$. This limits the maximum value of these entities to its upper bound plus one and them failing only if their instantaneous power value is just above the upper bound.

Constraint Set 6: To capture the power flow equations in equations 1-3 the following constraints are developed. Consider the equation $\sum_{e_{m}\in O_{e_{l}}}e_{m,c,t}=\sum_{e_{n}\in I_{e_{l}}}e_{n,c,t}-e_{l,c,t}$. Naively, this can be constructed as a non-linear constraint $\sum_{e_{m}\in O_{e_{l}}}(1-x_{m,t})\times y_{m,t}=\sum_{e_{n}\in I_{e_{l}}}(1-x_{n,t})\times y_{n,t}-(1-x_{l,t})y_{l,t+1}$. The constraint denotes that the instantaneous flow values of the transmission lines/load demand are taken into consideration if the the load bus is operational at the next time step (as failure due to IDR is reflected after $1$ unit of time). This constraint can be linearized as $\sum_{e_{m}\in O_{e_{l}}}(y_{m,t}-x_{m,t}\times e_{m,u})=\sum_{e_{n}\in I_{e_{l}}}(y_{n,t}-x_{n,t}\times e_{n,u})-(y_{l,t}-x_{l,t+1}e_{l,u})$. If a transmission line $e_{n}$ fails at time instant $t$ then its instantaneous power value is set to its upper bound (owing to constraint set $5$). This would equate the term $(y_{n,t}-x_{n,t}\times e_{n,u})$ corresponding to this transmission line to [math]. If the transmission line $e_{n}$ is operational then $x_{n,t}=0$ and hence $(y_{n,t}-x_{n,t}\times e_{n,u})$ would equate to $y_{n,t-1}$ thus being considered in the power flow equation. Similarly if the load bus is not operational the value of $(y_{l,t}-x_{l,t+1}e_{l,u})$ is set to [math]. These constraints are constructed for all time steps $0\leq t\leq|E|-2$ and similar constraints are generated for equations 1 and 3 as well.

*Constraint Set 7:*For each transmission line type entity $e_{a}\in E$ flowing out power from a bus type entity $e_{b}$ the constraint $x_{a,t}\leq x_{b,t}$ is added for each time step $1\leq t\leq|E|-1$. This captures the condition that if a bus type entity fails then all transmission lines to which it transmits power also fails.

It is to be noted that there won’t be any infeasibility in solution arising due to the constraints. The load and neutral buses can only be made non-operational through their dependency equations. Whereas, the transmission lines and generators can be only made non-operational through change in power flow/generation values (as they don’t have any dependency equations). The objective in (4) along with these set of constraints, finds the the set of $K$ entities whose initial failure at $t=0$ maximizes the number of entities failed at the end of cascading process. As this is a maximization problem the power flow and generation at each time step is set to values that maximizes the total number of entities failing at the steady state. Hence the MIP captures the notion of WCCP.

V-B Heuristic Solution

In this section we design a sub-optimal heuristic that finds a solution to the KCoL problem in polynomial time. Primarily we use the following two definitions that drive the heuristic.

• Definition:

Kill Set: For an entity $e_{j}\in E$ in an power network $P(E,B,C_{t},F)$, the Kill Set is denoted as $KS(e_{j})$. This set comprises of all the entities made non-operational including the entity itself at the end of the cascading process when entity $e_{j}$ is made non-operational initially.

• Definition:

Fractional Minterm Hit Value: For an entity $e_{j}\in E$ in an power network $P(E,B,C_{t},F)$ the Fractional Minterm Hit Value is denoted as $FMHV(e_{j})$. It is calculated as $FMHV(e_{j})=\sum_{i=1}^{m}\frac{c_{i}}{|s_{i}|}$ where $m$ are the minterms and for a given minterm $m_{i}$, $c_{i}$ are the number of entities that belong to $Kill(e_{j})$ and $|s_{i}|$ is its size. This metric provides an estimate of impact of other operational entities that can be made non-operational at future time steps if the entity $e_{j}$ is made non-operational.

Algorithm 2 returns a sub-optimal value of $E^{\prime}$ which when failed initially would greedily maximize (based on Kill Set and FMHV) the number of entities failed at the end of the cascade. The algorithm runs in $\mathcal{O}(n^{3})$ where $n=|E|$. It is to be note that the greedy failure maximization is done based on IDR. To get the actual number of entities failed when the set of entities $E^{\prime}$ fail initially we use the MIP. Essentially, we modify the constraint 1 such that only entities in $E^{\prime}$ fail at $t=0$ and see the number of entites failed at the final time step. This gives us a measure to compare the efficacy of the heuristic solution with respect to the MIP.

VI Experimental Results

We analyzed the run time performance and quality of the heuristic solution with respect to MIP for different bus systems. The quality of the solution is defined by the number of components reported to be non-operational for a given value of $K$. Specifically we used the $9,14,24,30,39,57,118,145,300$ and $2383$ Winter Polish bus systems available in MATPOWER. For a given bus system, we used the MATPOWER AC power solver. Using the data, the abstract power network $P(E,C_{0},F)$ was generated. On the constructed power network the MIP and heuristic solutions were executed. The implementation was done in Java and a student licensed version of IBM CPLEX optimizer was used to solve the MIP. A UNIX system with 8 GB of RAM and intel i5 processor was used for the execution. All the implementation along with data sets are made open source. The repository can be found in https://github.com/jbanerje1989/ContingencyList.

In Table II, a comparison between the MIP and the heuristic solution with respect to the number of entities in non-operational state for different bus systems with $K$ varied from $1$ to $5$ in steps of $1$ are provided. Additionally, the total number of entities (buses and transmission lines) for each bus system is mentioned. It is to be noted that the total number of entities and the number of entities in non-operational state include the entities constructed for generator buses with non-zero load demand (as mentioned in Section III-B). Table III reports the IDR generation time for each bus system along with the time taken to execute the MIP and heuristic for different values of $K$. Some insightful observations from the results are as follows — (a) The heuristic solution performs very nearly to that of MIP with respect to quality and have an almost same performance for $K=1$. (b) For almost all the cases, the maximum percent difference in the number of non operational entities in heuristic with respect to MIP is under $1\%$ with a maximum percent difference of $7\%$ for $57$ bus system at $K=3$. (c) It is observed that more than $50\%$ of the total entities in a given bus system will be non-operational if $K=1$. This implies that the power system are extremely vulnerable even if a single entity is attacked, (d) For almost all the bus systems from $9$ to $300$ the heuristic finds a solution to the KCol problem nearly $100$ times faster than the MIP, (e) for the $2383$ Winter Polish bus system the heuristic is $10$ to $20$ times faster. However, it is be noted that the comparison is done based on a serialized implementation of the heuristic. There exists scope to parallelize the heuristic to achieve even faster run time.

Hence, it can be reasonably argued that the Heuristic solves the KCol problem achieving near optimal solution at a much faster time compared to MIP. Thus the abstraction provided by the MIIR model along with the Heuristic solution can be used by a Power Network operator to obtain the $K$ Contingency List at any given time and make appropriate control decisions.

VII Conclusion

In this paper, we introduce the MIIR model to capture dependencies and analysis of cascading failure in power network. We used the model to solve the $K$ Contingency list problem that finds a set of $K$ entities which when made non-operational at a certain time step would cause the maximum number of entities to fail at the end of cascading process. As the problem is NP-complete we devise a MIP to obtain the optimal solution and a polynomial time sub-optimal heuristic. The time and quality performance of the heuristic with respect to the MIP are compared for different bus systems and it is shown that the heuristic outperforms the MIP with respect to running time and additionally produces near optimal solution.