Decreased resilience in power grids under dynamically induced vulnerabilities
Cristian Camilo Galindo-González, David Angulo-García, Gustavo Osorio
DDecreased resilience in power grids underdynamically induced vulnerabilities
C. C. Galindo-Gonz´alez , D. Angulo-Garcia and G. Osorio Grupo de Percepci´on y Control Inteligente (PCI). Departamento de Ingenier´ıaEl´ectrica, Electr´onica y Computaci´on. Universidad Nacional de Colombia - SedeManizales. Kil´ometro 9, v´ıa al aeropuerto La Nubia. Manizales, Colombia. Grupo de Modelado Computacional - Din´amica y Complejidad de Sistemas.Instituto de Matem´aticas Aplicadas. Universidad de Cartagena. Carrera 6 [email protected], [email protected],[email protected]
July 2020
Abstract.
In this paper, a methodology inspired on bond and site percolationmethods is applied to the estimation of the resilience against failures in power grids.Our approach includes vulnerability measures with both dynamical and structuralfoundations as an attempt to find more insights about the relationships betweentopology and dynamics in the second-order Kuramoto model on complex networks.As test cases for numerical simulations, we use the real-world topology of theColombian power transmission system, as well as randomly generated networks withspatial embedding. It is observed that, by focusing the attacks on those dynamicalvulnerabilities, the power grid becomes, in general, more prone to reach a state oftotal blackout, which in the case of node removal procedures it is conditioned by thehomogeneity of power distribution in the network.
Keywords : Kuramoto model, power grid, basin stability, resilience, stability, percolation.
1. Introduction
Emergence of synchronization in systems of agents, coupled through local interactions,is a widely studied phenomenon with multiple applications in physics, biology and socialsciences [1, 2, 3]. In particular, the optimal functioning of a power grid, which representsone of the most complex interacting systems in engineering, highly depends on its abilityto maintain a synchronous operation over time despite external disturbances. Losingthat synchrony, even locally, may lead to cascading failures and complete blackout of thenetwork [4, 5, 6]. Two main concerns have motivated the discussion around dynamicalanalysis and design of power grids in the last years, those are: the strong economicaland social impact that a power outage could cause in the highly electricity-dependant a r X i v : . [ n li n . AO ] J u l ecreased resilience in power grids under dynamically induced vulnerabilities ecreased resilience in power grids under dynamically induced vulnerabilities
2. Model and methods
In this paper, the framework of the synchronous motor representation is used to modelthe general structure and dynamics of power grids [11, 12], since it has been employedconstantly and successfully in the past years, as mentioned in the previous section. Thisapproach is described in the following.
As a simplified model of a real-world power grid, consider a system of N interactingsynchronous machines arranged in a connected graph G ( ϑ, ε ), with a set of nodes ϑ = { , , ..., 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) orconsumer machines (demand energy from the grid), thus ϑ = ϑ g ∪ ϑ c , being ϑ g the setof generators and ϑ c the set of consumers.Let the dynamical state of each machine be represented by its phase angle φ i andits phase velocity ˙ φ i := dφ i / dt , i ∈ ϑ . Under an appropriate operation of the powergrid, 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 , thenthe dynamical behavior of θ i and its velocity ˙ θ i can be described by the well-known ecreased resilience in power grids under dynamically induced vulnerabilities (a) Dense
SproutSparse
SproutProper
LeafRootBulkInner
Tree d k N o d e s (b) c k N o d e s (c) b k × N o d e s (d) e k × E d g e s (e) t θ (f) t ˙ θ (g) Figure 1: Colombian Power Grid: (a) Topology of the network with N = 102 and M = 158. Generating nodes N g = 28 (cross symbols) and N c = 74 consumer nodes(circle symbols). (b)-(e) Distribution of the node degree d k , clustering coefficient c k ,node betweenness centrality b k and edge betweenness centrality e k . (f)-(g) Time trace ofthe phases and phase velocities of each oscillator when the grid converges to a frequencysynchronized state (red lines: generators, blue lines: consumers).second-order Kuramoto model for coupled oscillators given by [12, 16]:¨ θ i ( t ) = P i − α i ˙ θ i ( t ) + K N (cid:88) j a ij sin( θ j ( t ) − θ i ( t ) ) (1)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 dampingcoefficient 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. Theconstant K is the maximum power transfer capacity between any pair of connectednodes, which depends on the impedance of transmission lines connecting them. Thisequation is derived when the ohmic losses on transmission lines are neglected, an equalimpedance level is assumed for every transmission line in the network and the deviationsfrom the phase velocity at any node as compared to the reference Ω are small, that is | ˙ θ i | (cid:28) Ω, ∀ i ∈ ϑ [12, 16].The dynamical system described by equation (1) is said to be in a phasesynchronized state if θ i ( t ) = θ j ( t ) , ∀ i, j ∈ ϑ , and in a frequency synchronized stateif ˙ θ i ( t ) = ˙ θ j ( t ) = ω s , ∀ i, j ∈ ϑ . Since we are interested in the dynamics around thesynchronization 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 NationalTransmission System [36], a power grid composed of 158 edges and 102 nodes (28 ofthem are generators and 74 are consumers). Power supplied by generators is set to ecreased resilience in power grids under dynamically induced vulnerabilities P k = 1 . ∀ k ∈ ϑ g , and the power drained by consumers is uniformly distributed suchthat the power balance condition in the network is fulfilled, that is: N (cid:88) i =1 P i = 0 . (2)This condition is a requirement in order to allow the existence of a locally-stableequilibrium point in the system [16]. Finally, the damping and maximum power transferare fixed to α i = 0 . K = 12 . θ ( t ) and ˙ θ ( t ) for every oscillator when the system convergesfrom an unsynchronized state to a frequency synchronized regime given the parametersabove. Aiming for a generalization of the results, random test cases were also usedand compared to the Colombian test case just described; to do that, synthetic powergrids were generated with the algorithm proposed in [39], which is specially useful asit creates spatially embedded networks that emulate connectivity distributions of real-world power grids. Synthetic networks were constructed by initializing a tree-like basearchitecture and then following a growth procedure conditioned by the optimizationof a redundancy-cost function that depends on both, geographic distance and numberof paths separating each pair of nodes in the graph. For a detailed description of thealgorithm, refer to Appendix A. The parameters used for the growth model algorithmwere chosen to match those used on [22], as presented in table 1. Note that the finalsize of the network is fixed to N = 102 so that it matches the Colombian power gridsize, nevertheless, geographic positions for nodes are chosen uniformly at random andan equal amount of generators and consumers is placed to reduce heterogeneity in thesamples (thus, P k = 1 . ∀ k ∈ ϑ g and P j = − . ∀ j ∈ ϑ c ).Table 1: Parameters for the construction of synthetic power grids. Parameter Value N : An initial amount of nodes, N ≥
1. 1 N : Final amount of nodes that the network will have N > N . 102 y = { y , y , ..., y N } : Geographical location of nodes. y k ∼ U [0 , p : Probability of constructing additional redundancy linksattached to new nodes (0 ≤ p ≤ / q : Probability of constructing additional redundancy linksbetween existing nodes (0 ≤ q ≤ / s : Probability of splitting an existing line (0 ≤ s ≤ / u : Cost-vs-redundancy trade-off parameter for equation A.1. / ecreased resilience in power grids under dynamically induced vulnerabilities A helpful indicator that quantifies the degree of synchronization is the complex-valuedorder parameter r defined as [40]: re i Ψ ( t ) = 1 N N (cid:88) j e iθ j ( t ) (3)Here, Ψ ( t ) is the average phase of the oscillators at time t . If all phases θ j areidentical, the magnitude of the order parameter is r ( t ) = 1. Conversely, if they are equallydistributed 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, inthe case of power grids in which equation (2) holds, when the system is synchronizedthe average phase is Ψ ( t ) = 0, and the real part of the order parameter IR[ r ] = 1, whileit oscillates around zero otherwise. While the dependence of the order parameter withtime is useful in transient analysis, we will focus on the steady state behavior, so wecan define the average order parameter in steady state as [16]: r ∞ = lim t →∞ (cid:20) lim t →∞ (cid:18) t (cid:90) t + t t r ( t ) dt (cid:19)(cid:21) (4)Another measure of frequency synchronization in power grids is the squared averagerotational speed v t ) defined as [16]: v t ) = 1 N N (cid:88) j ˙ θ j ( t ) (5)and the associated steady state value of the speed v ∞ : v ∞ = (cid:115) lim t →∞ (cid:20) lim t →∞ (cid:18) t (cid:90) t + t t v t ) dt (cid:19)(cid:21) . (6) Following the results in [4], it can be proven that, by applying a small-angleapproximation, the phases at the frequency synchronized equilibrium point, which aredenoted by θ ∗ ∈ IR N , can be estimated as: θ ∗ ≈ K L † P (7)where L † is the pseudo-inverse of the Laplacian matrix, which is computed throughthe adjacency matrix as L = diag ( (cid:80) nj =1 a ij ) − A . Let us now define the oriented incidencematrix B = { b ij } ∈ R N × M , with: b ij = i incides in edge j − j incides in node i ecreased resilience in power grids under dynamically induced vulnerabilities G is an undirected graph, a direction for each edge can beassumed without loss of generality for the following analysis, allowing for this definitionof B . Under these circumstances, the steady state phase difference between any pair ofconnected nodes ( i, j ) ∈ ε can be approximated by:∆ θ ≈ K B T L † P (9)with ∆ θ ∈ R M . Then, the system (1) is guaranteed to have a unique and stable solutionwith synchronized frequencies and cohesive phases if the next sufficient condition for ∆ θ holds: || ∆ θ || ∞ < sin( γ ) (10)with 0 ≤ γ < π/
2. In the limit γ → π/
2, equation (10) simply reduces to || ∆ θ || ∞ < π / [4]. This rather simple equation, readily connects topological and dynamical featuresof the network. While condition (10) is only sufficient, it is a rapid way to assesssynchronization without relying on the computationally expensive time evolution of alarge number of phase oscillators. Moreover, equation (10) provides an approximationto the critical value K c of the coupling strength when substituting from (9): K c ≈ || B T L † P || ∞ (11)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 realpart 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 ≈ .
51, thesystem achieves a state of phase cohesiveness characterized by IR[ r ∞ ] ≈ v ( t ) abruptly decreases to zero, meaningthat 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 goodagreement with the numerical value of the transition. This evidence supports the useof condition (10) as a synchronization requisite in the following analysis. Note that thechosen value of K = 12 . K > K c , thus it allowsfor a synchronized state to exist.As a measure of the vulnerability of each transmission line, we take advantage ofthe equation (9), given that the edge ( i, j ) is considered weak when | ∆ θ ( i, j ) | is close to1, while it is considered more stable when it is close to 0, regarding the synchronizationlimit imposed by (10). This approach is equivalent to that presented in [41, 42], where ecreased resilience in power grids under dynamically induced vulnerabilities K I R [ r ∞ ] K c ≈ . a ) t IR[ r ( t ) ] t IR[ r ( t ) ] K v ∞ K c ≈ . b ) t v ( t ) t v ( t ) Figure 2: Synchronization for the Colombian power grid: (a) Average steady statevalue of the real part of the order parameter as a function of the coupling strength K . Left (right) inset shows the time trace of IR[ r ( t ) ] before (after) critical coupling. (b)Average rotational speed as a function of K . As in (a), insets show time traces beforeand after critical coupling. In both cases, vertical dashed line depicts the value of K c calculated with equation (11).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 immediatelyfrom equation (9). 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 ∗ , notedby β , 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 ofthe system (1). Similarly, the likelihood of a randomly perturbed trajectory to returnback to β , known as the basin stability S B , can be defined as [20, 21]: S B ( β ) = (cid:90) Γ β ( x ) ρ ( x ) dx (12)where Γ β ( x ) is a function that indicates whether a state x belongs to the basin ofattraction of X ∗ or not, that isΓ β ( x ) = (cid:40) ∀ x ∈ β ∀ x / ∈ β (13)The function ρ ( x ) is the density of states to which the system can be pushed bysome non-local perturbation, such that (cid:82) 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 synchronousstate is unstable and S B = 1 when it is globally stable. ecreased resilience in power grids under dynamically induced vulnerabilities S B by randomlysampling a sufficiently high number of disturbed states I C (following a certaindistribution ρ ( x ) ) from a representative subspace Π ⊂ X , and using them as the initialconditions for time-evolution simulations. For this work, ρ ( x ) is chosen as an uniformdistribution in the restricted subspace Π, that is: ρ ( x ) = | Π | ∀ x ∈ Π0 ∀ x / ∈ Π (14)Finally, the amount F C of simulated trajectories that are found to approachasymptotically to the attractor is assumed to be proportional to the volume of the basinof 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 therobustness of a node against large perturbations applied to it, by means of the single-node basin stability (SNBS) [21, 23], given by: S ( i ) B = F C I C , i ∈ ϑ (15)where the I C initial conditions are drawn from perturbations applied to node i . Asillustration, 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 totaldatapoints. Throughout this work, the disturbances for each node are drawn from thesubspace Π = [ − π, π ] × [ − , Equations (9) and (15) define already a dynamical vulnerability for edges and nodes,respectively. Furthermore a vulnerability based merely on the connectivity features ofthe complex network could also be considered, since it is intuitive that by removinga highly connected node or a certain critically located edge could lead to a rapiddiminishing on power grid stability and performance. In that regard, let us definestructural vulnerability in terms of the following centrality measures: • Degree centrality d ( i ) k : Amount of nodes to which node i is connected, normalizedby the maximum possible degree in the network. It can be computed as: d ( i ) k = 1 N − N (cid:88) j a ij (16) • 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: c ( i ) k = T i d ( i ) k ( d ( i ) k −
1) (17) ecreased resilience in power grids under dynamically induced vulnerabilities i of the Colombian grid when random disturbancesare applied to θ i or ˙ θ i . Green crosses are initial conditions that return to the frequencysynchronized state, while the blue circles are out of the basin of attraction. The reddashed line denotes the subspace Π from which disturbances are chosen. This diagramwas created by sampling I C = 500 points. • Node betweenness centrality b ( i ) k : Sum of the shortest paths between every pair ofnodes ( s, t ) in the network that pass through node i : b ( i ) k = (cid:88) s (cid:54) = t (cid:54) = i ∈ ϑ σ s,t ( i ) σ s,t (18)where σ s,t is the amount of shortest paths between nodes s and t and σ s,t ( i ) is thenumber of those paths that pass through node i [38]. • Edge betweenness centrality e ( l ) k : In the same fashion as b k measures centrality fora node, e k does it for an edge; it is simply defined as the number of shortest pathsin the network that include edge l : e ( l ) k = (cid:88) s (cid:54) = t ∈ ϑ σ s,t ( l ) σ s,t , l ∈ ε (19)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. ecreased resilience in power grids under dynamically induced vulnerabilities (a) S B S B -1 N o d e s (b) S B S B -1 N o d e s (c) S B S B -1 N o d e s (d) S B S B -1 N o d e s (e) S B S B -1 N o d e s (f) S B S B -1 N o d e s (g) S B S B -1 N o d e s (h) S B S B -1 N o d e s Figure 4: Node removal algorithm under the S B -focused attacking scheme. Circle nodesindicate consumers while crosses represent generators. Node color is mapped with the S ( i ) B (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 showsthe histogram of the basin stability for the whole network. One realization is presentedin panels (a)-(d) and a different one in panels (e)-(h) in order to visualize two differentmechanisms of cluster vanishing.
3. Resilience measures
In order to assess the resilience of a power grid, we propose a percolation-basedalgorithm. To do so, we will follow the capability of the network to maintain itsfunctioning upon different node removal (site percolation) and edge removal (bondpercolation) rules. In particular, the algorithm that will be described below removes anode or an edge from the graph on each iteration by applying a random attack, wherethe element to be removed is chosen uniformly at random; or a focal attack, where themost vulnerable node or edge in the graph is removed first. Here, the most vulnerablenode 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 vulnerableedge 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 focalattack, 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. ecreased resilience in power grids under dynamically induced vulnerabilities K > K c , where K c is computedby the approximation (11). If it does not, the whole cluster is removed from thegraph.(vi) Measure the needed observables and return to step 1 until the whole power gridhas gone to blackout, that is, no operational cluster exists in the network.An important remark has to be done about the algorithm, regarding node removalunder the focal attack methodology based on SNBS. Computing this vulnerability indexis computationally expensive since it requires a large amount of time-domain simulationsin 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, whilefor each of the other cases, 500 simulations are performed. It is also worth noting thatthe fine structure of the basin of attraction is not captured by this approach, and itcould impose numerical problems on dynamical systems that posses either fractal basinboundaries or riddled or intermingled basins [43].For this elimination procedure, two main phase transitions will be analysed: • Transition T : This transition occurs at the iteration when the graph is dividedinto multiple functional subgraphs for the first time (turning from a connectedgraph to an unconnected graph). • Transition T : The final iteration of the algorithm, when the power grid goes tocomplete 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 clustersafter removing one node. Note that in the frame (c) there is only one generator inthe northern cluster, and interestingly, it has a very poor basin stability level, so whenthat node is removed from (c) to (d), the whole northern cluster goes to blackout. Adifferent realization is presented in panels (e)-(h); note that from (g) to (h) it suffices toremove one node to bring the bigger cluster to blackout even though there would stillremain connected subgraphs with both generators and consumers. This occurs becausethe remaining cluster would have a topology that does not sattisfy the synchronizationcondition (10) and thus it vanishes. Similarly, figure 5 shows the algorithm in actionunder ∆ θ -focused attacks. In the transition from panel (b) to (c), it can be appreciatedthat, after suppressing a very vulnerable transmission line, a significant cluster ofisolated consumers in the east side of the network vanishes. ecreased resilience in power grids under dynamically induced vulnerabilities (a) ∆ θ ∆ θ -1 E d g e s (b) ∆ θ ∆ θ -1 E d g e s (c) ∆ θ ∆ θ -1 E d g e s Figure 5: Edge removal algorithm under the focal attacking scheme. Edge color ismapped with the ∆ θ of the nodes it connects. The red arrow indicates the edge that isgoing to be supressed. Also, the inset shows the histogram of the phase differences forthe whole network. Figure 6(a) shows the fraction of nodes in the largest cluster of the network (also calledthe giant component of the graph and denoted by C ) subject to node removal with focaland random attacking policies. Similarly, figure 6(b) presents the fraction of removednodes after each iteration with respect to the original size (named removed componentand denoted by ( N (0) − N ( I ) ) / N (0) , being N (0) and N ( I ) the original size of G and its sizeafter I removal iterations, respectively), which allows to detect the phase transition T ,corresponding to the iteration at which ( N (0) − N ( I ) ) / N (0) = 1. It is thus observed that thegrid goes to complete blackout around 30th removal iteration for the case of the focalattacking 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 havethe 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 fastertransition to the total blackout state of the network than just removing nodes atrandom implies that there is indeed a relationship between the stability of a nodeagainst dynamical disturbances (related to the SNBS) and the structural connectivityand resilience to cascading failures of the graph.Furthermore, by focusing the attacks on nodes with the highest b k or d k , thedimension of the giant component goes down drastically during the first eliminationiterations and then it goes down at a slower rate, reaching total blackout in anintermediate level of attacks between the S B and random methodologies. c k on theother hand, produces a slower reduction of the graph size and a total blackout with themost iterations in average. This is related to the tree-like structure of the power grid thatcan be seen on figure 1. In this kind of topologies, the nodes with the higher degree andbetweenness are located at the bulk regions of the graph, therefore, their removal likelydivides the giant component instantly into multiple smaller clusters, as also confirmed bythe figure 6(d), where it is shown that these methodologies produce by far, the highest ecreased resilience in power grids under dynamically induced vulnerabilities Iteration C RandomFocused on d k Focused on c k Focused on b k Focused on S B (a) Giant component of the network. Iteration N ( ) − N ( I ) N ( ) RandomFocused on d k Focused on c k Focused on b k Focused on S B (b) Removed component of the network. Iteration C RandomFocused on d k Focused on c k Focused on b k Focused on S B (c) Second largest component of the network. Iteration | S G | RandomFocused on d k Focused on c k Focused on b k Focused on S B (d) Amount of operational clusters in the network. Figure 6: Evolution in the structure of the
Colombian power grid subject to the nodalelimination algorithm. The shading accounts for the standard deviation calculated over10 realizations with I C = 50 for the S B -focused methodology and 500 realizations forthe other attacking schemes.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 theneighbours that it is connected to, are also connected between them, hence, removingit is not going to divide the graph and generate new clusters. Eventually, after multiplenode 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 attackingstrategy behaves roughly in the same way as the random attacking one after a fewiterations.Figure 6(c) plots the size of the second largest component of the graph C , whichreveals the location of the phase transition T as the first iteration when C jumpsfrom 0 to a non-zero value. This figure shows, as expected, that the d k and b k attacking ecreased resilience in power grids under dynamically induced vulnerabilities Iteration C Focused on S B RandomFocused on d k Focused on c k Focused on b k (a) Giant component of the network. Iteration N ( ) − N ( I ) N ( ) Focused on S B RandomFocused on d k Focused on c k Focused on b k (b) Removed component of the network. Iteration C Focused on S B RandomFocused on d k Focused on c k Focused on b k (c) Second largest component of the network. Iteration | S G | Focused on S B RandomFocused on d k Focused on c k Focused on b k (d) Amount of operational clusters in the network. Figure 7: Evolution in the structure of the randomly grown synthetic power grids subject to the nodal elimination algorithm. The shading accounts for the standarddeviation calculated over 10 realizations with I C = 50 for the S B -focused methodologyand 500 realizations for the other attacking schemes.policies generate in a couple of iterations a new cluster composed by a significant amountof nodes (peaking roughly at 40% of the nodes), while random and c k focal methodsgenerate smaller subgraphs and the transition takes in average, a higher attacking effort.Interestingly, the transition to an unconnected graph composed of multiple functioningclusters for the S B case, takes several iterations (at least 15 iterations) and it occurscloser 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 strategygenerates very few clusters as compared to the others; this is likely related again tothe 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 ecreased resilience in power grids under dynamically induced vulnerabilities S B focalattacks, and as by definition they do not have a significant amount of connections, thegraph does not segregate on multiple clusters that easily.Table 2 summarizes the average value of the transitions for the node removalalgorithm applied to the Colombian power grid.Table 2: Average value of T and T for the nodal resilience testing algorithm in theColombian grid. Random Node b k c k d k S kB T T To test the generality of the results shown for the Colombian power grid, the process isrepeated for a set of synthetic power grids, which are generated as described in section2.1. Results on the nodal resilience for these synthetic power grids are presented in figure7. It is surprising that, besides the wider standard deviation in every trace (which wasto 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 randomattacking method is identical to that observed in the specific case of the Colombian powergrid in figure 6. This confirms the fact that the algorithm proposed in [39] is indeedcapable of reproducing real-world examples of networked power systems, but also revealsthat an homogeneous power distribution did not have, in average, a significant impactin 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 thedynamical vulnerabilities. The S B -focused attacking, however, shows a very differentevolution, similar to that of the random attacks, except in figure 7(d), where S B attackskeep producing the lowest amount of clusters, which is related to the same explanationgiven before: the nodes located at dead-tree arrangements usually possess low SNBSlevels, and their removal hardly generates new functional clusters. For reference, thevalues for the transitions T and T on the synthetic power grids are presented in table 3.Table 3: Average value of T and T for the resilience testing algorithm in the syntheticpower grids. Random Node b k c k d k S kB T T ecreased resilience in power grids under dynamically induced vulnerabilities S B -focused one, as can be alsoconfirmed by comparing tables 2 and 3. Next section focuses the discussion around theeffects of these power heterogeneity in the resilience profiles observed. Iteration N N g N c (a) Colombian power grid. Iteration N N g N c (b) Synthetic power grids. Figure 8: Number of consumers N c and generators N g existing in the graph on eachiteration of the S B -focused removal algorithm. The shading accounts for the standarddeviation over the 10 random realizations with I C = 50.Figure 8 shows the amount of generators N g and consumers N c after each iterationin both, the Colombian power grid and the synthetic power grids, and it allows to seehow 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 tobelieve that every new cluster created in the network preserves a comparable numberof both kind of nodes and that homogeneity makes the total destruction of said clusterless likely. This also explains why, in average, T for the S B -method is close to thatof the random case (table 3) and hints that the main reason why clusters are beingremoved from the graph is because the synchronization condition is not satisfied (step 5in the resilience testing algorithm), rather than the absence of generators and consumers(step 3 of the algorithm). Nevertheless, for the Colombian power grid the decreasingslope for N g is faster than that of N c for the first iterations. Since the initial amount ofgenerators is lower than the initial amount of consumers in this test case, after some ofthe nodes in the former group have been removed from the graph, the resulting clustersare 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 ecreased resilience in power grids under dynamically induced vulnerabilities P i S ( i ) B (a) Colombian power grid. P i S ( i ) B (b) Synthetic power grids. Figure 9: SNBS of each generator node as a function of the power it supplies to thenetwork. Red dashed line marks the initial set power P = 1 . T in the S B -focused method is the lowest value in table 2.In addition, as mentioned in the algorithm, the power supplied by each generatoris modified after each percolation step in order to maintain the power balance in thegrid. Figure 9 shows the SNBS of each generator node in the complex network asa function of P . For both study cases: the Colombian power grid and the randomlygrown 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 inthe node, that is, the energy demand that it needs to satisfy: having to provide energyfor a higher amount of consumer nodes, a generator becomes then more susceptible torandom perturbations, as indicated by the reduction on S ( i ) B . After performing a Pearsoncorrelation testing between S ( i ) B and P ( i ) , a correlation coefficient σ and a p-value p val was found as ( σ, p val ) = ( − . ,
0) for the Colombian grid and ( σ, p val ) = ( − . , S B , can be contrasted with the results presented in [16], where it was foundthat, by distributing the demand among multiple small power generators, instead ofa few big power plants, the value of the critical coupling K c usually diminishes, thussynchronization is favoured, but for disturbances applied to power consumption, themost robust performance was found when there is a mixture of both, small and bigpower generators. ecreased resilience in power grids under dynamically induced vulnerabilities Iteration η w (a) Colombian power grid. Iteration η w (b) Synthetic power grids. Figure 10: Percentage of weak nodes in the network as a function of the percolationiteration.
Let us now consider weak nodes as those that exhibit a SNBS lower than 0.4, then thefraction of weak nodes can be expressed as: η w ( I ) = N w ( I ) N ( I ) (20)where N w ( I ) and N ( I ) are respectively the number of weak nodes and the total amountof nodes in the graph at iteration I of the resilience testing algorithm. Figure 10 showsthe behaviour of η w after each iteration for the S B -focused attacking algorithm. Fromthis 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 thestability of the whole power grid. This implies that, although attacking nodes with apoor SNBS does not ensure a faster transition to blackout than other strategies, it doesmake the whole grid more prone to dynamical disturbances in oscillation frequency andphase angle. It is worth noting that, results on figure 10 are statistically less reliablein the last iterations due to the fact that each power grid reaches total blackout at adifferent iteration value in general, thus, last iterations are not averaging necessarilyover 10 realizations as the first iterations are.Another important analysis that can be performed is observing how the criticaliterations for which the phase transitions occur, are affected for the parameters in thesystem. In that regard, for the node resilience procedure over synthetic power grids, T and T 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 500simulations for the other methods. It should be noted that by decreasing K below 10, T for the SNBS case becomes lower than the random case in average, hinting that theeffect where weaker nodes are located mainly in inner-tree nodes is more noticeable when ecreased resilience in power grids under dynamically induced vulnerabilities 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 anode is increased when the coupling of the network is enhanced (this behavior howeveris not monotonic for all nodes, as pointed out in [44, 45]). For the random and thetopology-dependant elimination methods, no tendency is observed for T when varying K . T on the other hand shows an evident increase for any removal method whenthe coupling strength is increased. As shown in figure 11(b), random and topological-based elimination methods seem to reach a plateau when K is sufficiently high. Thisreduction in the variability of T is because for a higher K , the network becomes morerobust 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 singleelimination (step 1 and 2). The number of attacks required to reach blackout on eitherof those two cases depends mostly on the initial construction of the network, therefore itshould not vary significantly among multiple experiments. The SNBS removal dependson 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 theset point in any node becomes globally stable in Π. Continuing now with the analysis for the edge removal procedure, figure 12 presentsthe evolution of the Colombian graph under this methodology. The difference betweenthe random and focal schemes is remarkable on these simulations: by focusing attackson transmission lines with higher loads ∆ θ ( i, j ), the size of the power grid goes downrapidly and reaches total blackout faster than other alternatives.By focusing on lines with a higher e k , the sparsity of the graph grows rapidly andmany small isolated but functional clusters form, as appreciated in figure 12(d), wherethe average number of clusters for e k can even be twice the amount seen in the randomcase and almost three times the ∆ θ ( i, j ) case. Qualitatively, exactly the same behavioris observed on synthetic power grids on figure 13, and actually, the Colombian gridresults lay inside the standard deviation of the experiments for synthetic grids. Thevalues found for T and T on these simulations are summarized on table 4.Table 4: Average value of T and T for the edge resilience testing algorithm. Colombian grid Synthetic gridsRandom e k ∆ θ Random e k ∆ θ T T T on panel (a), no clear tendency ecreased resilience in power grids under dynamically induced vulnerabilities K T b k c k d k Random S B K s x (a) First phase transition. K T b k c k d k Random S B K s x (b) Second phase transition. Figure 11: Phase transitions in the node percolation process for synthetic power gridsas a function of K . Insets present the standard deviation s x for each computed datapoint.is observed, besides the expected fact that this transition is, in general, higher forthe random method than for other elimination strategies. Panel (b) however showssomething interesting for T : for any value of K , attacking the most vulnerable edgesbased on the dynamics (∆ θ ) produces the fastest transition to blackout in the syntheticpower grids, while attacking randomly produces in general the slowest transition.Furthermore, this transition also reaches the plateau when K is sufficiently large, forany elimination method, as observed for the node attacking results. Evidently, thedynamical-based elimination with ∆ θ flats so quickly ( K = 6) because increasing K strengthens the transmission lines and thus directly raises ∆ θ , as opposed to the nodeelimination based on SNBS, where the basin stability is known to be affected by K butnot in the same trivial and immediate way [44, 45]. ecreased resilience in power grids under dynamically induced vulnerabilities Iteration C RandomFocused on e k Focused on ∆ θ (a) Giant component of the network. Iteration N ( ) − N ( I ) N ( ) RandomFocused on e k Focused on ∆ θ (b) Removed component of the network. Iteration C RandomFocused on e k Focused on ∆ θ (c) Second largest component of the network. Iteration | S G | RandomFocused on e k Focused on ∆ θ (d) Amount of operational clusters in the network. Figure 12: Evolution in the structure of the
Colombian power grid subject to theedge elimination algorithm. The shading accounts for the standard deviation calculatedover 500 realizations.
4. Concluding remarks
In this paper, an algorithm was proposed to test the resilience of power grids againstfailures. Our approach differs from others found in the literature since it takes intoconsideration dynamical stability against random finite-size disturbances to describenodal vulnerability, as well as linear stability of phase differences to compute edgevulnerability. Sequential elimination of elements in the network based on thesedynamical vulnerabilities was also compared with the same procedure applied tostructural vulnerabilities, which were defined in terms of connectivity patterns in thetopology of the graph. This allowed to extract some interesting findings which will bediscussed in the following.First of all, focusing nodal attacks on redundantly connected nodes, such as thosewith the highest clustering coefficient, is obviously the least effective strategy, being ecreased resilience in power grids under dynamically induced vulnerabilities Iteration C RandomFocused on e k Focused on ∆ θ (a) Giant component of the network. Iteration N ( ) − N ( I ) N ( ) RandomFocused on e k Focused on ∆ θ (b) Removed component of the network. Iteration C RandomFocused on e k Focused on ∆ θ (c) Second largest component of the network. Iteration | S G | RandomFocused on e k Focused on ∆ θ (d) Amount of operational clusters in the network. Figure 13: Evolution in the structure of the randomly grown synthetic powergrids subject to the edge elimination algorithm. The shading accounts for the standarddeviation calculated over 500 realizations.comparable with the random removal; the blackout transition T for both is always thehighest, 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 dividingthe graph into multiple operational clusters is not likely until the redundant links areremoved from the network. Eventually, both methods become essentially the same, oncethe clustering coefficient for the whole grid is uniform.The fastest way to reduce the size of the giant component of the network is, asexpected, to attack the most central nodes, either those with the highest degree orbetweenness; when these bulk nodes are removed, the graph rapidly segregates intomultiple perfectly operational clusters. In other words, the transition to an unconnectedgraph T is always the lowest for these methods.For the case of node elimination based on basin stability it was found that the ecreased resilience in power grids under dynamically induced vulnerabilities K T Random∆ θ e k K s x (a) First phase transition. K T Random∆ θ e k K s x (b) Second phase transition. Figure 14: Phase transitions in the edge percolation process for synthetic power gridsas a function of K . Insets present the standard deviation s x for each computed datapoint.transition to blackout T in the randomly grown synthetic power grids did not differsignificantly to that found for the random removal approach. However, in the specificcase of the Colombian power grid, this method proved to be the most effective to bringthe whole network to total blackout. Essentially, the only key difference between bothtest cases is the power distribution: for the synthetic grids the power was initiallydistributed uniformly, with an equal number of generators and consumers, but for theColombian 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 thegraph, and this was confirmed when it was observed that for the Colombian grid, thefirst nodes to be removed are generators, leading to a higher stress in the remainingpower plants: a few generators will have to increase drastically the power they producein order to balance the demand. In addition, this power increase in the generators willreduce their basin stability, that is, it will make them more prone to dynamical randomdisturbances in oscillatory frequency and phase angle. This negative correlation betweenthe power parameter of a node and its non-linear stability is supported by the findings in ecreased resilience in power grids under dynamically induced vulnerabilities θ ) was found to be the most efficient inorder to bring the whole network to blackout, and remarkably, this behaviour can beobserved 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 causesthat are yielding the reduction of non-linear stability due to heterogeneous powerdistribution. Similarly, a more detailed study on the phase transitions T and T , aswell as the influence of structural parameters in these critical points, could potentiallypoint in the right direction to design the most appropriate connectivity in a graph andparameter distribution such that resilience against cascading failures is optimized. Moredynamical measures of vulnerability could also be incorporated in the removal algorithm;for instance, the diffusivity of perturbations [10] that could indicate regions of fasterfailure propagation in the grid, or the multi-node basin stability [23] computed over thenearest-neighbours, as it would include information about the most vulnerable clusters,rather than individual elements. This would however impose a higher computationaleffort to a task that is already computationally expensive.
5. Acknowledgements
We would like to pay gratitude to Jorge Fernando Guti´errez (deceased on september2019) for valuable ideas and advice at the initial stages of this research. We are alsothankful with Arthur Montanari, Luis A. Aguirre, Laura Lotero, Fabiola Angulo andJuan Gabriel Restrepo for fruitful discussions and the technical assistance of ReinelTabares and Rafael Ruiz. C. C. Galindo-Gonz´alez received financial support from
Universidad Nacional de Colombia - Manizales, Direcci´on de Investigaci´on y Extensi´on(DIMA) and
Ministerio de Ciencia, Tecnolog´ıa e Innovaci´on (Minciencias) under thecontracts FP44842 - 052 - 2016 and FP44842 - 505 - 2017. D. Angulo-Garc´ıa wasfinancially supported by “Vicerrectoria de Investigaciones - Universidad de Cartagena”through the project 085-2018. ecreased resilience in power grids under dynamically induced vulnerabilities Appendix A. Random growth model for power grids
This section describes the algorithm proposed in [39] which allows to generate complexnetworks that incorporate a spatial embedding (accounts for geographic position ofnodes) and that have similar structural properties to those found in real-world powergrids. For this algorithm, a redundancy cost function that will be subject to optimizationhas to be defined as: f ( i,j,G ) = ( d G ( i, j ) + 1) u d S ( y i , y j ) (A.1)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. Afterdefining the parameters from table 1, the construction algorithm is developed throughtwo phases: initialization and growth, which are performed as follows: Appendix A.1. Initialization (i) Initialize a minimum spanning tree for the initial N nodes such that it optimizesthe edges that have to be built between them based on the minimum Euclideandistance d S ( y i , y j ).(ii) Let m = (cid:98) N (1 − s )( p + q ) (cid:99) . For each h = { , , ..., m } find the pair of nodes ( i, j )that are not yet linked and for which f ( i,j,G ) is maximal and connect them. Appendix A.2. Growth (i) For each h = { , , ..., N − N } do:(a) Add a new node h to the graph.(b) Generate a random number ζ . If ζ ≤ − s then:– Set the position of the new node y h .– 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 ζ , if ζ < 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 ζ , if ζ < q draw a node h (cid:48) from the networkuniformly at random. Then find the node l (cid:48) that is not yet linked to h (cid:48) and for which f ( h,l,G ) is maximal and connect them.(c) Otherwise (if ζ > − s ), select an edge ( i, j ) of the network uniformly atrandom. The new geographic location of node h will be given by y h = ( y i + y j ) / ,so the edge ( i, j ) is removed from the graph and then the edges ( i, h ) and ( h, j )are added. ecreased resilience in power grids under dynamically induced vulnerabilities References [1] Pikovsky A, Rosenblum M and Kurths J 2001
Synchronization: A Universal Concept in NonlinearSciences
Cambridge Nonlinear Science Series (Cambridge University Press) URL https://doi.org/10.1017/CBO9780511755743 [2] Arenas A, D´ıaz-Guilera A, Kurths J, Moreno Y and Zhou C 2008
Physics Reports
93 – 153 ISSN 0370-1573 URL [3] Osipov G, Kurths J and Zhou C 2007
Synchronization in Oscillatory Networks (Springer Seriesin Synergetics) ISBN 978-3-540-71268-8 URL https://link.springer.com/book/10.1007/978-3-540-71269-5 [4] D¨orfler F, Chertkov M and Bullo F 2013
Proceedings of the National Academy of Sciences [5] Motter A E, Myers S A, Anghel M and Nishikawa T 2013
Nature Physics Preprint ) URL https://doi.org/10.1038/nphys2535 [6] Sakaguchi H and Matsuo T 2012
Journal of the Physical Society of Japan https://doi.org/10.1143/JPSJ.81.074005 [7] Buldyrev S V, Parshani R, Paul G, Stanley H E and Havlin S 2010 Nature http://dx.doi.org/10.1038/nature08932 [8] Brockway P E, Owen A, Brand-Correa L I and Hardt L 2019
Nature Energy https://doi.org/10.1038/s41560-019-0425-z [9] Xi K, Dubbeldam J L A, Lin H X and Schuppen J H V 2018 Automatica https://doi.org/10.1016/j.automatica.2018.02.019 [10] Tamrakar S, Conrath M and Kettemann S 2018 Scientific Reports Preprint ) URL http://dx.doi.org/10.1038/s41598-018-24685-5 [11] Nishikawa T and Motter A E 2015
New Journal of Physics ISSN 13672630 (
Preprint ) URL https://doi.org/10.1088/1367-2630/17/1/015012 [12] Filatrella G, Nielsen A H and Pedersen N F 2008
European Physical Journal B Preprint ) URL https://doi.org/10.1140/epjb/e2008-00098-8 [13] Kuramoto Y 1975 Self-entrainment of a population of coupled non-linear oscillators
InternationalSymposium on Mathematical Problems in Theoretical Physics ed Araki H (Berlin, Heidelberg:Springer Berlin Heidelberg) pp 420–422 ISBN 978-3-540-37509-8[14] Acebr´on J A, Bonilla L L, P´erez Vicente C J, Ritort F and Spigler R 2005
Rev. Mod. Phys. (1)137–185 URL https://link.aps.org/doi/10.1103/RevModPhys.77.137 [15] D¨orfler F and Bullo F 2011 SIAM Journal on Applied Dynamical Systems Preprint https://doi.org/10.1137/10081530X ) URL https://doi.org/10.1137/10081530X [16] Rohden M, Sorge A, Witthaut D and Timme M 2014
Chaos: An Interdisciplinary Journal ofNonlinear Science https://doi.org/10.1063/1.4865895 [17] Rohden M, Sorge A, Timme M and Witthaut D 2012 Phys. Rev. Lett. (6) 064101 URL https://link.aps.org/doi/10.1103/PhysRevLett.109.064101 [18] Olmi S, Navas A, Boccaletti S and Torcini A 2014
Phys. Rev. E (4) 042905 URL https://link.aps.org/doi/10.1103/PhysRevE.90.042905 [19] Xi K, Dubbeldam J L and Lin H X 2017 Chaos ISSN 10541500 (
Preprint ) URL http://dx.doi.org/10.1063/1.4973770 [20] Menck P J, Heitzig J, Marwan N and Kurths J 2013
Nature Physics http://dx.doi.org/10.1038/nphys2516 [21] Menck P J, Heitzig J, Kurths J and Joachim Schellnhuber H 2014 Nature Communications https://doi.org/10.1038/ncomms4969 [22] Nitzbon J, Schultz P, Heitzig J, Kurths J and Hellmann F 2017 New Journal of Physics https://doi.org/10.1088%2F1367-2630%2Faa6321 [23] Mitra C, Choudhary A, Sinha S, Kurths J and Donner R V 2017 Phys. Rev. E (3) 032317 URL ecreased resilience in power grids under dynamically induced vulnerabilities https://link.aps.org/doi/10.1103/PhysRevE.95.032317 [24] Mitra C, Kittel T, Choudhary A, Kurths J and Donner R V 2017 New Journal of Physics https://doi.org/10.1088%2F1367-2630%2Faa7fab [25] Wolff M F, Lind P G and Maass P 2018 Chaos: An Interdisciplinary Journal of Nonlinear Science https://doi.org/10.1063/1.5040689 [26] Schultz P, Heitzig J and Kurths J 2014 New Journal of Physics https://doi.org/10.1088%2F1367-2630%2F16%2F12%2F125001 [27] Plietzsch A, Schultz P, Heitzig J and Kurths J 2016 The European Physical Journal Special Topics https://doi.org/10.1140/epjst/e2015-50137-4 [28] Montanari A N, Moreira E I and Aguirre L A 2020
Communications in Nonlinear Science andNumerical Simulation https://doi.org/10.1016/j.cnsns.2020.105296 [29] Simonsen I, Buzna L, Peters K, Bornholdt S and Helbing D 2008 Physical review letters https://doi.org/10.1103/physrevlett.100.218701 [30] Rohden M, Jung D, Tamrakar S and Kettemann S 2016
Phys. Rev. E (3) 032209 URL https://link.aps.org/doi/10.1103/PhysRevE.94.032209 [31] Sch¨afer B, Witthaut D, Timme M and Latora V 2018 Nature Communications ISSN 20411723(
Preprint ) URL http://dx.doi.org/10.1038/s41467-018-04287-5 [32] Moreno A M 2017
An´alisis de estabilidad transitoria basado en teor´ıa de redes complejas yel fen´omeno de percolaci´on
Master’s thesis Universidad Nacional de Colombia URL https://repositorio.unal.edu.co/handle/unal/62236 [33] Fu T, Zou L, Li C and Zhao J 2017
Physics Letters, Section A: General, Atomic and Solid StatePhysics http://dx.doi.org/10.1016/j.physleta.2017.06.005 [34] Li D, Zhang Q, Zio E, Havlin S and Kang R 2015
Reliability Engineering and System Safety http://dx.doi.org/10.1016/j.ress.2015.05.021 [35] Huang Z, Wang C, Ruj S, Stojmenovic M and Nayak A 2013 Modeling cascading failuresin smart power grid using interdependent complex networks and percolation theory pp 1023–1028 URL https://ieeexplore.ieee.org/document/6566517 [36] UPME 2016
STN - Sistema de Transmisi´on Nacional de Energ´ıa El´ectrica, Colombia
URL http://sig.simec.gov.co/GeoPortal/Mapas/Mapas [37] Saram¨aki J, Kivel¨a M, Onnela J P, Kaski K and Kert´esz J 2007
Phys. Rev. E (2) 027105 URL https://link.aps.org/doi/10.1103/PhysRevE.75.027105 [38] Brandes U 2001 The Journal of Mathematical Sociology https://doi.org/10.1080/0022250X.2001.9990249 [39] Schultz P, Heitzig J and Kurths J 2014 The European Physical Journal Special Topics https://doi.org/10.1140/epjst/e2014-02279-6 [40] Kuramoto Y 1984
Chemical oscillations, waves, and turbulence (Springer, Berlin, Hei-delberg) ISBN 978-3-642-12600-0 URL https://link.springer.com/book/10.1007/978-3-642-69689-3 [41] Witthaut D and Timme M 2012
New Journal of Physics https://doi.org/10.1088%2F1367-2630%2F14%2F8%2F083036 [42] Rohden M, Witthaut D, Timme M and Meyer-Ortmanns H 2017 New Journal of Physics https://doi.org/10.1088%2F1367-2630%2Faa5597 [43] Schultz P, Menck P J, Heitzig J and Kurths J 2017 New Journal of Physics https://doi.org/10.1088%2F1367-2630%2Faa5a7b [44] Kim H, Lee S H and Holme P 2016 Phys. Rev. E (6) 062318 URL https://link.aps.org/doi/10.1103/PhysRevE.93.062318 [45] Kim H, Lee S H, Davidsen J and Son S W 2018 New Journal of Physics113006 URL