Decreased resilience in power grids under dynamically induced vulnerabilities

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 ${\dot {\phi }}_{i}{:=}\frac{\mathrm{d}{\phi }_{i}}{\mathrm{d}t}$, 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 ${\dot {\theta }}_{i}$ can be described by the well-known second-order Kuramoto model for coupled oscillators given by [12, 16]:

$$\begin{equation}{\ddot {\theta }}_{i\left(t\right)}={P}_{i}-{\alpha }_{i}{\dot {\theta }}_{i\left(t\right)}+K\sum _{j}^{N}{a}_{ij}\enspace \mathrm{sin}\left({\theta }_{j\left(t\right)}-{\theta }_{i\left(t\right)}\right)\end{equation} \tag{ 1 }$$

where P_i is equal to the injected (drained) power at the ith node up to a scaling factor, and it is positive (negative) for generators (consumers); α_i is related to the damping coefficient of the ith 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 $\vert {\dot {\theta }}_{i}\vert \ll {\Omega}$, ∀ 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 ${\dot {\theta }}_{i\left(t\right)}={\dot {\theta }}_{j\left(t\right)}={\omega }_{\mathrm{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 [46], 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:

$$\begin{equation}\sum _{i=1}^{N}{P}_{i}=0.\end{equation} \tag{ 2 }$$

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 [23]. Panels (b)–(e) present distributions for some classic centrality measures of the graph: node degree, clustering coefficient [47], node betweenness and edge betweenness centrality [48]. Panels (f) and (g) display the behavior of θ_(t) and ${\dot {\theta }}_{\left(t\right)}$ 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 [49], which is specially useful as it creates spatially embedded networks that emulate connectivity distributions of real-world 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 B. The parameters used for the growth model algorithm were chosen to match those used on [23], 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).

Zoom In Zoom Out Reset image size

Table 1. Parameters for the construction of synthetic power grids.

Parameter	Value
N0: an initial amount of nodes, N0 ⩾ 1.	1
N: final amount of nodes that the network will have N > N0.	102
y = {y1, y2, ..., yN }: geographical location of nodes.	yk ∼ U[0,1]
p: probability of constructing additional redundancy links attached to new nodes (0 ⩽ p ⩽ 1).	$\frac{1}{5}$
q: probability of constructing additional redundancy links between existing nodes (0 ⩽ q ⩽ 1).	$\frac{3}{10}$
s: probability of splitting an existing line (0 ⩽ s ⩽ 1).	$\frac{1}{10}$
u: cost-vs-redundancy trade-off parameter for equation (B.1).	$\frac{1}{3}$

A helpful indicator that quantifies the degree of synchronization is the complex-valued order parameter r defined as [50]:

$$\begin{equation}r\enspace {\text{e}}^{\text{i}{{\Psi}}_{\left(t\right)}}=\frac{1}{N}\sum _{j}^{N}{\text{e}}^{\text{i}{\theta }_{j\left(t\right)}}.\end{equation} \tag{ 3 }$$

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]:

$$\begin{equation}{r}_{\infty }=\underset{{t}_{1}\to \infty }{\mathrm{lim}}\left[\underset{{t}_{2}\to \infty }{\mathrm{lim}}\left(\frac{1}{{t}_{2}}{\int }_{{t}_{1}}^{{t}_{1}+{t}_{2}}{r}_{\left(t\right)}\mathrm{d}t\right)\right].\end{equation} \tag{ 4 }$$

Another measure of frequency synchronization in power grids is the squared average rotational speed ${v}_{\left(t\right)}^{2}$ defined as [16]:

$$\begin{equation}{v}_{\left(t\right)}^{2}=\frac{1}{N}\sum _{j}^{N}{\dot {\theta }}_{j\left(t\right)}^{2}\end{equation} \tag{ 5 }$$

and the associated steady state value of the speed v_∞:

$$\begin{equation}{v}_{\infty }=\sqrt{\underset{{t}_{1}\to \infty }{\mathrm{lim}}\left[\underset{{t}_{2}\to \infty }{\mathrm{lim}}\left(\frac{1}{{t}_{2}}{\int }_{{t}_{1}}^{{t}_{1}+{t}_{2}}{v}_{\left(t\right)}^{2}\mathrm{d}t\right)\right]}.\end{equation} \tag{ 6 }$$

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:

$$\begin{equation}{\theta }^{{\ast}}\approx \frac{1}{K}{L}^{{\dagger}}P\end{equation} \tag{ 7 }$$

where L^† is the pseudo-inverse of the Laplacian matrix, which is computed through the adjacency matrix as $L=\mathrm{diag}\left({\sum }_{j=1}^{n}{a}_{ij}\right)-A$. Let us now define the oriented incidence matrix $B=\left\{{b}_{ij}\right\}\in {\mathbb{R}}^{N{\times}M}$, with:

$$\begin{equation}{b}_{ij}=\begin{cases}1\quad \hfill & \mathrm{i}\mathrm{f}\enspace \mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\enspace i\enspace \mathrm{i}\mathrm{n}\mathrm{c}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\enspace \mathrm{i}\mathrm{n}\enspace \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}\enspace j\hfill \\ -1\quad \hfill & \mathrm{i}\mathrm{f}\enspace \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}\enspace j\enspace \mathrm{i}\mathrm{n}\mathrm{c}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\enspace \mathrm{i}\mathrm{n}\enspace \mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\enspace i\hfill \\ 0\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{cases}.\end{equation} \tag{ 8 }$$

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:

$$\begin{equation}{{\Delta}}_{\theta }\approx \frac{1}{K}{B}^{\mathrm{T}}{L}^{{\dagger}}P\end{equation} \tag{ 9 }$$

with ${{\Delta}}_{\theta }\in {\mathbb{R}}^{M}$. Then, the system described by equation (1) is guaranteed to have a unique and stable solution with synchronized frequencies and cohesive phases if the next sufficient condition for Δ_θ holds:

$$\begin{equation}{\Vert}{{\Delta}}_{\theta }{{\Vert}}_{\infty }{< }\mathrm{sin}\left(\gamma \right)\end{equation} \tag{ 10 }$$

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 $\frac{\pi }{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):

$$\begin{equation}{K}_{\mathrm{c}}\approx {\Vert}{B}^{\mathrm{T}}{L}^{{\dagger}}P{{\Vert}}_{\infty }.\end{equation} \tag{ 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 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.

Zoom In Zoom Out Reset image size

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 [51, 52], 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).

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₍₀₎ that asymptotically converge to X* [21]. 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 [21, 22]:

$$\begin{equation}{S}_{\mathrm{B}\left(\beta \right)}=\int {{\Gamma}}_{\beta \left(x\right)}{\rho }_{\left(x\right)}\mathrm{d}x\end{equation} \tag{ 12 }$$

where Γ_β(x) is a function that indicates whether a state x belongs to the basin of attraction of X* or not, that is

$$\begin{equation}{{\Gamma}}_{\beta \left(x\right)}=\begin{cases}1\quad \forall \enspace x\in \beta \quad \hfill \\ 0\quad \forall \enspace x\notin \beta \hfill \end{cases} .\end{equation} \tag{ 13 }$$

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:

$$\begin{equation}{\rho }_{\left(x\right)}=\begin{cases}\frac{1}{\left\vert {\Pi}\right\vert }\quad \hfill & \forall \enspace x\in {\Pi}\hfill \\ 0\quad \hfill & \forall \enspace x\notin {\Pi}\hfill \end{cases}\end{equation} \tag{ 14 }$$

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}_{\mathrm{B}}\approx \frac{{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 single-node basin stability (SNBS) [22, 24], given by:

$$\begin{equation}\begin{matrix}\hfill {S}_{\mathrm{B}}^{\left(i\right)}=\frac{{F}_{C}}{{I}_{C}},\quad i\in \vartheta \hfill \end{matrix}\end{equation} \tag{ 15 }$$

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].

Zoom In Zoom Out Reset image size

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:

• Degree centrality ${d}_{k}^{\left(i\right)}$: Amount of nodes to which node i is connected, which in an undirected graph is simply given by the sum of the elements in the ith row of the adjacency matrix A [53]. Here, a normalized version of this metric is used, so it is divided by the maximum possible degree that a node could have in the network, which is evidently N − 1 provided that self-loops are not allowed. It is thus computed as:
  $$\begin{equation}{d}_{k}^{\left(i\right)}=\frac{1}{N-1}\sum _{j}^{N}{a}_{ij}\end{equation} \tag{ 16 }$$
• Clustering coefficient ${c}_{k}^{\left(i\right)}$: Number of triangles in which node i is involved, normalized by the maximum possible amount of such triangles [47]. A triangle involves node i when two neighbors of such node, namely j and k, are also directly connected, so more formally, if (i, j), (j, k), (i, k) ∈ ɛ, then i, j and k form a triangle. Then the clustering coefficient is simply computed as:
  $$\begin{equation}{c}_{k}^{\left(i\right)}=\frac{2{W}_{i}}{{\hat{d}}_{k}^{\left(i\right)}\left({\hat{d}}_{k}^{\left(i\right)}-1\right)}\end{equation} \tag{ 17 }$$
  where W_i is the total number of triangles involving node i and ${\hat{d}}_{k}^{\left(i\right)}={\sum }_{j}^{N}{a}_{ij}$ is the degree centrality without the normalization factor.
• Node betweenness centrality ${b}_{k}^{\left(i\right)}$: Sum of the shortest paths between every pair of nodes (s, t) in the network that pass through node i. For clarity, let us define a path that joins s and t as a sequence of edges H_m = {(s, v₁), (v₁, v₂), (v₂, v₃)..., (v_j , t)}, (v_i , v_k ) ∈ ɛ where s and t are the beginning and end nodes, respectively. Then, by considering the length of the path as the number of edges contained in H_m , the shortest path is the sequence H_m that yields the minimum possible length. Node betweenness centrality is then given by [48, 54]:
  $$\begin{equation}{b}_{k}^{\left(i\right)}=\sum _{s\ne t\ne i\in \vartheta }\frac{{\lambda }_{s,t}\left(i\right)}{{\lambda }_{s,t}}\end{equation} \tag{ 18 }$$
  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.
• Edge betweenness centrality ${e}_{k}^{\left(l\right)}$: 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:
  $$\begin{equation}{e}_{k}^{\left(l\right)}=\sum _{s\ne t\in \vartheta }\frac{{\lambda }_{s,t}\left(l\right)}{{\lambda }_{s,t}},\quad l\in \varepsilon .\end{equation} \tag{ 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.

Figure 6(a) shows the fraction of nodes in the largest cluster of the network (also called the giant component of the graph and denoted by C) subject to node removal with focal and random attacking policies. Similarly, figure 6(b) presents the fraction of removed nodes after each iteration with respect to the original size (named removed component and denoted by $\frac{\left({N}_{\left(0\right)}-{N}_{\left(I\right)}\right)}{{N}_{\left(0\right)}}$, being N₍₀₎ and N_(I) the original size of G and its size after I removal iterations, respectively), which allows to detect the phase transition T₂, corresponding to the iteration at which $\frac{\left({N}_{\left(0\right)}-{N}_{\left(I\right)}\right)}{{N}_{\left(0\right)}}=1$. 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.

Zoom In Zoom Out Reset image size

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 neighbors 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 neighbors 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₂, which reveals the location of the phase transition T₁ as the first iteration when C₂ 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 [22], where nodes located on dead-tree arrangements (more specifically, inner tree nodes [23]), 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.

Table 2. Average value of T₁ and T₂ for the nodal resilience testing algorithm in the Colombian grid.

	Random node	bk	ck	dk	${S}_{\mathrm{B}}^{k}$
T1	14.41	3.00	28.77	3.18	20.60
T2	61.92	40.85	65.03	40.54	29.80

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. It is also worth mentioning that the same experiment was repeated for synthetic power grids with the same power distribution as the Colombian power grid (that is, 28 generators and 74 consumers), and the results achieved were virtually identical to those presented in figure 6, thus they are not presented here, yet this allows us to safely conclude that the unique difference between the specific Colombian case and the random networks is the power distribution. This confirms the fact that the algorithm proposed in [49] is indeed capable of reproducing real-world examples of networked power systems, but also reveals that a 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₁ and T₂ on the synthetic power grids are presented in table 3.

Zoom In Zoom Out Reset image size

Table 3. Average value of T₁ and T₂ for the resilience testing algorithm in the synthetic power grids.

	Random node	bk	ck	dk	${S}_{\mathrm{B}}^{k}$
T1	10.20	2.35	19.58	2.50	27.00
T2	61.62	39.91	63.22	39.02	54.30

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₂ 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₂ in the S_B-focused method is the lowest value in table 2.

Zoom In Zoom Out Reset image size

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}_{\mathrm{B}}^{\left(i\right)}$. After performing a Pearson correlation testing between ${S}_{\mathrm{B}}^{\left(i\right)}$ and P⁽ⁱ⁾, 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 [29]. 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 favored, 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.

Zoom In Zoom Out Reset image size

Let us now consider weak nodes as those that exhibit an SNBS lower than 0.4, then the fraction of weak nodes can be expressed as:

$$\begin{equation}{\eta }_{\mathrm{w}\left(I\right)}=\frac{{N}_{\mathrm{w}\left(I\right)}}{{N}_{\left(I\right)}}\end{equation} \tag{ 20 }$$

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 behavior 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.

Zoom In Zoom Out Reset image size

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₁ 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 500 simulations 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 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 [56, 57]). For the random and the topology-dependant elimination methods, no tendency is observed for T₁ when varying K.

Zoom In Zoom Out Reset image size

T₂ 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 topological-based elimination methods seem to reach a plateau when K is sufficiently high. This reduction in the variability of T₂ 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 Π.

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.

Zoom In Zoom Out Reset image size

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 giant component reduction profile observed for random and e_k methods is in agreement with the results presented in [40]. The values found for T₁ and T₂ on these simulations are summarized on table 4.

Zoom In Zoom Out Reset image size

Finally, for the line resilience procedure over synthetic power grids, the behavior of the transitions is observed on figure 14. For T₁ 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₂: 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 [56, 57].

Zoom In Zoom Out Reset image size

• (a)  
  Initialize a minimum spanning tree for the initial N₀ 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 ).
• (b)  
  Let m = ⌊N₀(1 − s)(p + q)⌋. For each h = {1, 2, ..., 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.

• (a)  
  For each h = {1, 2, ..., N − N₀} do:

  • 1.  
    Add a new node h to the graph.
  • 2.  
    Generate a random number ζ₁. If ζ₁ ⩽ 1 − s then:

    • (i)  
      Set the position of the new node y_h .
    • (ii)  
      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.
    • (iii)  
      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.
    • (iv)  
      Generate a random number ζ₃, if ζ₃ < 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 which f_(h,l,G) is maximal and connect them.
  • 3.  
    Otherwise (if ζ₁ > 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}=\frac{\left({y}_{i}+{y}_{j}\right)}{2}$, so the edge (i, j) is removed from the graph and then the edges (i, h) and (h, j) are added.