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

The synchrony of electric power systems is important in order to maintain stable electricity supply. Recently, the measure basin stability was introduced to quantify a node's ability to recover its synchronization when perturbed. In this work, we focus on how basin stability depends on the coupling strength between nodes. We use the Chilean power grid as a case study. In general, basin stability goes from zero to one as coupling strength increases. However, this transition does not happen at the same value for different nodes. By understanding the transition for individual nodes, we can further characterize their role in the power-transmission dynamics. We find that nodes with an exceptionally large transition window also have a low community consistency. In other words, they are hard to classify to one community when applying a community detection algorithm. This also gives an efficient way to identify nodes with a long transition window (which is computationally time consuming). Finally, to corroborate these results, we present a stylized example network with prescribed community structures that captures the mentioned characteristics of basin stability transition and recreates our observations.


Introduction
Electric power systems are often modeled as oscillators on networks [1,2,3]. The rotational motion of power generators for alternating currents needs to be synchronized with a rated frequency of an electric power grid [4,5,6,7,8,9]. Stable synchrony of electric power grids' nodes is important to prevent cascading failures. When a node gets an external perturbation leading the node away from a synchronous state, then the node either returns to the synchronous state or escapes from its basin of attraction, typically to a different limit cycle. Basin stability is calculated as the fraction of the possible phase values a node can be perturbed to and still recover synchronicity [6]. It has become the standard way to quantify a power-grid node's robustness to large point perturbations. The basin stability is determined by the magnitude of the perturbation and network characteristics of the node. It also depends on the coupling strength between the nodes in the dynamic system representing the electric power flow. For instance, Refs. [7,8] investigate topological characteristics of particularly stable or unstable nodes. The authors conclude that, in general, dead ends weaken the basin stability, while detours strengthen it. Such simple indirect estimates are helpful because basin stability is very computer intensive to calculate.
These previous studies, however, study basin stability for one specific coupling strength [6,7,8]. However, basin stability is sensitive to the coupling strength. In a power grid, the transmission capacity, which determines to the coupling strength, is decided based on the physical conditions such as the air temperature, the length of transmission line, and the amount of transmitting electricity. When the amount of electricity transmitted exceeds the transmission capacity, it may cause the breakdown of a large part of the system. Therefore, it is important to understand the functional dependence of basin stability on the coupling strength, not only its behavior at a fixed and intermediate coupling strength. We use the real electric power grid in Chile as our prime example, but also discuss how the results generalize. To be more specific, we investigate the network-structural factors [10,11] behind a node's basin stability as a function of the coupling strength. Based on the Chilean case study, we argue that there is a connection between the coupling strength dependence and community structure-the way the network can be decomposed into communities that are densely connected within and sparsely connected between each other [12,13]. Other authors have noticed that those communities easily synchronize internally. In a similar spirit, we investigate how the membership strength of a node within a community decomposition (its community consistency) of the network relates to the coupling strength dependence, specifically the width in parameter space of the basin stability's transition from zero to one.

The dynamical model of electricity transmission
The synchronization dynamics between power-grid nodes is commonly modeled as the set of nonlinear oscillators using a Kuramoto-type model [1,2,4,5,6,7,8,9]: 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 i j = 1 if there is an edge between node i and j, and a i j = 0 otherwise; θ i is the phase of node i; ω i =θ 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ω < 5 × 10 −2 and the angular frequency ω < 5 × 10 −2 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ω < 5×10 −2 and ω < 5×10 −2 after the 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 Refs. [6,7,8], we set α i = 0.1 for all i.

The Chilean power grid
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 powergrid 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 Eq. (1)) or net consumers (P i < 0) based on the real data [18]. After determining the sign of P i values, we first set P i = −1 for P i < 0 and P i = P pos > 0 so that i P i = 0 (Kirchhoff's circuit laws). The inactive nodes such as towers or net zero traders (P i = 0 in Eq. (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 the 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.

Basin stability transition
Basin stability is a stability measure for a network [6,7,8]. For each node, an initial phase value is selected uniformly at random. Then we run the dynamics until convergence (according to the criterion presented in Sec. 2.1). The fraction of initial phase values that leads to convergence defines the basin stability. We sample 100 points uniformly at random from the intervals −100 < ω < 100 and −π < θ < π. Assuming a uniformly random sampling is done for the sake of simplicity (following Refs. [6,7,8]). 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 fig. 1. For low K, all nodes have zero basin stability. As K increases, a few nodes start to reach a non-zero basin stability ( fig. 1(a)). When the coupling strength reaches a certain threshold, almost all the nodes suddenly turn into a stable state. Note that 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 Sec. 3, we address the transition of basin stability and specifically introduce basin stability transition window.

Basin stability transition window
Basin stability changes, as mentioned, with the coupling strength. Figure 2 shows a schematic diagram of the basin stability transition pattern that we call transition window, for two nodes of a hypothetical network. At K = K 0 , node 1 has larger basin stability than node 2. On the other hand, when the coupling strength increases to K = K 1 , 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 on the transition window (the shaded rectangles of fig. 2).
To locate the transition windows of the nodes, we use the bisection method [24]. 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 the 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 fig. 3(a) undergoes the 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 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 fig. 3(b)) and the range of average [K low , K high ] for each bin (the right panel of fig. 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 1 -10 2 ) is more than 13 times as large as that for the bin corresponding to the minimum K range (10 −3 -10 −2 ). In general, the wider of the transition window a node is, the larger K value in the middle ((K low + K high )/2, henceforth denoted as K mid ) the node has.
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 the similar ∆K 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 started to decrease, the node suddenly loses the 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 the wider window are better for an early detection of system failure.

Chilean power-grid analysis
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] predict that dead-end nodes have small basin stability, while Schultz et al. [8] predict that detour nodes have large basin stability.
The basin stability estimators in Refs. [7] and [8] assume a specific K-value. If K is not precisely known or varying, it can only identify extreme cases (very stable or unstable nodes). The width of 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 on 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 on networks [25,26,27]. 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 [28], a variant of conventional Louvain community detection algorithm [29] 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 [30,31,32]. We define community consistency of individual nodes as the nodes' degree of certainty of the assigned community. Previously, Refs. [14,15] calculate community consistency (or "consensus clustering") among assigned communities in multiple runs of each stochastic community detection algorithm, where they focus on the community consistency of each algorithm. In this paper, however, we measure 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 the we use the same resolution parameter [30,31,32]γ = 0.0323. We chose this value to get around three communities which matches our visual impression of the network. As hinted in fig. 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 fig. 4). In terms of power-grid dynamics, such nodes could, be believe, be influenced by perturbations from different directions (communities).
To measure the community consistency of individual nodes, we first observe a node i is perfectly community consistent if, and only if, it gets always classified into a cluster with the same other nodes. Let φ i j 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 φ i j is either 1 or 0 for all j. Since φ i j ∈ [0, 1], we can obtain a metrics for community consistency by measuring the average distance to φ i j = 1/2. But rather than the linear distance (|φ i j − 1/2|), we sum the square distance. First, this conforms our measure to standard measures of deviation or spread like the 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 a factor 4 to the outlined measure to get a value in the unit interval. In summary, the community consistency of i is given by The correlation between community consistency and ∆K of basin stability transition window of nodes is shown in fig. 5 and Table 1, where the larger 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 is colored, which provides a visual evidence of our result of 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 basin stability transition window. Moreover, compared to the actual numerical integration for calculating the basin stability, the calculating community consistency is much faster (at least for reasonably fast community detection algorithms).
To check if there is any other (possibly simpler) network measures that can predict ∆K, we measure other network metrics. 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 −3 ), 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.

Example network analysis
In order to verify our results from the Chilean power grid that the basin stability transition is closely related to the community membership (figs. 8(a) and (b)) we construct a simple example network with prescribed community structures depicted in fig. 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 Eq. (1)) and the circular nodes are producers (P i = 1 in Eq. (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 fig. 8(a)). We measure basin stability transition window of the nodes in the example network and show representative cases in fig. 8(d). Most characteristics of the transition patterns observed from various nodes in the Chilean power-grid nodes ( fig. 8(b)) are captured by some representative nodes in the example network ( fig. 8(d)), e.g., concave versus convex transitions and the staircase pattern.
To be more specific, the basin stability transition curve becomes convex or shows the 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 Chilean power grid and our example network. Such outlier nodes seem to have a tendency of reaching a stable state (the basin stability equals to unity) at larger values of K compared to the nodes with stronger community membership. In other words, our example network can be considered as the 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 (figs. 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 to the terminal node, the transition curve changes its shape form from concave down to concave up. However, the actual terminal node G in the Chilean power grid does not follow this trend and reaches 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 ( fig. 8(d)) shows a similar pattern. Those nodes are commonly located at the terminus of the interface branch, and we speculate that this pattern is related to 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 ( fig. 8(d)). For example, nodes 1 and 6 in the example network (figs. 8(c) and (d)) has 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 synchronization on networks [3].

Summary and discussions
We have studied the basin stability transition throughout the coupling strength parameter space, focusing on the basin stability transition window as a new metrics for characterizing the contribution of a node to the stability of a power grid. As previous works focus on the stability for specific coupling strength values, our approach complements these. While a narrow transition window implies a sudden change in the stability, a wide transition window can provide an early signal of danger when the coupling strength is gradually weakened. By comparing the mesoscopic network property of community consistency with the transition window, we found that the former is a good predictor of the latter (or vice versa), signifying the importance of community structures on the synchronization dynamics again. On a practical side, such community-consistency based predictions provide a proxy of the actual time-consuming simulations of calculating basin stability, especially for very large systems. Furthermore, once we assign the communities and the strength of nodes' memberships by community consistency, we are able to analyze the dynamics in the unit of such communities instead of the entire node set. We would also like to emphasize that the network example with given community structures and the bridge effectively captures the basin stability transition properties observed in the Chilean power grid.
Our approach can be improved by considering more realistic situations such as assigning different parameters of real fluctuations of power input and output, different edge weights based on admittance of the transmission line, etc. We see our work as a starting point, and anticipate more studies in the future. Finally, we would like to emphasize that our communityconsistency measure is not limited to the power-grids, but could be used in all kinds of community detection problems.