Community consistency determines the stability transition window of power-grid nodes

The synchronization dynamics between power-grid nodes are commonly modeled as a set of nonlinear oscillators using a Kuramoto-type model [1, 2, 4–9]:

$$\begin{eqnarray}&&{\ddot{\theta }}_{i}={\dot{\omega }}_{i}={P}_{i}-{\alpha }_{i}{\dot{\theta }}_{i}-K\displaystyle \sum _{j\ne i}{a}_{{ij}}\mathrm{sin}({\theta }_{i}-{\theta }_{j}),\end{eqnarray} \tag{ 1 }$$

where P_i is the net power generation (positive) or consumption (negative) at node i; α_i is a dissipation constant; K is the coupling strength between nodes i and j; the adjacency matrix a_ij = 1 if there is an edge between node i and j, and a_ij = 0 otherwise; θ_i is the phase of node i; ${\omega }_{i}={\dot{\theta }}_{i}$ is the angular velocity of node i. A power grid is considered to be stable and synchronized when ω_i vanishes for all of the nodes so that the system maintains the desirable constant frequency. For numerical integration, we use the Runge–Kutta method [16] with the convergence criterion [7] of the time derivative of angular frequency $\dot{\omega }\lt 5\times {10}^{-2}$ and the angular frequency ω < 5 × 10⁻² based on the actual fluctuations, i.e., we consider the initial phase and angular velocity (θ_i, ω_i) of i as belonging to the stable basin if the system eventually converges to $\dot{\omega }\lt 5\times {10}^{-2}$ and ω < 5 × 10⁻² after i is perturbed. We perturb i by initializing ω_i and θ_i to random values in the intervals [−100, 100] and [−π, π). The other nodes have their values of the synchronous state. Following [6–8], we set α_i = 0.1 for all i.

We construct the Chilean power grid, our prime example, from the operational records published by the National Energy Committee of Chile [17] and the data regarding the power-grid transmission lines published by Centro de Despacho Económico de Carga del Sistema Interconectado Central (CDEC-SIC)—the major electric power company in Chile [18]. The Chilean power grid consists of power plants and substations as nodes, and the transmission lines between them as edges [10, 11]. A power plant generates electricity and a substation distributes it to the final consumers. We aggregate the power grid using Zhukov's aggregation scheme [19], resulting in nodes being either net generators (P_i > 0 in equation (1)) or net consumers (P_i < 0) based on the real data [18]. After determining the sign of the P_i values, we first set P_i = −1 for P_i < 0 and P_i = P_pos > 0 so that ${\displaystyle \sum }_{i}{P}_{i}=0$ (Kirchhoff's circuit laws). By these parameter values, we assume that all of the power plants produce the same amount of electricity and distribute it evenly for all substations with the same transmission loss on each transmission line. In other words, we focus only on the heterogeneity caused by the network structure, as in [6–8]. The inactive nodes such as towers or net zero traders (P_i = 0 in equation (1)) do not intervene with the perturbation propagation and are removed. We apply the so-called Kron reduction [20, 21] procedure. It is a node elimination technique widely used in electrical circuit design. For a given network configuration, Kron reduction removes a node and adds new edges between neighbors of the eliminated node at each step [22, 23]. Consequently, as more nodes are Kron-reduced, the denser (in terms of mean degree) the network becomes. However, even though the network structure is changed, the new edges are placed so that the network maintains the relationship between nodes in terms of electricity flow. For the Chilean power grid, we eliminate 46 inactive nodes using Kron reduction. As a result, the final network consists of 291 net consumers, 129 net generators, and 573 edges. We present more details on the network construction in [24].

To calculate the basin stability of a node i, ω_j and θ_j are initialized to random variables in the intervals [−100, 100] and [−π, π), respectively, for i and zero for the other nodes. Then we run the dynamics until convergence (according to the criterion presented in section 2.1). We repeat this procedure 100 times and estimate the basin stability as the fraction of runs that converge to the same state. The fraction of initial phase values that leads to convergence defines the basin stability. The precision we can achieve is rather limited by the time-consuming basin-stability calculation. Random sampling of ω_i and θ_i is done for the sake of simplicity (following [6–8]). The standard error for the single-node basin stability in our simulation is about 10⁻². Note that the basin stability of a node could be compared to other nodes in the same network, but not necessarily to nodes in other networks.

For the so-called infinite bus-bar model (basically a mean-field approximation of the problem we study), the basin stability increases gradually with the transmission capacity and suddenly becomes unity (stable) [6]. This jump appears when the transmission capacity is large enough for all perturbations or initial shocks from the entire phase space to be absorbed into the system. However, in the more realistic multi-node model, the basin stability and its transition form are different for different nodes. So for a complete picture, every node needs to be analyzed separately.

We show the results for the distribution of basin stabilities of the Chilean power grid's nodes as functions of the coupling strength K in figure 1. For low K, all nodes have zero basin stability. As K increases, a few nodes start to reach a non-zero basin stability (figure 1(a)). When the coupling strength reaches a certain threshold, almost all the nodes suddenly turn into a stable state. The basin stability shows a nonlinear transition as a function of coupling strength. The dramatically different distribution of basin stabilities for the different K reinforces the importance of considering not only a certain basin stability at a K value but also the transition of basin stability corresponding to different K values. In section 3, we address the transition of basin stability and specifically introduce a basin stability transition window.

Zoom In Zoom Out Reset image size

Basin stability changes, as mentioned, with the coupling strength. Figure 2 shows a schematic diagram of the basin stability transition pattern that we call the transition window, for two nodes of a hypothetical network. At K = K₀, node 1 has a higher basin stability than node 2. On the other hand, when the coupling strength increases to K = K₁, node 2 becomes more stable than node 1. These transitions of basin stability of nodes 1 and 2 are represented as a function of K in the transition window (the shaded rectangles of figure 2). Note that the transition window is the range of coupling strength values where the basin stability is between zero and one, whereas the basin stability is the phase space volume of the basin of attraction.

Zoom In Zoom Out Reset image size

To locate the transition windows of the nodes, we use the bisection method [25]. We start by finding the lower bound of the transition window K_low, which is the minimum value of K that makes the basin stability larger than zero. We do this by tracking the basin stability from K_min = K_min,init and K_max = K_max,init, and measure the basin stability of a node at the midpoint K_midpoint = (K_min + K_max) /2. When the basin stability at K_midpoint is non-zero, K_max is set to the current K_midpoint value and the new K_midpoint value is recalculated from the new K_max value. Similarly, when the basin stability at K_midpoint is zero, the new K_min is set to the current K_midpoint. This iteration continues until K_max − K_min < K_threshold and K_low is set to the final values of (K_min + K_max)/2. The upper bound K_high of the transition window (the maximum value of K that makes the basin stability smaller than unity) is also determined in the same fashion. As a result, we obtain the transition width ΔK = K_high − K_low that represents the width of the transition window. We use K_min,init = 0, K_max,init = 20, and K_threshold = 0.01.

ΔK can vary much between the nodes. For instance, node A (one of the nodes with the narrowest ΔK values) of the Chilean power grid shown in figure 3(a) undergoes a sudden basin stability jump at K ≃ 1.3 and reaches the completely stable state at a smaller K value than other nodes, such as nodes B and C as comparison. The distribution of ΔK in the Chilean power grid is highly right-skewed (with a mean of 0.154 and a median of 0.022). Figure 3(b) shows various statistics of the transition window, where we divide the ΔK values into five different logarithmic bins and present the fraction of nodes belonging to each bin (the left panel of figure 3(b)) and the range of average [K_low, K_high] for each bin (the right panel of figure 3(b)). The vast majority (about 96%) of nodes have very small values of ΔK < 0.1 (about 35% of nodes even have ΔK < 0.01). However, there are a few nodes with relatively very large ΔK values, e.g., the average range [K_low, K_high] for the bin corresponding to the maximum K range (10¹–10²) is more than 13 times as large as that for the bin corresponding to the minimum K range (10⁻³–10⁻²). In general, the wider the transition window of a node is, the larger K value in the middle ((K_low + K_high)/2, henceforth denoted as K_mid) the node has.

Zoom In Zoom Out Reset image size

It is interesting to note that the transition curves themselves seem to have different shapes for different nodes, even for the nodes with similar ΔK values. For example, ΔK of nodes B and C have similar values but the former (latter) node shows roughly concave (convex) curves at the transition. The different range of ΔK and the shape can be interpreted in several ways. If a node has small ΔK, the location of K_mid is also small according to our observations, so the synchrony of the node tends to be stable even for small values of coupling strength. In other words, the node stays functional and is hardly affected by the accidental performance drop of the transmission line connected to it. However, at the same time, once the basin stability value starts to decrease, the node suddenly loses its stability. Consequently, it makes the detection of system failure difficult. When the transition of the basin stability is gradual, on the other hand, a system disturbance can be detected by a drop in the basin stability. In summary, narrow transition windows could be good for keeping stability, while wider windows are better for early detection of system failure.

One drawback of basin stability as a robustness measure of nodes in power-grids is that it is computationally demanding. Therefore, it is desirable to find less computationally restrictive, indirect measures to estimate the basin stability. Some previous studies proposed such topological indicators for basin stability—Menck et al [7] predicted that dead-end nodes have low basin stability, while Schultz et al [8] predicted that detour nodes have high basin stability.

The basin stability estimators in [7] and [8] assume a specific K-value. If K is not precisely known or varying, they can only identify extreme cases (very stable or unstable nodes). The width of the transition window is thus a more appropriate measure to capture the role of a node in the transmission dynamics. To understand why some nodes have larger ΔK than others, we try to find the most appropriate explanatory measure. This is far from trivial. Indeed phase synchronization in a network is a complex consequence of the propagation of phase difference and recovery due to the interactions between nodes [3].

It is known that the mesoscopic properties of networks, such as community structure, play a significant role in the synchronization of networks [26–28]. A community is a subnetwork that is more strongly connected within than to the rest of the network [12, 13]. Different community detection methods could divide the same network in different ways, but there would typically be pairs of nodes that always are classified as belonging to the same community. Some community detection methods are non-deterministic and could also give a different community decomposition between different runs of the same algorithm. One could therefore define nodes often classified to the same community as having high community consistency with respect to this algorithm. We use GenLouvain [29], a variant of conventional Louvain community detection algorithm [30] known for its computational scalability, to identify communities in networks. We can also control the resolution parameter in the modularity function used as the objective function to maximize, to detect communities with different scales [31–33]. We define the community consistency of individual nodes as the nodes' degree of certainty of the assigned community. Previously, [14, 15] calculated community consistency (or 'consensus clustering') among assigned communities in multiple runs of each stochastic community detection algorithm, where they focussed on the community consistency of each algorithm. In this paper, however, we measure the community consistency of individual nodes to relate it to the nodes' basin stability transition width ΔK.

Figure 4 shows four different community detection results of the GenLouvain method on the Chilean power grid, where we use the same resolution parameter [31–33]—γ = 0.0323. We chose this value to get around three communities which match our visual impression of the network. As hinted in figure 4, most nodes are assigned to the same community for all iterations. We call such nodes community consistent. Some nodes, however, are assigned to different communities in different runs (the magnified region in figure 4). In terms of power-grid dynamics, such nodes could, we believe, be influenced by perturbations from different directions (communities).

Zoom In Zoom Out Reset image size

To measure the community consistency of individual nodes, we first observe that a node i is perfectly community consistent if, and only if, it is always classified into a cluster with the same other nodes. Let ϕ_ij be the fractions of runs of the community detection algorithm when i and j are classified into the same community, then i is perfectly community consistent if, and only if, ϕ_ij is either 1 or 0 for all j. Since ${\phi }_{{ij}}\in [0,1],$ we can obtain a metric for community consistency by measuring the average distance to ϕ_ij = 1/2. But rather than the linear distance ($| {\phi }_{{ij}}-1/2|$), we sum the square of the distance. First, this conforms our measure to standard measures of deviation or spread like root-mean-square, variance, and radius of gyration. Second, the non-linearity of the parabola accentuates values close to 0 or 1 and tones down the mid-interval values. This turns out to be practical since most ϕ-values are close to 0 or 1. Finally, we multiply the outlined measure by a factor of 4 to get a value in the unit interval. In summary, the community consistency of i is given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{j\ne i}{(1-2{\phi }_{{ij}})}^{2}.\end{eqnarray} \tag{ 2 }$$

The correlation between community consistency and ΔK of basin stability transition window of nodes is shown in figure 5 and table 1, where the larger the values of community consistency a node has, the narrower ΔK the node has. In other words, if a node's community consistency is weak, so that the node is assigned to different communities for each different realization, the node tends to have a wide basin stability transition window. Figure 6 illustrates the nodes where ΔK and community consistency are colored, which provides visual evidence on our result on the correlation between ΔK and community consistency. It seems that the basin stability transition occurs at some characteristic values of coupling strength for different communities of nodes as a unit, and the transition threshold is somehow smeared out for those nodes with weak or inconsistent community membership, which makes ΔK large. There are several nodes in the northern part of the power grid with lower community consistency, but we believe that the effect of the nodes is not strong enough to affect the basin stability transition window. Moreover, compared to the actual numerical integration for calculating the basin stability, the calculation of community consistency is much faster (at least for reasonably fast community detection algorithms).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a comparison, we also test the Clauset–Newman–Moore (CNM) community detection method [34] in addition to GenLouvain. Although CNM is deterministic in principle, there is a reported effect of particular node ordering on the result, so that isomorphic graphs can get different partitions [14]. We therefore shuffle the order of nodes in the Chilean power grid network and use 100 different realizations of the edge list to obtain different sets of community detection results. We find that the CNM community consistency is also negatively correlated with the basin stability transition window, notwithstanding the weaker correlation value (r = −0.290 and p-value < 10⁻⁸) than that from GenLouvain. The two methods give the same qualitative results, but GenLouvain is inherently stochastic, while the CNM stochasticity is imposed from the random node labeling. Therefore, we advocate the GenLouvain method for better prediction of the basin stability transition window of nodes.

To check if there is any other (possibly simpler) network measure that can predict ΔK, we measure other network metrics (figure 7). The Pearson correlation coefficients for ΔK versus degree [10, 11] and versus community consistency are about 0.033 (p-value around 0.5) and −0.581 (p-value less than 10⁻³), respectively. For the correlation coefficient, values for other representative network centralities and community consistency are shown in table 1, indicating that only the mesoscopic measure of community consistency is significantly correlated with ΔK, in contrast to other conventional measures such as degree, clustering coefficient [11], and current flow betweenness centralities (supposedly a more relevant type of betweenness in our power-grid case than the ordinary one) [11]. The two former centralities are microscopic or local, while the latter is macroscopic as it deals with all of the possible pathways for an entire network. We tried other measures than the ones presented in table 1, but we could not find anything statistically significantly correlated with ΔK.

In order to verify our results from the Chilean power grid that the basin stability transition is closely related to the community membership (figures 8(a) and (b)), we construct a simple example network with prescribed community structures depicted in figure 8(c). The example network consists of 18 nodes and most nodes are separated into two communities, each of which has six fully connected nodes and a single attached dead-end node, except for four nodes branched from the bridge between the communities. The square nodes are assigned as consumers (P_i = −1 in equation (1)) and the circular nodes are producers (P_i = 1 in equation (1)). The example network symbolizes a network structure with two well-defined communities and outlier nodes on an interface branch. The interface branch represents a connected subgraph without clear community membership, just like the nodes F and G in the Chilean power grid (shown in figure 8(a)). We measure the basin stability transition window of the nodes in the example network and show representative cases in figure 8(d). Most characteristics of the transition patterns observed from various nodes in the Chilean power-grid nodes (figure 8(b)) are captured by some representative nodes in the example network (figure 8(d)), e.g., concave versus convex transitions and the staircase pattern.

Zoom In Zoom Out Reset image size

To be more specific, the basin stability transition curve becomes convex or shows a staircase pattern when the node is located far from a prescribed community and near to the branch from the bridge (thus with a weak community membership), in both the Chilean power grid and our example network. Such outlier nodes seem to have a tendency to reach a stable state (the basin stability equals unity) at larger values of K compared to the nodes with stronger community membership. In other words, our example network can be considered as a simplified version of the Chilean power grid at least in terms of the basin stability transition. In particular, we find some special transition patterns in both networks. Node G in the Chilean power grid (figures 8(a) and (b)) shows a stability at a smaller value of K compared to other nodes in the bridge. When a node is located near the terminal node, the transition curve changes its shape from concave down to concave up. However, the actual terminal node G in the Chilean power grid does not follow this trend and reaches a medium level of basin stability at relatively small K, and maintains the plateau until the basin stability suddenly reaches unity. Node 6 in our example network (figure 8(d)) shows a similar pattern. These nodes are commonly located at the terminus of the interface branch, and we speculate that this pattern is related to the fact that the nodes cannot spread external perturbation from the neighbor nodes to the rest of the community. There are also some differences between the Chilean network and the example network, for example the spikes from the plateau to the maximum basin stability (figure 8(d)). For example, nodes 1 and 6 in the example network (figures 8(c) and (d)) have maximum basin stability at K ≈ 5 and 10 (node 1) and K ≈ 13 and 17 (node 6), and they become unstable again. We find this phenomenon for the small example network too (not shown). Why this happens is beyond the scope of this paper, but not surprising in the light of the instabilities in the synchronization of networks [3].