Decreased resilience in power grids under dynamically induced vulnerabilities

In this paper, a methodology inspired on bond and site percolation methods is applied to the estimation of the resilience against failures in power grids. Our approach includes vulnerability measures with both dynamical and structural foundations as an attempt to find more insights about the relationships between topology and dynamics in the second-order Kuramoto model on complex networks. As test cases for numerical simulations, we use the real-world topology of the Colombian power transmission system, as well as randomly generated networks with spatial embedding. It is observed that, by focusing the attacks on those dynamical vulnerabilities, the power grid becomes, in general, more prone to reach a state of total blackout, which in the case of node removal procedures it is conditioned by the homogeneity of power distribution in the network.


Introduction
Emergence of synchronization in systems of agents, coupled through local interactions, is a widely studied phenomenon with multiple applications in physics, biology and social sciences [1,2,3]. In particular, the optimal functioning of a power grid, which represents one of the most complex interacting systems in engineering, highly depends on its ability to maintain a synchronous operation over time despite external disturbances. Losing that synchrony, even locally, may lead to cascading failures and complete blackout of the network [4,5,6]. Two main concerns have motivated the discussion around dynamical analysis and design of power grids in the last years, those are: the strong economical and social impact that a power outage could cause in the highly electricity-dependant modern society [7], and the slow transition to renewable energy sources that is being promoted all around the world [8,9], since it is known that renewable and small power producers have negligible inertia, which imposes instability risks in the power grids [10]. Multiple models have been proposed for the analysis of synchronization in power grids, and they differ mainly in the way loads, power generators and transmission lines are modeled [11]. Probably the most common model is the one called synchronous motor, where both, generators and consumers are represented by second-order synchronous machines [12]. This model is essentially equivalent to the celebrated Kuramoto model when inertia terms are considered [13,14,15]. Self-synchronization behavior in this model occurs when the interactions between connected oscillators are sufficiently strong to overcome the dissimilarity in the ensemble [4,16,17,18].
Recently, the scientific community has put a lot of effort in trying to determine the topological or structural features of those interactions that can enhance or undermine non-linear stability of a power grid, that is, its capability to reject finitesize disturbances, a matter of paramount importance in the design of the future smart grids. This has been approached, for instance, by the means of energy barrier functions [19] and the basin stability concept [20,21,22,23,24,25]. Some interesting findings that deserve to be mentioned, as they provide great insights about the relationship between topology and dynamical stability include: the poor basin stability usually detected on dead-tree arrangements [20], the strong stability found on triangle-shaped motifs [26], the enhanced non-linear stability achieved by increasing global redundancy in the connections [27] or by adding small cyclic motifs [19] and the lower basin stability exhibited in general by high-power generator nodes in [28]. Other studies have focused on estimating the resilience of the complex network against cascading failures [29,6,30], the robustness in the transient behavior after localized disturbances [31] and their diffusion through the interconnected system [10].
In the present work, a percolation inspired method is presented in order to evaluate the resilience of a power grid against random and focalized disturbances. By defining a vulnerability measure in terms of the dynamical properties of both, nodes and edges, it is found that, attacks that focus on the weakest components of the network, can easily provoke cascading failures that lead to total blackout faster than random attacks in some specific cases. Note that, in the following, the term total blackout refers to the situation when no operational composition of at least one generator and one consumer exists in the network. Test cases chosen for this study include randomly generated networks with spatial embedding that emulate properties of real-world power grids, as well as the particular case of the Colombian power transmission network. It is also found that, by attacking the most dynamically-loaded transmission lines in the graph, the power system reaches blackout faster than any other bond-percolation criterion, while, for the sitepercolation methods, it is produced by attacking the most dynamically vulnerable nodes, but only if the distribution of power generation and consumption is not homogeneous, which often occurs in real-world architectures.
A similar work can be found in [32], where the authors use centrality measures, namely, the node degree, clustering coefficient and betweenness centrality, to estimate the vulnerability of a node and then proceed to remove weakest nodes to observe the evolution of the giant component in the network as a function of the number of remaining nodes in the whole graph. This structural approach to quantify the node vulnerability will be also used in this work to compare results of removal methods based on dynamical measurements. Other relevant works have to be mentioned; for instance, based on the statistical physics description of networks, some basic percolation properties for different networks have been analysed using generating functions [33], the lifetime and reliability of networks has been estimated by setting some probability distributions for node failures [34] and the propagation of cascading failures through the power transmission or communication architecture of smart grids has also been modeled [35]. The rest of this paper is organized as follows: Section 2 introduces the dynamical model of a power grid and the structure of the specific test cases, as well as some theoretical tools that will be applied in further analysis of synchronization and vulnerability measurements. Section 3 presents an algorithm for resilience assessment of power grids and results regarding line and node removal procedures over synthetic power grids and the Colombian transmission network. Finally, some concluding remarks are discussed in section 4.

Model and methods
In this paper, the framework of the synchronous motor representation is used to model the general structure and dynamics of power grids [11,12], since it has been employed constantly and successfully in the past years, as mentioned in the previous section. This approach is described in the following.

Power grid model
As a simplified model of a real-world power grid, consider a system of N interacting synchronous machines arranged in a connected graph G(ϑ, ε), with a set of nodes ϑ = {1, 2, ..., N } and a set of edges ε ⊂ ϑ × ϑ, such that the amount of edges is |ε| = M . Nodes can be labeled as generator machines (supply energy to the grid) or consumer machines (demand energy from the grid), thus ϑ = ϑ g ∪ ϑ c , being ϑ g the set of generators and ϑ c the set of consumers.
Let the dynamical state of each machine be represented by its phase angle φ i and its phase velocityφ i := dφ i /dt, i ∈ ϑ. Under an appropriate operation of the power grid, it is expected that every machine will be rotating at some reference frequency Ω = 2πf (where by convention f is either 50 Hz or 60 Hz), so let θ i(t) = φ i(t) − Ωt be the phase deviation of the machine i with respect to the reference angle Ωt, then the dynamical behavior of θ i and its velocityθ i can be described by the well-known second-order Kuramoto model for coupled oscillators given by [12,16]: where P i is equal to the injected (drained) power at the i-th node up to a scaling factor, and it is positive (negative) for generators (consumers); α i is related to the damping coefficient of the i-th machine, a ij are the elements of the adjacency matrix A, thus, a ij = 1 if nodes i and j are connected (that is, (i, j) ∈ ε) and a ij = 0 otherwise. The constant K is the maximum power transfer capacity between any pair of connected nodes, which depends on the impedance of transmission lines connecting them. This equation is derived when the ohmic losses on transmission lines are neglected, an equal impedance level is assumed for every transmission line in the network and the deviations from the phase velocity at any node as compared to the reference Ω are small, that is |θ i | Ω, ∀i ∈ ϑ [12,16]. The dynamical system described by equation (1) is said to be in a phase synchronized state if θ i(t) = θ j(t) , ∀i, j ∈ ϑ, and in a frequency synchronized state ifθ i(t) =θ j(t) = ω s , ∀i, j ∈ ϑ. Since we are interested in the dynamics around the synchronization frequency Ω, it can be assumed without loss of generality that ω s = 0, which amounts to consider the dynamics in the co-rotating reference frame Ωt (see [4]). In the following, the topology of G(ϑ, ε) is taken from the Colombian National Transmission System [36], a power grid composed of 158 edges and 102 nodes (28 of them are generators and 74 are consumers). Power supplied by generators is set to P k = 1.0, ∀k ∈ ϑ g , and the power drained by consumers is uniformly distributed such that the power balance condition in the network is fulfilled, that is: This condition is a requirement in order to allow the existence of a locally-stable equilibrium point in the system [16]. Finally, the damping and maximum power transfer are fixed to α i = 0.1 and K = 12.0 respectively. Figure 1(a) shows the topology of the network with the tree-like node classification introduced in [22]. Panels (b)-(e) present distributions for some classic centrality measures of the graph: node degree, clustering coefficient [37], node betweenness and edge betweenness centrality [38]. Panels (f) and (g) display the behavior of θ (t) andθ (t) for every oscillator when the system converges from an unsynchronized state to a frequency synchronized regime given the parameters above. Aiming for a generalization of the results, random test cases were also used and compared to the Colombian test case just described; to do that, synthetic power grids were generated with the algorithm proposed in [39], which is specially useful as it creates spatially embedded networks that emulate connectivity distributions of realworld power grids. Synthetic networks were constructed by initializing a tree-like base architecture and then following a growth procedure conditioned by the optimization of a redundancy-cost function that depends on both, geographic distance and number of paths separating each pair of nodes in the graph. For a detailed description of the algorithm, refer to Appendix A. The parameters used for the growth model algorithm were chosen to match those used on [22], as presented in table 1. Note that the final size of the network is fixed to N = 102 so that it matches the Colombian power grid size, nevertheless, geographic positions for nodes are chosen uniformly at random and an equal amount of generators and consumers is placed to reduce heterogeneity in the samples (thus, P k = 1.0, ∀k ∈ ϑ g and P j = −1.0, ∀j ∈ ϑ c ). Table 1: Parameters for the construction of synthetic power grids.

Parameter
Value N 0 : An initial amount of nodes, N 0 ≥ 1. 1 N : Final amount of nodes that the network will have N > N 0 . 102 y = {y 1 , y 2 , ..., y N }: Geographical location of nodes.

Synchronization measures
A helpful indicator that quantifies the degree of synchronization is the complex-valued order parameter r defined as [40]: Here, Ψ (t) is the average phase of the oscillators at time t. If all phases θ j are identical, the magnitude of the order parameter is r (t) = 1. Conversely, if they are equally distributed around the unit circle r (t) = 0 and the system is said to be desynchronized. Intermediate values of r correspond to partially synchronized states. Additionally, in the case of power grids in which equation (2) holds, when the system is synchronized the average phase is Ψ (t) = 0, and the real part of the order parameter IR[r] = 1, while it oscillates around zero otherwise. While the dependence of the order parameter with time is useful in transient analysis, we will focus on the steady state behavior, so we can define the average order parameter in steady state as [16]: Another measure of frequency synchronization in power grids is the squared average rotational speed v 2 (t) defined as [16]: and the associated steady state value of the speed v ∞ :

Edge vulnerability
Following the results in [4], it can be proven that, by applying a small-angle approximation, the phases at the frequency synchronized equilibrium point, which are denoted by θ * ∈ IR N , can be estimated as: where L † is the pseudo-inverse of the Laplacian matrix, which is computed through the adjacency matrix as L = diag( n j=1 a ij )−A. Let us now define the oriented incidence matrix B = {b ij } ∈ R N ×M , with: Note that although G is an undirected graph, a direction for each edge can be assumed without loss of generality for the following analysis, allowing for this definition of B. Under these circumstances, the steady state phase difference between any pair of connected nodes (i, j) ∈ ε can be approximated by: with ∆ θ ∈ R M . Then, the system (1) is guaranteed to have a unique and stable solution with synchronized frequencies and cohesive phases if the next sufficient condition for ∆ θ holds: with 0 ≤ γ < π/2. In the limit γ → π/2, equation (10) simply reduces to ||∆ θ || ∞ < 1.
In other words, equation (10) states that, in order to achieve synchronization, the largest difference between the phases of connected pairs in the network has to be smaller than π /2 [4]. This rather simple equation, readily connects topological and dynamical features of the network. While condition (10) is only sufficient, it is a rapid way to assess synchronization without relying on the computationally expensive time evolution of a large number of phase oscillators. Moreover, equation (10) provides an approximation to the critical value K c of the coupling strength when substituting from (9): In figure 2(a), it is shown the behavior of the steady state order parameter IR[r ∞ ] as a function of the coupling strength K. As seen in the figure and the insets, the real part of the order parameter oscillates around zero when the coupling strength is small, and therefore its average value is zero. At a critical coupling strength K c ≈ 1.51, the system achieves a state of phase cohesiveness characterized by IR[r ∞ ] ≈ 1 indicating a synchronized regime. The same transition can be observed by looking at the average velocity in figure 2(b). In particular, it can be seen that before the synchronization transition, a finite and different than zero, steady-state value for the angular velocity is obtained, indicating that the phases of the oscillators increase constantly and no equilibrium is reached. After the transition, v (t) abruptly decreases to zero, meaning that the phases of the oscillators remain in a steady state.
In both figures, it is possible to see as well the accuracy of the sufficient condition (11), where the critical value is indicated by the vertical dashed line, showing a good agreement with the numerical value of the transition. This evidence supports the use of condition (10) as a synchronization requisite in the following analysis. Note that the chosen value of K = 12.0 mentioned before is justified here, since K > K c , thus it allows for a synchronized state to exist.
As a measure of the vulnerability of each transmission line, we take advantage of the equation (9), given that the edge (i, j) is considered weak when |∆ θ (i, j)| is close to 1, while it is considered more stable when it is close to 0, regarding the synchronization limit imposed by (10). This approach is equivalent to that presented in [41,42], where the weakness of an edge is measured as the power transferred between the linked nodes (line load), except that here, the steady state of the system is approximated immediately from equation (9).

Nodal vulnerability
Consider a multistable dynamical system that evolves in the state space X and let X * ⊂ X be the set of desirable attracting states. The basin of attraction of X * , noted by β, is defined as the set of all initial conditions x (0) that asymptotically converge to X * [20]. For the purpose of this work, X * will be the frequency synchronization manifold of the system (1). Similarly, the likelihood of a randomly perturbed trajectory to return back to β, known as the basin stability S B , can be defined as [20,21]: where Γ β(x) is a function that indicates whether a state x belongs to the basin of attraction of X * or not, that is The function ρ (x) is the density of states to which the system can be pushed by some non-local perturbation, such that X ρ (x) dx = 1. Note that under this definition, the basin stability is a number between 0 and 1, being S B = 0 when the synchronous state is unstable and S B = 1 when it is globally stable.
Monte Carlo simulations can be performed in order to estimate S B by randomly sampling a sufficiently high number of disturbed states I C (following a certain distribution ρ (x) ) from a representative subspace Π ⊂ X, and using them as the initial conditions for time-evolution simulations. For this work, ρ (x) is chosen as an uniform distribution in the restricted subspace Π, that is: Finally, the amount F C of simulated trajectories that are found to approach asymptotically to the attractor is assumed to be proportional to the volume of the basin of attraction (restricted to the subspace Π), thus S B is approximated by S B ≈ F C /I C . Now, making use of the concept of basin stability we can define an indicator of the robustness of a node against large perturbations applied to it, by means of the singlenode basin stability (SNBS) [21,23], given by: where the I C initial conditions are drawn from perturbations applied to node i. As illustration, figure 3 shows the phase plane of some node i subject to large disturbances; its SNBS is, loosely speaking, the number of crosses divided by the number of total datapoints. Throughout this work, the disturbances for each node are drawn from the subspace Π = [−π, π] × [−100, 100].

Structural vulnerability
Equations (9) and (15) define already a dynamical vulnerability for edges and nodes, respectively. Furthermore a vulnerability based merely on the connectivity features of the complex network could also be considered, since it is intuitive that by removing a highly connected node or a certain critically located edge could lead to a rapid diminishing on power grid stability and performance. In that regard, let us define structural vulnerability in terms of the following centrality measures: k : Amount of nodes to which node i is connected, normalized by the maximum possible degree in the network. It can be computed as: • Clustering coefficient c (i) k : Number of triangles (T i ) in which node i is involved, normalized by the maximum possible amount of such triangles [37], that is: Figure 3: Phase space of a node i of the Colombian grid when random disturbances are applied to θ i orθ i . Green crosses are initial conditions that return to the frequency synchronized state, while the blue circles are out of the basin of attraction. The red dashed line denotes the subspace Π from which disturbances are chosen. This diagram was created by sampling I C = 500 points.
• Node betweenness centrality b (i) k : Sum of the shortest paths between every pair of nodes (s, t) in the network that pass through node i: where σ s,t is the amount of shortest paths between nodes s and t and σ s,t (i) is the number of those paths that pass through node i [38].
• Edge betweenness centrality e (l) k : In the same fashion as b k measures centrality for a node, e k does it for an edge; it is simply defined as the number of shortest paths in the network that include edge l: So in this work, a node will be considered structurally weak if d k , c k or b k is high. Similarly, an edge will be considered structurally weak if its e k is high. (computed with I C = 50) and the size of the node is proportional to the power P i . The red arrow indicates the node that is going to be suppressed. Also, the inset shows the histogram of the basin stability for the whole network. One realization is presented in panels (a)-(d) and a different one in panels (e)-(h) in order to visualize two different mechanisms of cluster vanishing.

Resilience measures
In order to assess the resilience of a power grid, we propose a percolation-based algorithm. To do so, we will follow the capability of the network to maintain its functioning upon different node removal (site percolation) and edge removal (bond percolation) rules. In particular, the algorithm that will be described below removes a node or an edge from the graph on each iteration by applying a random attack, where the element to be removed is chosen uniformly at random; or a focal attack, where the most vulnerable node or edge in the graph is removed first. Here, the most vulnerable node is the one with the lowest basin stability (S (i) B ) or the highest degree centrality (d k ), clustering coefficient (c k ) or node betweenness centrality (b k ); while the most vulnerable edge is the one for which phase difference (∆ θ (i, j)) or edge betweenness centrality (e k ) is the highest.
The percolation-based algorithm then proceeds as such: (i) Remove a node (edge) by applying a random or a focal attack. Note that in a focal attack, multiple nodes (edges) could share the same vulnerability level, in that case, the attacked node is chosen randomly among them.
(ii) Compute the existing disconnected clusters in the new graph.
(iii) For each cluster, check if it contains at least one generator node and one consumer node, if it does not, remove the whole cluster from the graph.
(iv) For each remaining cluster the power has to be compensated, so as a way to simulate a real-world power control (where the supplied power level is increased or reduced according to the demand of the consumers), the power of all generators in the cluster is modified uniformly such that equation (2) remains true.
(v) For each remaining cluster, check whether or not it is synchronizable, that is, its topology satisfies the synchronization condition that K > K c , where K c is computed by the approximation (11). If it does not, the whole cluster is removed from the graph.
(vi) Measure the needed observables and return to step 1 until the whole power grid has gone to blackout, that is, no operational cluster exists in the network.
An important remark has to be done about the algorithm, regarding node removal under the focal attack methodology based on SNBS. Computing this vulnerability index is computationally expensive since it requires a large amount of time-domain simulations in order to estimate reasonably well the equation (15). Thus, in the following, for the S B -focused percolation, only 10 simulations are considered for averaging purposes, while for each of the other cases, 500 simulations are performed. It is also worth noting that the fine structure of the basin of attraction is not captured by this approach, and it could impose numerical problems on dynamical systems that posses either fractal basin boundaries or riddled or intermingled basins [43].
For this elimination procedure, two main phase transitions will be analysed: • Transition T 1 : This transition occurs at the iteration when the graph is divided into multiple functional subgraphs for the first time (turning from a connected graph to an unconnected graph).
• Transition T 2 : The final iteration of the algorithm, when the power grid goes to complete blackout. Figure 4 illustrates the algorithm for nodal resilience measurement under focal attacks. The transition from the frame (a) to (b) shows the formation of 2 independent clusters after removing one node. Note that in the frame (c) there is only one generator in the northern cluster, and interestingly, it has a very poor basin stability level, so when that node is removed from (c) to (d), the whole northern cluster goes to blackout. A different realization is presented in panels (e)-(h); note that from (g) to (h) it suffices to remove one node to bring the bigger cluster to blackout even though there would still remain connected subgraphs with both generators and consumers. This occurs because the remaining cluster would have a topology that does not sattisfy the synchronization condition (10) and thus it vanishes. Similarly, figure 5 shows the algorithm in action under ∆ θ -focused attacks. In the transition from panel (b) to (c), it can be appreciated that, after suppressing a very vulnerable transmission line, a significant cluster of isolated consumers in the east side of the network vanishes.  Figure 5: Edge removal algorithm under the focal attacking scheme. Edge color is mapped with the ∆ θ of the nodes it connects. The red arrow indicates the edge that is going to be supressed. Also, the inset shows the histogram of the phase differences for the whole network.
3.1. Node removal resilience for the Colombian power grid It is thus observed that the grid goes to complete blackout around 30th removal iteration for the case of the focal attacking on S B , and around the 62nd iteration for the random attacking. For d k and b k the transition occurs at the 40th iteration and for the random and c k cases, we have the less effective removal technique since this transition occurs around the 60th attack. The fact that attacking the most vulnerable nodes based on S B produces a faster transition to the total blackout state of the network than just removing nodes at random implies that there is indeed a relationship between the stability of a node against dynamical disturbances (related to the SNBS) and the structural connectivity and resilience to cascading failures of the graph.
Furthermore, by focusing the attacks on nodes with the highest b k or d k , the dimension of the giant component goes down drastically during the first elimination iterations and then it goes down at a slower rate, reaching total blackout in an intermediate level of attacks between the S B and random methodologies. c k on the other hand, produces a slower reduction of the graph size and a total blackout with the most iterations in average. This is related to the tree-like structure of the power grid that can be seen on figure 1. In this kind of topologies, the nodes with the higher degree and betweenness are located at the bulk regions of the graph, therefore, their removal likely divides the giant component instantly into multiple smaller clusters, as also confirmed by the figure 6(d), where it is shown that these methodologies produce by far, the highest     amount of operational clusters (quantity denoted by |S G |) during the whole procedure. Correspondingly, a node with a high clustering coefficient, by definition, implies that the neighbours that it is connected to, are also connected between them, hence, removing it is not going to divide the graph and generate new clusters. Eventually, after multiple node eliminations, redundant links between said neighbours will be suppressed and c k will become uniform over the whole system. That is the reason why this attacking strategy behaves roughly in the same way as the random attacking one after a few iterations. Figure 6(c) plots the size of the second largest component of the graph C 2 , which reveals the location of the phase transition T 1 as the first iteration when C 2 jumps from 0 to a non-zero value. This figure shows, as expected, that the d k and b k attacking     policies generate in a couple of iterations a new cluster composed by a significant amount of nodes (peaking roughly at 40% of the nodes), while random and c k focal methods generate smaller subgraphs and the transition takes in average, a higher attacking effort. Interestingly, the transition to an unconnected graph composed of multiple functioning clusters for the S B case, takes several iterations (at least 15 iterations) and it occurs closer to the total blackout transition (about 30 iterations) than any other methodology. This behavior is further complemented with figure 6(d), here, the S B focal strategy generates very few clusters as compared to the others; this is likely related again to the tree-like topology of the Colombian power grid and to the observed phenomenon in [21], where nodes located on dead-tree arrangements (more specifically, inner tree nodes [22]), were found to exhibit, in general, a lower SNBS than the rest of the nodes. As these are some of the nodes that are more likely to be removed first under S B focal attacks, and as by definition they do not have a significant amount of connections, the graph does not segregate on multiple clusters that easily. Table 2 summarizes the average value of the transitions for the node removal algorithm applied to the Colombian power grid.

Node removal resilience for synthetic power grids
To test the generality of the results shown for the Colombian power grid, the process is repeated for a set of synthetic power grids, which are generated as described in section 2.1. Results on the nodal resilience for these synthetic power grids are presented in figure  7. It is surprising that, besides the wider standard deviation in every trace (which was to be expected considering the randomness in the generation of the complex network), every other behavior related to the elimination based on b k , c k , d k and even the random attacking method is identical to that observed in the specific case of the Colombian power grid in figure 6. This confirms the fact that the algorithm proposed in [39] is indeed capable of reproducing real-world examples of networked power systems, but also reveals that an homogeneous power distribution did not have, in average, a significant impact in the transition to blackout for these specific attacking schemes. This is not surprising, as b k , c k and d k are structure-based vulnerabilities and completely independent of the dynamical vulnerabilities. The S B -focused attacking, however, shows a very different evolution, similar to that of the random attacks, except in figure 7(d), where S B attacks keep producing the lowest amount of clusters, which is related to the same explanation given before: the nodes located at dead-tree arrangements usually possess low SNBS levels, and their removal hardly generates new functional clusters. For reference, the values for the transitions T 1 and T 2 on the synthetic power grids are presented in table 3. It seems quite natural to think that, as the random networks have more generator nodes and they are randomly distributed over the graph, the iteration for which the transition to blackout occurs would, for any of the attacking methods, be higher than the Colombian case since after each removal step the probability to obtain non-operational clusters would be reduced, yet it only happens for the S B -focused one, as can be also confirmed by comparing tables 2 and 3. Next section focuses the discussion around the effects of these power heterogeneity in the resilience profiles observed.   Figure 8 shows the amount of generators N g and consumers N c after each iteration in both, the Colombian power grid and the synthetic power grids, and it allows to see how the node elimination is working: for the random synthetic power grids in average, both, consumers and generators are being eliminated at the same rate, this leads to believe that every new cluster created in the network preserves a comparable number of both kind of nodes and that homogeneity makes the total destruction of said cluster less likely. This also explains why, in average, T 2 for the S B -method is close to that of the random case (table 3) and hints that the main reason why clusters are being removed from the graph is because the synchronization condition is not satisfied (step 5 in the resilience testing algorithm), rather than the absence of generators and consumers (step 3 of the algorithm). Nevertheless, for the Colombian power grid the decreasing slope for N g is faster than that of N c for the first iterations. Since the initial amount of generators is lower than the initial amount of consumers in this test case, after some of the nodes in the former group have been removed from the graph, the resulting clusters are more likely to be consumers isolated without a connected generator, subsequently, these clusters go to blackout. This implies that the main reason why clusters are being   suppressed from the Colombian grid is the absence of generators in them (step 3 of the algorithm), and explains why T 2 in the S B -focused method is the lowest value in table 2.

Influence of power heterogeneity
In addition, as mentioned in the algorithm, the power supplied by each generator is modified after each percolation step in order to maintain the power balance in the grid. Figure 9 shows the SNBS of each generator node in the complex network as a function of P . For both study cases: the Colombian power grid and the randomly grown synthetic power grids, it can be observed that, as the power of the node increases, its stability diminishes; this comes from the fact that P acts as a measure of the stress in the node, that is, the energy demand that it needs to satisfy: having to provide energy for a higher amount of consumer nodes, a generator becomes then more susceptible to random perturbations, as indicated by the reduction on S (i) B . After performing a Pearson correlation testing between S (i) B and P (i) , a correlation coefficient σ and a p-value p val was found as (σ, p val ) = (−0.605, 0) for the Colombian grid and (σ, p val ) = (−0.221, 0) for the synthetic grids, confirming the existence of a negative relationship between both variables; which is in accordance with the general behavior observed in the recent results shown in [28]. This result, as well as the previous observation on the resilience to blackout for S B , can be contrasted with the results presented in [16], where it was found that, by distributing the demand among multiple small power generators, instead of a few big power plants, the value of the critical coupling K c usually diminishes, thus synchronization is favoured, but for disturbances applied to power consumption, the most robust performance was found when there is a mixture of both, small and big power generators.

Distribution of vulnerable nodes
Let us now consider weak nodes as those that exhibit a SNBS lower than 0.4, then the fraction of weak nodes can be expressed as: where N w(I) and N (I) are respectively the number of weak nodes and the total amount of nodes in the graph at iteration I of the resilience testing algorithm. Figure 10 shows the behaviour of η w after each iteration for the S B -focused attacking algorithm. From this figure, a rather counterintuitive observation emerges: in general, for both test cases, removing the most vulnerable nodes in the network is causing an overall reduction of the stability of the whole power grid. This implies that, although attacking nodes with a poor SNBS does not ensure a faster transition to blackout than other strategies, it does make the whole grid more prone to dynamical disturbances in oscillation frequency and phase angle. It is worth noting that, results on figure 10 are statistically less reliable in the last iterations due to the fact that each power grid reaches total blackout at a different iteration value in general, thus, last iterations are not averaging necessarily over 10 realizations as the first iterations are. Another important analysis that can be performed is observing how the critical iterations for which the phase transitions occur, are affected for the parameters in the system. In that regard, for the node resilience procedure over synthetic power grids, T 1 and T 2 were calculated over variations in the coupling strength K as shown in figure 11, where each data point is the average of 10 simulations for the SNBS method and 500 simulations for the other methods. It should be noted that by decreasing K below 10, T 1 for the SNBS case becomes lower than the random case in average, hinting that the effect where weaker nodes are located mainly in inner-tree nodes is more noticeable when K is somewhat large, while for lower K, the weak nodes are spread all around the graph. This is understandable since, in general, it can be expected that the basin stability of a node is increased when the coupling of the network is enhanced (this behavior however is not monotonic for all nodes, as pointed out in [44,45]). For the random and the topology-dependant elimination methods, no tendency is observed for T 1 when varying K.
T 2 on the other hand shows an evident increase for any removal method when the coupling strength is increased. As shown in figure 11(b), random and topologicalbased elimination methods seem to reach a plateau when K is sufficiently high. This reduction in the variability of T 2 is because for a higher K, the network becomes more robust to lose clusters due to lack of synchronization (step 5 in the algorithm) and thus, clusters are removed merely by the lack of generators and consumers (step 3) or by single elimination (step 1 and 2). The number of attacks required to reach blackout on either of those two cases depends mostly on the initial construction of the network, therefore it should not vary significantly among multiple experiments. The SNBS removal depends on the system dynamics, thus its curve does not flat as the others for the observed range, however it is expected to behave in the same way for a significantly higher K, when the set point in any node becomes globally stable in Π.

Edge removal resilience
Continuing now with the analysis for the edge removal procedure, figure 12 presents the evolution of the Colombian graph under this methodology. The difference between the random and focal schemes is remarkable on these simulations: by focusing attacks on transmission lines with higher loads ∆ θ (i, j), the size of the power grid goes down rapidly and reaches total blackout faster than other alternatives.
By focusing on lines with a higher e k , the sparsity of the graph grows rapidly and many small isolated but functional clusters form, as appreciated in figure 12(d), where the average number of clusters for e k can even be twice the amount seen in the random case and almost three times the ∆ θ (i, j) case. Qualitatively, exactly the same behavior is observed on synthetic power grids on figure 13, and actually, the Colombian grid results lay inside the standard deviation of the experiments for synthetic grids. The values found for T 1 and T 2 on these simulations are summarized on table 4. Finally, for the line resilience procedure over synthetic power grids, the behavior of the transitions is observed on figure 14. For T 1 on panel (a), no clear tendency   is observed, besides the expected fact that this transition is, in general, higher for the random method than for other elimination strategies. Panel (b) however shows something interesting for T 2 : for any value of K, attacking the most vulnerable edges based on the dynamics (∆ θ ) produces the fastest transition to blackout in the synthetic power grids, while attacking randomly produces in general the slowest transition. Furthermore, this transition also reaches the plateau when K is sufficiently large, for any elimination method, as observed for the node attacking results. Evidently, the dynamical-based elimination with ∆ θ flats so quickly (K = 6) because increasing K strengthens the transmission lines and thus directly raises ∆ θ , as opposed to the node elimination based on SNBS, where the basin stability is known to be affected by K but not in the same trivial and immediate way [44,45].

Concluding remarks
In this paper, an algorithm was proposed to test the resilience of power grids against failures. Our approach differs from others found in the literature since it takes into consideration dynamical stability against random finite-size disturbances to describe nodal vulnerability, as well as linear stability of phase differences to compute edge vulnerability. Sequential elimination of elements in the network based on these dynamical vulnerabilities was also compared with the same procedure applied to structural vulnerabilities, which were defined in terms of connectivity patterns in the topology of the graph. This allowed to extract some interesting findings which will be discussed in the following. First of all, focusing nodal attacks on redundantly connected nodes, such as those with the highest clustering coefficient, is obviously the least effective strategy, being  comparable with the random removal; the blackout transition T 2 for both is always the highest, meaning that they require more attacks in order to destroy the whole network. Similarly, the giant component of the graph reduces at the slowest rate, since dividing the graph into multiple operational clusters is not likely until the redundant links are removed from the network. Eventually, both methods become essentially the same, once the clustering coefficient for the whole grid is uniform. The fastest way to reduce the size of the giant component of the network is, as expected, to attack the most central nodes, either those with the highest degree or betweenness; when these bulk nodes are removed, the graph rapidly segregates into multiple perfectly operational clusters. In other words, the transition to an unconnected graph T 1 is always the lowest for these methods.
For the case of node elimination based on basin stability it was found that the   transition to blackout T 2 in the randomly grown synthetic power grids did not differ significantly to that found for the random removal approach. However, in the specific case of the Colombian power grid, this method proved to be the most effective to bring the whole network to total blackout. Essentially, the only key difference between both test cases is the power distribution: for the synthetic grids the power was initially distributed uniformly, with an equal number of generators and consumers, but for the Colombian case only 27% of the initial nodes are generators, and the rest consumers. This power homogeneity in the synthetic grids then is enhancing the resilience of the graph, and this was confirmed when it was observed that for the Colombian grid, the first nodes to be removed are generators, leading to a higher stress in the remaining power plants: a few generators will have to increase drastically the power they produce in order to balance the demand. In addition, this power increase in the generators will reduce their basin stability, that is, it will make them more prone to dynamical random disturbances in oscillatory frequency and phase angle. This negative correlation between the power parameter of a node and its non-linear stability is supported by the findings in [28]. Interestingly, the number of weak nodes in the power grid (nodes with a low basin stability level), was found to increase during the node removal procedure focused on basin stability for both test cases. This implies that this attacking strategy is useful to undermine the general dynamical stability of the power grid and make it more prone to random disturbances, even though, it does not necessarily propagate structural failures to larger proportions of the network (it depends on the amount of generators, as already explained).
Regarding the line resilience testing performed, it was found that for either the Colombian power grid or the synthetic power grids, attacking highly central edges, as those with a high edge betweenness centrality, the graph segregates rapidly into multiple clusters and the giant component reduces abruptly. The transition to blackout, however, is not faster than random elimination. It is, as expected, a very similar behaviour to that found in node elimination through node betweenness. The elimination based on the dynamical measure of phase difference (∆ θ ) was found to be the most efficient in order to bring the whole network to blackout, and remarkably, this behaviour can be observed for any value tested of the maximum transfer capacity in transmission lines.
Some work is still left open to research such as further exploring the hidden causes that are yielding the reduction of non-linear stability due to heterogeneous power distribution. Similarly, a more detailed study on the phase transitions T 1 and T 2 , as well as the influence of structural parameters in these critical points, could potentially point in the right direction to design the most appropriate connectivity in a graph and parameter distribution such that resilience against cascading failures is optimized. More dynamical measures of vulnerability could also be incorporated in the removal algorithm; for instance, the diffusivity of perturbations [10] that could indicate regions of faster failure propagation in the grid, or the multi-node basin stability [23] computed over the nearest-neighbours, as it would include information about the most vulnerable clusters, rather than individual elements. This would however impose a higher computational effort to a task that is already computationally expensive.

Appendix A. Random growth model for power grids
This section describes the algorithm proposed in [39] which allows to generate complex networks that incorporate a spatial embedding (accounts for geographic position of nodes) and that have similar structural properties to those found in real-world power grids. For this algorithm, a redundancy cost function that will be subject to optimization has to be defined as: where d G (i, j) is the length of the shortest path between nodes i and j in the graph G, d S (y i , y j ) is the spatial distance (Euclidean distance) between the nodes positions y i and y j , and u is a parameter that controls the cost-vs-redundancy trade-off. After defining the parameters from table 1, the construction algorithm is developed through two phases: initialization and growth, which are performed as follows: Appendix A.1. Initialization (i) Initialize a minimum spanning tree for the initial N 0 nodes such that it optimizes the edges that have to be built between them based on the minimum Euclidean distance d S (y i , y j ). -Add a link between the new node h and the closest node j geographically, that is, the node j for which d S (y h , y j ) is minimal. -Generate a random number ζ 2 , if ζ 2 < p add a link between the new node h and some node l for which f (h,l,G) is maximal. -Generate a random number ζ 3 , if ζ 3 < q draw a node h from the network uniformly at random. Then find the node l that is not yet linked to h and for whichf (h,l,G) is maximal and connect them. (c) Otherwise (if ζ 1 > 1 − s), select an edge (i, j) of the network uniformly at random. The new geographic location of node h will be given by y h = (y i +y j ) /2, so the edge (i, j) is removed from the graph and then the edges (i, h) and (h, j) are added.