Structural Vulnerability of Power Grids to Disasters: Bounds, Adversarial Attacks and Reinforcement
aa r X i v : . [ c s . S Y ] S e p Structural Vulnerability of Power Grids to Disasters:Bounds, Adversarial Attacks and Reinforcement
Deepjyoti Deka,
Student Member, IEEE and Sriram Vishwanath,
Senior Member, IEEE
Abstract —Natural Disasters like hurricanes, floods or earth-quakes can damage power grid devices and create cascadingblackouts and islands. The nature of failure propagation andextent of damage is dependent on the structural features of thegrid, which is different from that of random networks. Thispaper analyzes the structural vulnerability of real power gridsto impending disasters and presents intuitive graphical metricsto quantify the extent of damage. Two improved graph eigen-value based bounds on the grid vulnerability are developedand demonstrated through simulations of failure propagation onIEEE test cases and real networks. Finally this paper studiesadversarial attacks aimed at weakening the grid’s structuralresilience and presents two approximate schemes to determinethe critical transmission lines that may be attacked to minimizegrid resilience. The framework can be also be used to designprotection schemes to secure the grid against such adversarialattacks. Simulations on power networks are used to compare theperformance of the attack schemes in reducing grid resilience.
I. I
NTRODUCTION
The topological structure of the power grid is an importantfeature that affects the delivery of electricity [1]. From aneconomic perspective, the capacity of transmission lines andthe graphical properties (whether tree-like or loopy) affect thelocational marginal electricity prices as well as the conver-gence of Optimal power flow algorithms [2]. Insufficient linecapacities and network structure can lead to highly fluctua-tion prices and even negative prices [3]. From an reliabilityperspective, the grid structure influences the extent of damagefollowing an natural or man-made disaster. In particular, itaffects the propagation of failures after an initial breakdownof equipment and thus in turn affects the formation of islandsand loss of service. Over the years, natural disasters likeearthquakes, floods, hurricanes have caused extensive poweroutages due to damage of grid equipment and loss of networkconnectivity [4], [5]. Such disasters also affect co-located in-terdependent transportation and communication infrastructureas well. Thus there is a greater need to quantify the effect ofthe grid structure on failure propagation in the grid followinga natural disaster and to incorporate the insights gained intotransmission planning techniques to improve grid resilience.There is existing work that studies the impact of structure ongrid reliability. Reference [6] describes a widely accepted andrealistic DC model for the propagation of equipment (nodesand links) failures in cascading outages in the power grid.
Deepjyoti Deka and Sriram Vishwanath are with the Department ofElectrical and Computer Engineering, The University of Texas at Austin,Austin, TX 78712 USA. This work is supported by the Defense ThreatReduction Agency (DTRA) through grant
Here, the propagation process begins with an initial failureof a network node or link which leads to a redistributing ofthe power flows for optimal dispatch. This can lead to sometransmission lines running above their prescribed capacitiesand subsequently tripping due to overheating. Subsequentlyin [7], the authors incorporate recovery mechanism for thetripped transmission lines into the failure model and showthat the size of blackout has a power-law distribution as seenin reality. Reference [8] analyzes the problem of finding theoptimal k lines in the grid that can be used for interdictionin the grid to create failures. Power Flow based analysis hasbeen done to analyze geographically correlated failures in [9].Similarly, an interdiction based analysis on grid resilienceconsidering short term impacts is discussed in [25]. However,these models include solving a Optimal Power Flow (OPF)or similar optimization problem to study the propagation offailures. Such an approach is harder to analyze. In particular,it seldom leads to a closed form expression of a metric orparameter that can quantify the resilience of the power gridto failures. In a separate line of work, efforts have been madeto study probabilistic failure propagation in power grids andrelated networks using techniques from percolation theoryand random graph theory. In this approach, initial failuresare supposed to propagate probabilistically from source toneighboring nodes and edges in the grid graph before termi-nating. The final state will often include greater number offailures than the initial state and efforts are made to studythe effect of the grid structure in influencing the spread.References [10], [30] study the effect of removing nodesfrom power grid graphs based on their centrality and degreemeasures and its effect on network connectivity. The authors of[11] analyze the propagation of structural failures in complexrandom networks based on similar neighborhood propagationrules. This approach has been extended to study interdependentnetworks failures as well [12] where nodes of two differentnetworks depend on one another for survivability. Similarly,references [14], [15] have analyzed node and edge percolationbased techniques to understand failure propagation in randomgraphs generated by stochastic geometry. An interaction graphbased model is presented in [19] where power flow basedcascading information is used to generate an interaction graphfor the network to study the inter-nodal dependencies oncascades. A good review of works pertaining to grid resilienceto natural disasters can be found in [26].It is worth mentioning that parallel analytical techniques arealso used in studying social, biological and cyber-networks forinformation dissemination and spread of viruses [20], [21].However, the accuracy of percolation based techniques and further of the use of random graphs in modeling power gridcascading failures is debatable [16]. Existing work [22], [23],[24] has demonstrated that the structure of real power gridsas well as their finite sizes create significant deviations inobserved graph parameters from those predicted in randomgraphs. This is because popular random network models likeErdos-Renyi, Barabasi-Albert, small world and configurationmodels [13] do not accurately capture the specific natureof the spatio-temporal evolution of power grids. Further,sharp breakdown thresholds emerge in analysis of topologicalfailure models on random graphs that are seldom observedin simulations of failures on real grid graphs and IEEE testcases [27]. Such thresholds arise due to the absence of localloops and locally tree-like nature of random graph modelsthat are not encountered in real grids. Hence, it is fair tosuggest that analysis of random graphs to study the failurepropagation (both probabilistic and power flow based) willnot extend directly to real grids.In this work, we focus on power grid failures induced bylarge natural disasters like hurricanes and earthquakes thatcreate disconnected islands in the grid and loss of connectivity.We study the size of the largest connected component in thepost-disaster grid and provide justification for using this asa valid metric for grid vulnerability in modern power gridsand micro-grids that have non-trivial fraction of renewableand other distributed generation resources. Note that priorliterature includes the use of the largest connected graphcomponent in simulation based studies of grid failures [18].This is distinct from failure propagation models where afailed node is considered to affect neighboring nodes with adegree or physical characteristic based probability. We extendprobabilistic analysis previously used for random graphs onknown real grid graphs and popular IEEE test cases to de-termine computable graphical parameters (Eg. eigenvalues ofthe grid adjacency matrix) that can be used to quantize theresilience of grids to such natural disasters. More importantly,we present a modified graph construction based on the truegrid graph and use it to develop improved bounds on theextent of damage created in the network by the disaster [1].The efficacy of the graphical parameter based bounds andsoft thresholds are demonstrated by simulations of failures onpublicly available grid data-sets. We then use the graphicalmetrics on grid resilience to identify critical transmission lines(graph edges) that maximally affect the grid resilience. Inparticular, we study attack on grid resilience by an adversarythat aims to damage a set of transmission lines to maximize theexpected damage to network connectivity following a naturaldisaster. As this problem is NP-hard in general, we presenttwo approximate algorithms to determine the optimal edges inthe adversary’s target set. The first algorithm is based on per-turbation based analysis of eigen-values of the grid adjacencymatrix while the second algorithm is based on greedy traceminimization of a higher power of the adjacency matrix. Theperformance of our algorithms for attack design in reducinggrid resilience is demonstrated through simulations and alsocompared with other techniques in literature, notably attackson nodes with high betweenness [30] or random attacks.From the system operator or grid controller’s perspective, these algorithms can be used to determine the critical linesthat need to be protected to build resilience and preventfurther degradation of grid resilience before any impendingnatural disaster. To summarize, our work presents a analyticalframework to quantify the resilience of real power grid graphsto natural disasters and develops two algorithms to determinethe critical transmission lines that need to be protected toimprove grid resilience and prevent adversarial deterioration.The rest of the paper is organized as follows. In the nextsection, we develop our intuitive graph theoretic quantificationof network resilience that is reasonable in the presence of localgeneration. Next, we analyze network failures and resiliencein actual grid graphs without employing any assumptions fromrandom graph theory in Section III. In Section IV, we presentour novel modified graph construction and use it to developimproved bounds on size of the network damage along withsimulation results on IEEE test cases and real power grids.We study adversarial attacks on transmission lines aimed atweakening grid resilience to natural disasters in Section Vand present our approximate greedy methods to determine thecritical transmission lines. Simulation results on our designedalgorithms and comparison with existing work is presentedin Section VI. Finally, we discuss the insights gained andprospective future work in Section VII.II. F AILURE M ODEL IN P OWER G RIDS
We begin by describing the power grid model and itsfeatures.
Network Model:
We consider a modern power grid (or mi-cro grid) in this paper that has distributed generation resourcesavailable on interior buses. Such generation may be providedby renewables (solar, wind etc.) or by conventional resources.We denote the grid by a graph G = ( V, E ) , where sets V and E represent the nodes/buses and the undirected edges/linesrespectively. Let the total number of buses in the system be N .We assume that under normal operating conditions, the lineshave sufficient transmission capacity to transfer power fromone part of the network to another. We denote the adjacencymatrix of the graph G by A G that is assumed to be known andnot generated by a probabilistic model. Each edge ( ij ) in E is represented by a value of for A ( i, j ) and A ( j, i ) in thebinary adjacency matrix. As an example, the IEEE bus testsystem [27] is given in Figure 1. Failure Model:
As described in the Introduction, we con-sider a natural disaster that causes equipment failures in thegrid. We assume that the natural disaster produces a proba-bilistic failure on all nodes in the system, with independent initial probability of failure denoted by p . In this work, weonly consider cases where the initial probability of failure onall nodes is the same. However, the entire analysis can beextended to cases where different nodes suffer distinct prob-abilities of initial failure. The initial failure rate p dependson the nature of the natural disaster (Eg. earthquake scale,wind speed of hurricane etc.) as well as on the geographicalplacement of the grid (Eg. topography of the land will affectthe failure rate). Such probabilities, in practice, are computedby agencies like the National Hurricane Center and used to Fig. 1. IEEE 14-bus test system [27]Fig. 2. Initial Node Failure ( j ) and Secondary Node Failure ( i ) predict the scale of damage and help in planning for evacuationstrategies [28]. As we are concerned with connectivity in thenetwork, we consider secondary failures in surviving nodesthat get separated from the rest of the network due to initialfailures in all of their neighboring nodes. This is shown inFigure 2. Next we describe our measure of network damagefollowing the natural disaster. Network Damage : As mentioned earlier, we assume thattransmission line capacities are sufficient to satisfy all load inthe system provided enough generating resources are onlineand connected. Let N s denote the size of the largest connectedcomponent of surviving nodes after the disaster. We consider N − N s to be the measure of network damage caused bythe natural disaster. It is worth mentioning that outside of thelargest component, smaller groups of nodes can be functionalas well if enough generation is available to satisfy the totalload in the node group. We select N − N s as the measureof network damage as it has a key characteristic as describednext.Let ∆ Pi denote the net power capacity (generating capacityminus load) of node i in the network, where nodal generatingcapacity and load are both random variables. Consider the casewhere each ∆ Pi is an independent Gaussian random variablewith mean µ ≥ and variance σ . If node i is disconnectedfrom the rest of the network, the probability q i that its localload is served is then given by: q i = P (∆ Pi ≥
0) = Z ∞ √ πσ e − ( x − µ ) /σ dx ⇒ q i = . Z µ √ πσ e − ( x − µ ) /σ dx as µ ≥ (1) On the other-hand, if node i is connected to a group of N k nodes, the probability q iN k that node i satisfies its load is givenby: q iN k = P ( N k X i =1 ∆ Pi ≥
0) = P ( N k X i =1 ∆ Pi /N k ≥ . Z µ p πσ /N k e − ( x − µ )2 σ /Nk dx (2)where Eq. (2) follows from the fact that P N k i =1 ∆ Pi /N k is aGaussian ( µ, σ /N k ) random variable. Note that its variancedecreases with increase in N k , the size of the connected setthat node i belongs to. For ≤ N k ≤ N s , using properties ofthe exponential function it follows that Z µ p πσ /N k e − ( x − µ )2 σ /Nk dx ≤ Z µ p πσ /N s e − ( x − µ )2 σ /Ns dx (3)Using Eqs. (1), (2) and (3), we have q i ≤ q iN k ≤ q iN s .Thus, the probability of a nodal load being served increaseswith an increase in the size of the connected component thatthe node belongs to. Thus the largest component shows thehighest group of nodes whose cumulative loads are satisfiedwith highest probability. This justifies our usage of N s (sizeof the largest connected component of surviving nodes) toquantify the functional network and correspondingly of N − N s to measure the scale of network damage. In the next section,considering the largest component as the surviving network,we analyze the effects of network structure on the extentof failures and determine a preliminary upper bound on theprobability of initial failure p beyond which the networkfragments.III. F AILURE A NALYSIS AND P RELIMINARY B OUND
As shown in Fig. 2, we consider initial failures and sec-ondary failures in the grid and analyze their creation in discretesteps. Let λ Vt denote the vector of survival probabilities of all N nodes in set V at step t . For node i , we have λ V ( i ) = 1 − p where p is the initial probability of failure. According to thefailure model, node i survives at step t if it did not fail at step t = 0 and did not get disconnected between steps and t − .In other words, at least one of its neighbors did not fail bystep t − . We express this mathematically as λ Vt ( i ) = (1 − p ) P [ [ j :( ij ) ∈ E { node j survives at t − } ] ⇒ λ Vt ( i ) ≤ (1 − p ) X j :( ij ) ∈ E λ Vt − ( j ) (4) ⇒ λ Vt ≤ (1 − p ) A G λ Vt − (5)Here Eq. (4) follows from the Union Bound for probabilities.Note that for a general graph G , this gives an inequality asagainst an equality that is obtained for a random graph model[12] where failure propagation from each distinct neighboris independent. This is a crucial distinction as real worldpower grid graphs are not locally tree-like and have correlatedfailure pathways unlike random graphs. Let β A be the largest Fig. 3. Each undirected edge results in two survival probabilities, one ineach direction eigenvalue of the adjacency matrix A G . Using relation (5), wehave (1 − p ) β A < , then λ V ∞ → (6)Thus, p > − /β A provides a upper bound on the thresholdon initial probability of random failures ( p ) beyond which thegrid fragments. In contrast, random graph analysis leads to anexact threshold and not an upper bound. It is worth noting thatthe current formulation does not specify the extent of damagein the region p < − /β A . In the next section, we presenta novel modified graph construction that overcomes this andhelps generate tighter bounds.IV. M ODIFIED G RAPH FOR I MPROVED B OUNDS
Note that in our failure model, the survival of any nodedepends on the existence of edges connecting it to the largestconnected component. Indeed we can analyze a node’s sur-vivability by considering the probability of it being connectedthrough operational edges in the grid graph. To motivate thisapproach better, consider two connected neighboring nodes i and j as shown in Fig. 3. Let B Et ( ij ) be the event that node i is connected to the surviving nodes in the largest componentthrough edge ( ij ) at step t . Let the probability of B Et ( ij ) be denoted by λ Et ( ij ) . Note that this event can be defined forevery neighboring node of node i and i ’s survivability requiresat least one event to be true. Thus, the probability of node i surviving at step t is given by: λ Vt ( i ) = (1 − p ) P [ [ j :( ij ) ∈ E [ B Et ( ij )] (7)Here, the (1 − p ) arises from the probability of node i sur-viving an initial failure, while the remaining terms correspondto the survivability due to a connecting edge. In a similarway, we define event B Et ( ji ) of probability λ Et ( ji ) for node j surviving through the edge with node i . This event is thereciprocal of B Et ( ij ) . Thus, every edge gives rise to twoprobabilities of survival, one along each direction as shownin Figure 3. If nodes i , j and k are connected as shown inFig. 3, the event B Et ( ij ) (node i surviving via edge ( ij ) )depends on B Et ( jk ) (node j surviving through edge ( jk ) ). Interms of their probabilities, λ Et ( ij ) and λ Et − ( jk ) are related.Extending to other nodes in the system, we write this relation Fig. 4. Modified graph G E formation from G mathematically ∀ i, j such that ( ij ) ∈ E as λ Et ( ij ) = (1 − p ) P [ [ k :( jk ) ∈ E,k = i B Et − ( jk )] ⇒ λ Et ( ij ) ≤ (1 − p ) X k :( jk ) ∈ E,k = i λ Et − ( jk ) (Union Bound)(8)In Eq. (8), the (1 − p ) comes from the fact that if node j failsinitially, then node i cannot survive through edge ( ij ) . Eq. (8)is the motivation behind our construction of a modified graphbased on the power grid graph to better estimate the scope ofnetwork failures. Modified Graph Construction:
We create modified di-rected graph G E from the original graph G as follows • Each edge ( ij ) in E gives rise to two nodes v ij and v ji in G E . • Two nodes v ab and v cd in G E are connected by an edgedirected from v ab towards v cd if b = c and a = d (seeFigure 4).Note that G E is similar in structure to a line graph of gridgraph G . However, there is a slight difference in that G E haslesser number of edges as a pair of nodes v ab and v ba are not neighbors in G E , though they would be in a standard line graphconstruction. The number of nodes in G E is equal to twice thenumber of edges in G . The size of A E , the adjacency matrixof G E , is | E | × | E | . We now write Eq. (8) in a vector formsimilar to Eq. (5) using the adjacency matrix of the modifiedgraph G E as follows: λ Et ≤ (1 − p ) A E λ Et − (9)Here each node in the modified graph is associated with oneprobability of failure that represents directional connectivityby the corresponding edge in the original graph G . Extendingthe analysis in the previous section to Eq. (9), it follows that if β A E (the largest eigenvalue of A E ) satisfies p > − /β A E ,then the grid disintegrates as λ E ∞ → and λ V ∞ → . Thus, − /β A E provides a second upper bound on the thresholdon p , beyond which the grid fragments. We compare the firstupper bound − /β A in (6) and the second lower bound − /β A E of IEEE test-cases and real grids and observe thatin all cases the second upper bound is smaller in magnitudeand hence provides an improved tighter bound on the failurethreshold. The comparisons are noted in simulation resultslater in this section.We now use the modified graph G E and Eq. (7) to analyzethe extent of damage in the grid. In particular, we are interested in bounding the size of number of failures over the entirerange of initial probability of failure p , even below the boundsderived earlier. For this, we use the basic failure definitionwhere a node fails eventually if either it fails initially or if itdoes not get connected to the largest component through anyof its edges. In other words, the probability of not surviving(given by − λ Vt ( ij ) for node i in graph G ) depends on theprobability of initial failure and on the probability that none of the neighbors of node v ij in G E survives. Let [ B Et ( ij )] c denote the event of node i not surviving through edge ( ij ) .Mathematically, we have − λ Vt ( i ) = p + (1 − p ) P [ \ j :( ij ) ∈ E [ B Et − ( ij )] c ] (10) ≤ p + (1 − p ) min j :( ij ) ∈ E (1 − P [ B Et − ( ij )]) (11) ≤ p + (1 − p ) X j :( ij ) ∈ E − λ Et − ( ij ) d i (12)where d i is the number of neighbors of node i in G . Eq. (10)follows from the failure definition where p denotes the initialfailure probability for original node i . Eq. (11) follows from P [ A ∩ B ] ≤ min( P [ A ] , P [ B ]) while Eq. (12) follows fromthe fact that the minimum of a set of numbers is less thantheir average. Likewise, we express ( − λ Et ( ij ) ) for failureprobabilities in the modified graph as below. − λ Et ( ij ) = p + (1 − p ) P [ \ k :( jk ) ∈ E,k = i [ B Et − ( jk )] c ] (13) ≤ p + (1 − p ) X k :( jk ) ∈ E,k = i − λ Et − ( jk ) d ij (14)where d ij is the number of neighbors of node v ij in modifiedgraph G E . Writing it in vector form for t → ∞ , we get − λ E ∞ = p ( I | E | − (1 − p ) D − E A E ) − (15)where I | E | is the identity matrix of dimension | E | (numberof nodes in G E ). D E is the diagonal matrix of node degrees( d ij ) in G E , while A E is its adjacency matrix. The upperbound on the expected number of node failures in the gird( N f ) is then given by: N f = N X i =1 (1 − λ V ∞ ( i )) ≤ p N + (1 − p ) T D − G A G ( − λ E ∞ ) ≤ p N + (1 − p ) N X i =1 X j :( ij ) ∈ E − λ E ∞ ( ij ) d i (16)Using Eqs. (15) and (16), the upper bound on node failurescan be computed. If there are nodes in G with unit degree(a common observation in power grids), the upper bound onnumber of failures will be non-trivial. To demonstrate theperformance of our bounds, we present simulations of randomfailures on known power grid graphs. A. Comparison of Bounds through simulations
The performances of the upper bound on network failureover all values of p and the two upper bounds on critical value of p beyond which the network disintegrates are shownthrough simulations on the IEEE and bus test systems[27] in Figs. 5 and 6 respectively. Subsequently we alsoconsider publicly available power grid topologies pertaining tothe Western US grid and the grid under the Union for Coordi-nation of Transmission of Electricity (UCTE) in Europe. TheWestern US grid has nodes and edges [31] whilethe power grid of the UCTE has buses and lines[32]. Failure propagation simulation and determined boundsfor these networks are shown in Figs. 7 and 8. e x pe c t ed nu m be r o f node s f a il ed Uppper BoundSimulation1−1/ β A E β A N −2*log(N)
Fig. 5. Upper bounds for number of failed nodes and bounds on p in IEEE bus system e x pe c t ed nu m be r o f node s f a il ed Upper boundSimulation1−1/ β A E
1− 1/ β A N −2* log(N)
Fig. 6. Upper bounds for number of failed nodes and bounds on p in IEEE bus system In all the power grid cases, we consider the grid to havefragmented if the size of the largest connected component ofsurviving nodes is less than N , where N is the initial sizeof the network. Note that the second upper bound ( − /β A E )on the threshold on p (initial probability of failure) derivedfrom the modified graph is lower and hence tighter than thefirst upper bound ( − /β A ) derived using the original graph.Further, it needs to be pointed out that the upper bound on thenumber of failures given by Eqs. (15) and (16) is tighter forhigher values of p compared to smaller values. The original e x pe c t ed nu m be r o f node s f a il ed Upper bound1−1/ β A E β A N − 2*log(N)Simulation
Fig. 7. Upper bounds for number of failed nodes and bounds on p inWestern US power grid e x pe c t ed nu m be r o f node s f a il ed Upper bound1−1/ β A E β A N − 2*log(N)
Fig. 8. Upper bounds for number of failed nodes and bounds on p in UCTEpower grid graph does not provide a non-trivial upper bound on failuresover the entire range of p , hence the modified graph has twodistinct advantages in failure analysis.In the next section, we discuss how our measures of networkreliability based on eigenvalues can be used to determinecritical transmission lines that may be attacked by an adversaryinterested in weakening the grid resilience to natural disasters.V. C RITICAL L INES FOR A DVERSARIAL A TTACK ON G RID R ESILIENCE
We consider an adversary that aims to maximally weakenthe grid structure to make it more vulnerable to failures duringnatural disasters. The adversary does so by attacking andremoving a fixed number ( k max ) of transmission lines inthe grid. As mentioned in prior sections, relations involvingthe eigenvalues of adjacency matrix of the grid graph or themodified graph provide upper bounds on the probability offailure beyond which grid connectivity diminishes greatly. Inthe remainder of this section, we focus our attention on the first upper bound (1 − /β A ) in Eq. (6), where β A is the largesteigenvalue of the adjacency matrix A G of the grid graph G .We thus formulate the adversary’s objective as damaging k max edges in the graph to minimize β A to reduce the first upperbound. Techniques based on the second upper bound based onthe modified graph in Eq. (9) will be the focus of our futurework in this area.Let normalized eigenvector u correspond to the largesteigenvalue β A of adjacency matrix A G . By definition u satisfies max k x k =1 x T A G x = u T A G u = β A (17)The Perron-Frobenius theorem [34] states that eigenvector u is a positive vector. Using this fact with Eq. (17), it is clear thatremoving edges from graph G (or deleting s from A G ) alwaysleads to a reduction in the magnitude of its largest eigenvalue.Hence, the adversary’s attack consists of determining thecritical edges that will enable the maximal reduction in thelargest eigenvalue of A G . This is a known NP-hard problem.Here, we present two approximate techniques to determine thecritical lines that will be included in the adversary’s target set.
Eigen-Perturbation based Attack Design:
In this attackscheme, we use perturbation analysis [34] to approximatethe change in eigenvalues of adjacency matrix A G followingremoval of edges and subsequently to determine the optimaltransmission lines to attack. Let the new adjacency matrix ofthe grid after removal of edges be given by A G − ∆ A G . Thechange in adjacency matrix ∆ A G has the following structure: ∆ A G ( i, j ) = ( if edge ( ij ) ∈ E is removed , otherwise, (18)Let the largest eigenvalue of new adjacency matrix be β ∆ A .From Eq. (17), we have β ∆ A = max k x k =1 x T ( A G − ∆ A G ) x ⇒ β ∆ A ≥ u T ( A G − ∆ A G ) u = β A − u T ∆ A G u ⇒ ∆ β A = β A − β ∆ A ≈ X removed ( ij ) u ( i ) u ( j )( using (18)) (19)where ∆ β A denotes the change in the maximum eigen-valuefollowing the removal of lines. The optimal transmission line whose removal approx-imately minimizes the maximum eigenvalue of the adja-cency matrix is thus given by maximizing the expression inEq. (19).
To determine the optimal k max lines, the adversaryiteratively computes u , the eigenvector corresponding to thelargest eigenvalue, removes line by maximizing Eq. (18) andrecomputes the adjacency matrix and its principal eigenvector.As shown later, the iterative scheme provides a far greaterreduction in grid resilience than selecting the k max lines tomaximize Eq. (18) all at once, though it leads to an increasein computational complexity. Complexity:
The computation of the eigenvector u takes O ( N ) steps via Singular Value Decomposition (SVD) of A G [33]. Given that the maximization of Eq. (19) takes | E | stepswhich is less than O ( N ) , computing the k max critical lines by iteratively computing the eigenvector for the largest eigenvaluehas a complexity of O ( k max N ) [33]. If all lines are selectedbased on a single computation of u , the complexity is O ( N ) due to SVD. Next we describe another technique for attackdesign that depends on trace minimization. Trace Minimization based Attack Design:
The trace ofa matrix refers to the sum of its diagonal elements and isequal to the sum of its eigenvalues [34]. Consider an even r th power of the trace of the adjacency matrix A G . As theeigenvalues of A r G are the r th powers of the eigenvalues( β = β A , β , ..., β N ) of A G , we have the following relationfor the trace trace ( A r G ) = N X i =1 β ri = β rA (1 + ( β β A ) r + ... + ( β N β A ) r ) (20) ⇒ trace ( A r G ) /β rA ≈ as r → ∞ (21)Here β A = β i is the largest eigenvalue of the adjacency matrix.Note that as | β β A | < , if we take higher values of r , the ratio ofeigenvalues becomes smaller in Eq. (20). Thus for extremelylarge values of r , the largest eigenvalue and trace of the r th power of A G are approximately equal as noted in Eq. (21). Inthis approach, thus we focus on reducing the trace of A r G byremoving lines instead of minimizing the largest eigenvalue β A or its higher power. Finding the optimal set of k max edgesto minimize the trace of A r G is computationally hard as well,however using trace minimization has certain advantages aswe discuss now. Theorem 1.
The trace of A r G , where A G is the adjacencymatrix of grid graph G is a supermodular function of theconstituent edges in the graph. A real-valued function f defined over set S is supermodular[35] if f ( A S C ) ≥ f ( B S C ) for B ⊂ A and A, B, C aresubsets of S . In other words, the returns due to addition of C are not diminishing. Proof:
Note that the i th diagonal element in A r G is equalto the number of cycles of length r that begin and end at node i . This can be shown by direct checks or by mathematicalinduction. Here, cycle of length r refers to a graph pathwith r hops (repetition allowed) that begins and ends at thesame node. Thus, the trace (sum of the diagonal elementsof A r G ) is given by the total number of cycles of length r that can be formed on all nodes in the grid graph. To showsupermodularity of trace of A r G as a function of graph edges,it is sufficient to show that the increase in the number of cyclesof length r in graph G after adding a new edge ( ij ) is lessthan the increase observed if edge ( ij ) is added after inclusionof another edge ( lm ) . This increase is indeed true as presenceof an edge ( lm ) prior to the addition of edge ( ij ) will permitthe existence of additional cycles that includes both edges ( lm ) and ( ij ) , and cannot exist without ( lm ) . Hence trace of higherpower of the adjacency matrix is a supermodular function ofthe graph edges.It is a known property [35] that greedy minimization of asupermodular function is equivalent to greedy maximizationof a submodular function and is provably at least − /e ( ≈ ) close to the optimal solution. Thus, the adversary’sattack policy in this scheme is to greedily remove k max edgesthat minimizes the trace of r th power of the adjacencymatrix of the grid graph.Complexity: The r th power of the symmetric adjacencymatrix is computed efficiently using Singular Value Decom-position (SVD) as A r G = U β rA U T where columns of U arethe eigenvectors and β rA is the diagonal matrix with r th powers of the eigenvalues. Note that matrix multiplicationand SVD are computed in O ( N ) while computing β rA takescomplexity O ( N log r ) . Since we greedily minimize the trace,the selection of one edge takes O ( | E | ( N + N log r )) . Theoverall complexity of computing k max optimal edges by thisscheme is thus O ( k max | E | ( N + N log r )) . This expressionimplies that increasing r to improve the accuracy of thisapproach will at most lead to a logarithmical increase in thecomplexity. Resilience:
From the grid controller’s perspective, these twotechniques can be used to determine the critical transmissionlines for enhancing security and reinforcement to preventadversarial manipulation aimed at disrupting grid resilienceto natural disasters. In the next section, we look at theperformance of these two approaches as an adversarial tool.VI. S
IMULATION R ESULTS OF A DVERSARIAL A TTACKS
We consider both approaches (eigen perturbation and traceminimization) for determining the optimal k max edges to min-imize the largest eigenvalue of the adjacency matrix of the gridgraph and thereby reduce the resilience of the grid to naturaldisasters. For comparison, we consider two alternate schemes,one where an adversary removes edges randomly, and anotherwhere an adversary removes edges in the decreasing order oftheir betweenness centralities [30]. We plot our results for theIEEE and bus test systems and the UCTE power gridnetwork in Figs. 9, 10 and 11 respectively. Note that bothalgorithms outperform random and betweenness based attacksto reduce the eigenvalues. It can also be noted that iterativeeigen perturbation reduces the largest eigenvalue further thanedge removal based on a single perturbation computation asmentioned in the previous section. Further, it can be observedfrom Figs. 9 and 10 that increasing the value of r , thepower of the adjacency matrix, leads to an improvement inthe trace minimization based scheme as it approximates thelargest eigenvalue better as noted in Eq. (21).VII. C ONCLUSION
We analyze topological vulnerability of power grids toprobabilistic failures introduced by natural disasters in thispaper. We present intuitive evidence that in modern grids andmicro-grids where distributed generation resources are present,a reasonable metric of damage is given by the size of thelargest connected component in the post-event grid graph.We analyze the evolving failure process that originates atnodes with initial failures. Based on the largest eigenvalueof the adjacency matrix of the grid, we present an upperbound on the critical probability of node failures beyond whichthe grid fragments. Further, we present the construction of a La r ge s t e i gen − v a l ue o f ad j a c en cy m a t r i x single eigen−perturbationiterative eigen perturbationtrace min., 4th power trace min., 8th power trace min., 16th power random removalbetweenness based Fig. 9. Comparison of adversarial schemes to reduce maximum eigenvalueof grid adjacency matrix in IEEE bus test system. La r ge s t e i gen − v a l ue o f ad j a c en cy m a t r i x single eigen−perturbationiterative eigen−perturbationtrace min., 4th powertrace min., 8th powertrace min., 16th powerrandom removalbetweenness based Fig. 10. Comparison of adversarial schemes to reduce maximum eigenvalueof grid adjacency matrix in IEEE bus test system. modified graph to analyze the probabilistic failures and use itto generate a tighter upper bound on the critical probability.This modified graph construction also enables us to derivenew non-trivial upper bounds on the expected number oftotal failures for all values of the initial failure probability.We present the performance of our derived analytical boundsthrough simulations on two IEEE test cases and two realgrid data sets. Finally, we discuss adversarial attacks on thepower grid aimed at damaging transmission lines to minimizethe grid’s resilience to natural disasters. We develop twoapproximate algorithms to identify the critical lines that willenable such adversarial attacks. The first algorithm is basedon perturbation analysis of the eigenvalues of the adjacencymatrix and the second algorithm is based on greedy minimiza-tion of the trace of a higher power of the adjacency matrix. Weanalyze both algorithms and their complexity and demonstratetheir performance against random and centrality based attacksstudied in literature through simulations. Potential areas offuture work include improving the bounds and developing aframework to incorporate topological analysis into power flow La r ge s t e i gen − v a l ue o f ad j a c en cy m a t r i x single eigen−perturbationiterative eigen−perturbationtrace min., 4th powerrandom removalbetweenness based Fig. 11. Comparison of adversarial schemes to reduce maximum eigenvalueof grid adjacency matrix in UCTE power grid. based studies on grid vulnerability to enhance its practicalcontribution. R
EFERENCES[1] D. Deka, R. Baldick, and S. Vishwanath, “Structural Vulnerability ofPower Grids to Disasters: Bounds and Reinforcement Measures”,
IEEEConf. on Innovative Smart Grid Tech. , 2015.[2] J. Lavaei and S. Low, “Zero duality gap in optimal power flow problem”,
IEEE Trans. Power Systems , vol. 27, 2012.[3] W. Chi-Keung, I. Horowitz, J. Moore, and A. Pacheco, “The impact ofwind generation on the electricity spot-market price level and variance:The Texas experience”,
Energy Policy , vol. 39, 2011.[4] A. Kwasinski, W. W. Weaver, P. L. Chapman, and P.T. Krein, “Telecom-munications Power Plant Damage Assessment for Hurricane Katrina SiteSurvey and Follow-Up Result”,
IEEE Systems Journal , vol.3, 2009.[5] L. Kantha, “Time to Replace the Saffir-Simpson Hurricane Scale?”,
Eos ,vol. 87, 2006.[6] I. Dobson, B.A Carreras, V.E. Lynch,and D.E. Newman, “An initial modelfo complex dynamics in electric power system blackouts”,
Proc. 34thHawaii Int. Conf on System Sciences , 2001.[7] I. Dobson, B.A Carreras, V.E. Lynch,and D.E. Newman, “Complexsystems analysis of series of blackouts: Cascading failure, critical points,and self-organization”,
Chaos , vol. 17, 2007.[8] D. Bienstock and A. Verma, “The N k problem in power grids: Newmodels, formulations, and numerical experiments”,
SIAM J. Optim. , vol.20, 2010.[9] A. Bernstein, D. Hay, M. Uzunoglu, and G. Zussman, “PowerGrid Vulnerability to Geographically Correlated Failures: Analy-sis and Coontrol Implications”, arxiv preprint , 2012. Avaible at:http://arxiv.org/abs/1206.1099.[10] R. Albert, H. Jeong, and A.-L. Barabasi, “Error and attack tolerance ofcomplex networks”,
Nature , vol. 406, 2000.[11] J. Wang and L. Rong, “Cascade-based Attack Vulnerability on the USPower Grid”,
Safety Science , volume 47, 2009.[12] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin,“Catastrophic cascade of failures in interdependent networks”,
Nature ,vol. 464, 2010.[13] G. A. Pagani and M. Aiello, “The power grid as a complex network: asurvey”,
Physica A: Statistical Mechanics and its Applications , 2013.[14] Z. Kong and E. M. Yeh, “Resilience to degree-dependent and cascadingnode failures in random geometric networks”,
IEEE Trans. Inf. Theory ,vol. 56, 2010.[15] H. Xiao and E. M. Yeh, “Cascading link failure in the power grid:A percolation-based analysis”
Proc. IEEE Int. Work. on Smart GridCommunication , 2011.[16] P. Hines, E. Cotilla-Sanchez, and S. Blumsack, “Do topological modelsprovide good information about electricity infrastructure vulenrablity?”,
Chaos , vol. 20, no. 3, p. 033122, Sept. 2010.[17] R. Durrett, “Random graph dynamics”,
Cambridge University Press ,2006. [18] C. M. Schneider, A. A. Moreira, J. S. Andrade, S. Havlin, and H. J.Herrmann, “Mitigation of malicious attacks on networks”,
Proc. NationalAcademy of Sciences , 2011.[19] J. Qi, K. Sun, and S. Mei, “An interaction model for simulation andmitigation of cascading failures”,
IEEE Trans. Power Systems , 2015.[20] J. C. Miller and J. M. Hyman, “Effective vaccination strategies forrealistic social networks”,
Physica A , vol. 386, 2007.[21] M. Boguna, R. Pastor-Satorras and A. Vespignan, “Epidemicspreading in complex networks with degree correlations”, arXiv:cond-mat/0301149v1 [cond-mat.stat-mech] , 2003.[22] Z. Wang, A. Scaglione, and R. Thomas, “Generating statistically correctrandom topologies for testing smart grid communication and controlnetworks”,
IEEE Trans. Smart Grid , vol. 1, 2010.[23] D. Deka and S. Vishwanath, “Generative Growth Model for PowerGrids”,
Int. Conf. on Signal-Image Tech and Internet-Based Systems ,2013.[24] D. Deka and S. Vishwanath, “Analytical Models for Power Networks:The case of the Western US and ERCOT grids”, arxiv preprint , 2015.Available at: http://arxiv.org/abs/1204.0165.[25] Y. Wang and R. Baldick, “Interdiction Analysis of Electric Gridscombining cascading outage and medium-term impacts”,
IEEE Trans.Power Systems , 2014.[26] Y. Wang, C. Chen, J. Wang, and R. Baldick, “Research on Resilience ofPower Systems Under Natural DisastersA Review”,
IEEE Trans. PowerSystems
IEEE INTELEC , 2014.[29] Q. Zhou and J. W. Bialek, “Approximate Model of European Intercon-nected System as a Benchmark System to Study Effects of Cross-BorderTrades”,
IEEE Transactions on Power Systems , Vol. 20, No. 2, May 2005.[30] A. E. Motter and Y. -C. Lai, “Cascade-based attacks on complexnetworks”,
Phys. Rev. E , vol. 66, 2002.[31] D. J. Watts and S. H. Strogatz, “Collective dynamics of ’small-world’networks”,
Nature ∼ jbialek/Europe load flow/.[33] G. H. Golub and C. F. V. Loan, Matrix computations , Vol. 3, JHU Press,2012.[34] G.W. Stewart and J. Sun, “Matrix perturbation theory”, 1990.[35] G. Nemhauser, L. Wolsey, and M. Fisher, “An analysis of the ap-proximations for maximizing submodular set functions”,