Research ArticleAPPLIED SCIENCES AND ENGINEERING

The key player problem in complex oscillator networks and electric power grids: Resistance centralities identify local vulnerabilities

See allHide authors and affiliations

Science Advances  22 Nov 2019:
Vol. 5, no. 11, eaaw8359
DOI: 10.1126/sciadv.aaw8359

Abstract

Identifying key players in coupled individual systems is a fundamental problem in network theory. We investigate synchronizable network-coupled dynamical systems such as high-voltage electric power grids and coupled oscillators on complex networks. We define key players as nodes that, once perturbed, generate the largest excursion away from synchrony. A spectral decomposition of the coupling matrix gives an elegant solution to this identification problem. We show that, when the coupling matrix is Laplacian, key players are peripheral in the sense of a centrality measure defined from effective resistance distances. For linearly coupled systems, the ranking is efficiently obtained through a single Laplacian matrix inversion, regardless of the operational synchronous state. The resulting ranking index is termed LRank. When nonlinearities are present, a weighted Laplacian matrix inversion gives another ranking index, WLRank. LRank provides a faithful ranking even for well-developed nonlinearities, corresponding to oscillator angle differences up to approximately Δθ ≲ 40°.

INTRODUCTION

Because of growing electric power demand, increasing difficulties with building new lines, and the emergence of intermittent new renewable energy sources, electric power systems are more often operated closer to their maximal capacity (1, 2). Accordingly, their operating state, its robustness against potential disturbances, and its local vulnerabilities need to be assessed more frequently and precisely. Furthermore, because electricity markets become more and more integrated, it is necessary to perform these assessments over geographically larger areas. Grid reliability is commonly assessed against n − 1 feasibility, transient stability, and voltage stability, by which one means that a grid is considered reliable if (i) it still has an acceptable operating state after any one of its n components fails, (ii) that acceptable state is reached from the original state following the transient dynamics generated by the component failure, and (iii) the new operating state is robust against further changes under operating conditions such as changes in power productions and loads. This n − 1 contingency assessment is much harder to implement in real-time for a power grid loaded close to its capacity where the differential equations governing its dynamics become nonlinear—the fast, standardly used linear approximation breaks down as the grid is more and more heavily loaded. Nonlinear assessment algorithms have significantly longer runtimes, which makes them of little use for short-time evaluations. In the worst cases, they sometimes even do not converge. Briefly, heavily loaded grids need more frequent and more precise reliability assessments, which are however harder to obtain, precisely because the loads are closer to the grid capacities.

Developing real-time procedures for n − 1 contingency assessment requires new, innovative algorithms. One appealing avenue is to optimize contingency ranking (3) to try and identify a subset of ns < n grid components containing all the potentially critical components. The n − 1 contingency assessment may then focus on that subset only, with a substantial gain in runtime if ns ≪ n. Identifying such a subset requires a ranking algorithm for grid components, following some well-chosen criterion. Procedures of this kind have been developed in network models for social and computer sciences, biology, and other fields, in the context of the historical and fundamental problem of identifying the key players (4–8). They may be, for instance, the players who, once removed, lead to the biggest changes in the other player’s activity in game theory or to the biggest structural change in a social network. That problem has been addressed with the introduction of graph-theoretic centrality measures (9, 10), which order nodes from the most “central” to the most “peripheral”—in a sense that they themselves define. A plethora of centrality indices has been introduced and discussed in the literature on network theory (9, 10), leading up to PageRank (11). The latter ranks nodes in a network according to the stationary probability distribution of a Markov chain on the network; accordingly, it gives a meaningful ranking of websites under the reasonable assumption that web surfing is a random process. Their computational efficiency makes PageRank and other purely graph-theoretic indicators very attractive to identify key players on complex networks. It is thus quite tempting to apply purely graph-theoretic methods to identify fast and reliably key players in network-coupled dynamical systems.

Processes such as web crawling for information retrieval are essentially random diffusive walks on a complex network, with no physical conservation law beyond the conservation of probability. The situation is similar for disease (12) or rumor (13) spreading and for community formation (14) where graph-theoretic concepts of index, centrality, betweenness, coreness, and so forth have been successfully applied to identify tightly bound communities.

Coupled dynamical systems such as complex supply networks (15), electric power grids (16), consensus algorithm networks (17), or more generally network-coupled oscillators (18, 19) are, however, fundamentally different. There, the randomness of motion on the network giving, e.g., the Markov chain at the core of PageRank, is replaced by a deterministic dynamics supplemented by physical conservation laws that cannot be neglected. Pure or partially extended graph-theoretic methods have been applied in vulnerability investigations of electric power grids (20–22) and investigations of cascades of failures in coupled communication and electric power networks (23, 24). They have, however, been partially or totally invalidated by investigations on more precise models of electric power transmission that take fundamental physical laws into account (in this case, Ohm’s and Kirchhoff’s laws) (25, 26). It is therefore doubtful that purely topological graph-theoretic descriptors are able to identify the potentially critical components in deterministic, network-coupled dynamical systems. Purely graph-theoretic approaches need to be extended to account for physical laws (20). The influence of the dynamics on transient performance for regular graphs on d-dimensional tori has been emphasized in (27).

Here, we give an elegant solution to the key player problem for a family of deterministic, network-coupled dynamical systems related to the Kuramoto model (18, 19). While we focus mostly on high-voltage electric power grids whose swing dynamics, under the lossless line approximation, is given by a second-order version of the Kuramoto model (16, 28), we show that our approach also applies to other generic models of network-coupled oscillators. Key players in these systems can be defined in various ways. For instance, they can be identified by an optimal geographical distribution of system parameters such as inertia, damping, or natural frequencies or alternatively as those whose removal leads to the biggest change in operating state. Here, we define key players as those nodes where a local disturbance leads to the largest short-time transient network response. In the context of electric power grids, transient stability is the ability of the grid to maintain synchrony under relatively large disturbances such as loss or fluctuations of power generation or of a large load (29). If under such a fault, the system remains in the vicinity of its original state, it has maintained synchrony. There are different measures to quantify the magnitude of the transient excursion, such as nadir and maximal rate of change of the network-averaged frequency (30, 31) or other dynamical quantities such as network susceptibilities (32) and the wave dynamics following disturbances (33). Here, we quantify the total transient excursion through performance measures that are time-integrated quadratic forms in the system’s degrees of freedom (see Materials and Methods and the Supplementary Materials). Transient excursions typically last 10 to 20 s in large, continental power grids, which sets the time scales we are interested in.

Anticipating on results to come, Fig. 1 illustrates the excellent agreement between analytical theory and numerical calculations for these performance measures. Particularly interesting is that in both asymptotic limits of quickly and slowly decorrelating noisy disturbance, the performance measures are simply expressed in terms of the resistance centrality (34, 35), which is a variation of the closeness centrality (9) based on resistance distances (36). This is shown in the insets of Fig. 1. Our main finding is that the resistance centrality is the relevant quantity to construct ranking algorithms in network-coupled dynamical systems.

Fig. 1 Comparison between theoretical predictions and numerical results for both performance measures P1 and P2 defined in Eq. 3.

Each point corresponds to a noisy disturbance on a single node of the European electric power grid sketched in Fig. 2A (see Materials and Methods and the Supplementary Materials) and governed by Eq. 1 with constant inertia and damping parameters. The time-dependent disturbance δPi(t) is defined by an Ornstein-Uhlenbeck noise of magnitude δP0 = 1 and correlation time γτ0 = 4 × 10−5 (red crosses), 4 × 10−4 (cyan), 4 × 10−3 (green), 4 × 10−2 (purple), 4 × 10−1 (black), and 4 (blue). Time scales are defined by the ratio of damping to inertia parameters γ = di/mi = 0.4 s−1, which is assumed constant with di = 0.02 s. The insets show 𝒫1 and 𝒫2 as a function of the resistance distance–based graph-theoretic predictions of Eq. 6 valid in both limits of very large and very short noise correlation time τ0. The limit of short τ0 for 𝒫2 gives a node-independent result (Eq. 6).

RESULTS

We consider network-coupled dynamical systems defined by sets of differential equations of the formmiθ¨i+diθ̇i=Pi−∑jbijsin (θi−θj),i=1,…,n(1)

The coupled individual systems are oscillators with a compact angle degree of freedom θi ∈ ( − π, π). Their uncoupled dynamics are determined by natural frequencies Pi (37), inertia parameters mi, and damping parameters di. Because the degrees of freedom are compact, the coupling between oscillators needs to be a periodic function of angle differences, and here, we only keep its first Fourier term. The coupling between pairs of oscillators is defined on a network whose Laplacian matrix has elements Lij(0)=−bij if i ≠ j and Lii(0)=∑k≠ibik. Without inertia, mi = 0 ∀i, Eq. 1 gives the celebrated Kuramoto model on a network with edge weights bij > 0, ∀i, j (18, 19). With inertia on certain nodes, it is an approximate model for the swing dynamics of high-voltage electric power grids in the lossless line approximation (16, 28, 38). The latter is justified in high-voltage transmission grids, where the resistance is smaller than the reactance typically by a factor of 10 or more. Applied to high-voltage grids, Eq. 1 describes the transient behavior of power grids on time scales of up to roughly 10 to 20 s. Over these time intervals, voltage amplitudes of high-voltage power grids are almost constant; accordingly, it is justified to consider only the dynamics of voltage angles (29). Here, we are interested in that transient time regime and accordingly focus on the voltage angle dynamics given by Eq. 1. When angle differences are small, a linear approximation sin(θi − θj) ≃ θi − θj is justified, giving first-order (without) or second-order (with inertia) consensus dynamics (17).

When the natural frequencies Pi are not too large, synchronous solutions exist that satisfy Eq. 1 with θ¨i=0 and θ̇i=ω0, ∀i. Without loss of generality, one may consider Eq. 1 in a frame rotating with the synchronous angular frequency ω0, in which case, these states correspond to stable fixed points with θ̇i=0. We consider a fixed point with angle coordinates θ(0)=(θ1(0),…,θn(0)) corresponding to natural frequencies P(0)=(P1(0),…,Pn(0)), to which we add a time-dependent disturbance, Pi(t)=Pi(0)+δPi(t). In the case of electric power grids, we will consider fixed points that are solutions to an optimal power flow problem. These solutions account for physical grid constraints such as thermal (i.e., capacity) limits of the lines and technical limitations of the power plants as well as economic constraints following from different production costs for different power plant types (see Materials and Methods and the Supplementary Materials) (39). Linearizing the dynamics about that solution, Eq. 1 becomesmi δθ¨i+di δθ̇i=δPi(t)−∑jbijcos(θi(0)−θj(0))(δθi−δθj),i=1,…,n(2)where δθi(t)=θi(t)−θi(0). This set of coupled differential equations governs the small signal response of the system corresponding to weak disturbances. The couplings are defined by a weighted Laplacian matrix Lij(θ(0))=−bijcos(θi(0)−θj(0)) if i ≠ j and Lii(θ(0))=∑kbikcos(θi(0)−θk(0)), which contains information about both the topology of the network and the operational state of the system. This weighted Laplacian matrix significantly differs from the network Laplacian L(0) when angle differences between coupled nodes are large.

We assess the nodal vulnerability of the system defined in Eq. 1 via the magnitude of the transient dynamics determined by Eq. 2 under a time-dependent disturbance δPi(t). We take the latter as an Ornstein-Uhlenbeck noise on the natural frequency of a single node, with vanishing average, δPi(t)¯=0, variance δP02, and correlation time τ0, δPi(t1)δPj(t2)¯=δik δjk δP02exp[−∣t1−t2∣/τ0]. It is sequentially applied on each of the k = 1, …, n nodes. This noisy test disturbance is designed to investigate network properties on different time scales by varying τ0 and identify the set of most vulnerable nodes, i.e., the key players, as the nodes where the system’s response to δPk(t) is largest. Besides being a probe to test nodal vulnerabilities, these noisy disturbances alternatively model fluctuating renewable energy sources in electric power grids. In this latter case, however, the correlation time τ0 is no longer a free parameter and is typically of the order of a minute or more, i.e., larger than any dynamical time scale in the system, as we discuss below. We quantify the magnitude of the response to the disturbance with the following two performance measures (40)𝒫1=limT→∞T−1∑i∫0T∣δθi(t)−Δ(t)∣2¯ dt(3A)𝒫2=limT→∞T−1∑i∫0T∣δθ̇i(t)−Δ̇(t)∣2¯ dt(3B)

They are similar to quadratic performance measures based on L2 or ℋ2 norms previously considered in the context of electric power grids, networks of coupled oscillators, or consensus algorithms (30, 40–46) but differ from them in two respects. First, here, we subtract the averages Δ(t) = n−1∑j δθj(t) and Δ̇(t)=n−1∑jδθ̇j(t), because the synchronous state does not change under a constant angle shift. Without that subtraction, artificially large performance measures may be obtained, which reflect a constant angle drift of the synchronous operational state but not a large transient excursion. Second, we divide 𝒫1,2 by T before taking T → ∞ because we consider a noisy disturbance that is not limited in time and that would otherwise lead to diverging values of 𝒫1,2.

Here, we calculate 𝒫1,2 for the network-coupled dynamical system defined in Eq. 1 when (i) both inertia and damping parameters are constant, mi ≡ m0 and di ≡ d0, (ii) the inertia vanishes, mi ≡ 0, (iii) the ratio γ ≡ di/mi is constant, and (iv) both inertia and damping vary independently. In cases (i) to (iii), 𝒫1,2 can be analytically expressed in terms of resistance centralities that will be introduced in the next section (see Materials and Methods and the Supplementary Materials). The next paragraphs focus on case (i), following which we present numerical data for case (iv), which illustrate the general applicability of these results for not too short noise correlation time.

The performance measures 𝒫1,2 can be computed analytically from Eq. 2 via Laplace transforms (see Materials and Methods and the Supplementary Materials), for homogeneous damping and inertia, i.e., di = d = γmi, ∀i. In the two limits of long and short noise correlation time τ0, they can be expressed in terms of the resistance centrality of the node k on which the noisy disturbance acts and of graph topological indices called generalized Kirchhoff indices (36, 40). Both quantities are based on the resistance distance, which gives the effective resistance Ωij between any two nodes i and j on a fictitious electrical network where each edge is a resistor of magnitude given by the inverse edge weight in the network defined by the weighted Laplacian matrix. One obtainsΩij(θ(0))=Lii†(θ(0))+Ljj†(θ(0))−Lij†(θ(0))−Lji†(θ(0))(4)where L† denotes the Moore-Penrose pseudo-inverse of L (36). The resistance centrality of the kth node is then defined as C1(k) = [n−1∑jΩjk]−1. It measures how central is the node kth in the electrical network, in terms of its average resistance distance to all other nodes—more central nodes have smaller C1(k). A network descriptor, the Kirchhoff index, is further defined as (36)Kf1≡∑i<jΩij(5)

Generalized Kirchhoff indices Kfp and resistance centralities Cp(k) can be defined analogously from the pth power of the weighted Laplacian matrix, which is also a Laplacian matrix (see Materials and Methods and the Supplementary Materials). In terms of these quantities, the performance measures defined in Eq. 3 depend on the value of the noise correlation time τ0 relative to the different time scales in the system. The latter are the ratios d/λα of the damping coefficient d with the nonzero eigenvalues λα, α = 2, …, n, of L(θ(0)) and the inverse ratio γ−1 = m/d of damping to inertia parameters. In high-voltage power grids, they are approximately given by d/λα < 1 s and m/d ≅ 2.5 s. Performance measures Eq. 3 can be obtained for any correlation time τ0 (see Materials and Methods and the Supplementary Materials). However, it is interesting to consider the specific cases where τ0 is the smallest (τ0 ≪ d/λα, γ−1) or the largest (τ0 ≫ d/λα, γ−1, appropriate for noisy power injections from new renewables) time scale in the probed system. The performance measures particularly take the asymptotic values𝒫1={(δP02τ0/d)(C1−1(k)−n−2Kf1), τ0≪d/λα,γ−1δP02(C2−1(k)−n−2Kf2), τ0≫d/λα,γ−1(6A)𝒫2={(δP02τ0/dm)(n−1)/n, τ0≪d/λα,γ−1(δP02/dτ0)(C1−1(k)−n−2Kf1), τ0≫d/λα,γ−1(6B)in the two limits when τ0 is the smallest or the largest time scale in the system. After averaging over the location k of the disturbed node, C1,2−1¯=2Kf1,2/n2, and one recovers the results in (40, 42, 43) for the global robustness of the system.

These results are remarkable: They show that the magnitude of the transient excursion under a local noisy disturbance is given by either of the generalized resistance centralities C1(k) or C2(k) of the perturbed node and the generalized Kirchhoff indices Kf1,2. The latter are global network descriptors and are therefore fixed in a given network with fixed operational state. One concludes that perturbing the less central nodes—those with largest inverse centralities C1,2−1(k)—generates the largest transient excursion. In a given network, key players are therefore nodes with smallest resistance centralities. It is important to keep in mind, however, that these centralities correspond to the weighted Laplacian defined above, where internodal couplings are normalized by the cosine of voltage angle differences. Accordingly, these centralities are dependent on the initial operating state. The asymptotic analytical results of Eq. 6 are corroborated by numerical results in the insets of Fig. 1, obtained directly from Eq. 1, i.e., without the linearization of Eq. 2. The validity of the general analytical expressions for any τ0 (see Materials and Methods and the Supplementary Materials) is further confirmed in the main panel of Fig. 1 and by further numerical results obtained for different networks shown in Materials and Methods and the Supplementary Materials.

The generalized resistance centralities and Kirchhoff indices appearing in Eq. 6 depend on the operational state via the weighted Laplacian L(θ(0)). For a narrow distribution of natural frequencies Pi ≪ ∑jbij, ∀i, angle differences between coupled nodes remain small, and the weighted Laplacian is close to the network Laplacian, L(θ(0))≃L(0). The resistance centralities C1(0) and C2(0) for the network Laplacian of the European electric power grid (see Materials and Methods and the Supplementary Materials) are shown in Fig. 2. For both centralities, the less central nodes are dominantly located in the Balkans and Spain. In addition, for C1(0), nodes in Denmark and Sicily are also among the most peripheral. The general pattern of these most peripheral nodes looks very similar to the pattern of most sensitive nodes numerically found in (47) and includes particularly many, but not all dead ends, which have been numerically found to undermine grid stability (48).

Fig. 2 Synchronous high-voltage power grid of continental Europe.

(A) Topology of the European electric power grid (see Materials and Methods and the Supplementary Materials) and location of the 10 test nodes listed in Table 1. Normalized generalized resistance centralities C1(0)(i) (B) and C2(0)(i) (C) for the network Laplacian matrix of the European electric power grid.

The asymptotic results of Eq. 6, together with the numerical results of Fig. 1, make a strong point that nodal sensitivity to fast or slowly decorrelating noise disturbances can be predicted by generalized resistance centralities. One may wonder at this point how generalized resistance centralities differ in that prediction from other, more common centralities such as geodesic centrality, nodal degree, or PageRank. Table 1 compares these centralities to each other and to the performance measures corresponding to slowly decorrelating noisy disturbances acting on the 10 nodes shown in Fig. 2A. As expected from Eq. 6, 𝒫1 and 𝒫2 are almost perfectly correlated with the inverse resistance centralities C2−1 and C1−1, respectively, but with no other centrality metrics. For the full set of nodes of the Europen electric power grid, we found Pearson correlation coefficients ρ(P1,C2−1)=0.997 and ρ(P2,C1−1)=0.975 fully corroborating the prediction of Eq. 6.

Table 1 Centrality metrics and performance measures 𝒫1,2 for the European electric power grid (see Materials and Methods and the Supplementary Materials) with noisy disturbances with large correlation time τ0 applied on the nodes shown in Fig. 2A.

The performance measures 𝒫1 and 𝒫2 are almost perfectly correlated with the inverse resistance centralities C2−1 and C1−1, respectively, but neither with the geodesic centrality, nor the degree, nor PageRank.

View this table:

DISCUSSION

Once a one-to-one relation between the generalized resistance centralities C1(k) and C2(k) of the disturbed node k and the magnitude of the induced transient response is established, ranking of nodes from most to least critical is tantamount to ranking them from smallest to largest C1 or C2. From Eq. 6, which of these two centralities is relevant depends on whether one is interested (i) in the transient response under fast or slowly decorrelating noise or (ii) in investigating transient behaviors for angles (using the performance measure 𝒫1) or frequencies (𝒫2). While this gives a priori four different rankings, Eq. 6 leads to only two rankings, based on either C1−1 or C2−1, which can be obtained through the performance measure 𝒫1 only, in asymptotic limit of either very fast (shortest time scale τ0) or very slowly (largest τ0) decorrelating noise. From here on, we therefore focus on the angle performance measure 𝒫1 of Eq. 3a and consider the two asymptotic limits in Eq. 6a.

We therefore define WLRank1 and WLRank2 (49) as two rankings, which order nodes from smallest to largest C1 and C2, respectively. Smallest WLRank1,2 therefore identify the most vulnerable nodes in a given network. Figure 3 shows that they differ very significantly. In particular, a number of nodes are among the most critical according to WLRank1 but not to WLRank2 and vice versa. This discrepancy means that nodes are not central in an absolute sense, instead, their centrality and hence how critical they are depend on details of the disturbance—in the present case, the correlation time τ0—and the performance measure of interest. One should therefore choose to use one or the other centrality measure, according to the network sensitivity one wants to check.

Fig. 3 Comparison of the two nodal rankings WLRank1 and WLRank2 obtained from the generalized resistance centralities C1 and C2, respectively, for the 3809 nodes of the European electric power grid sketched in Fig. 2A (see Materials and Methods and the Supplementary Materials).

Blue dots correspond to a moderate load during a standard winter weekday, and red dots correspond to a significantly heavier load corresponding to the exceptional November 2016 situation with a rather large consumption and 20 French nuclear reactors shut down.

The resistance centralities in Eq. 6 correspond to the network defined by the weighted Laplacian L(θ(0)) defined in Eq. 2. They therefore depend on the unperturbed, operating state θ(0), and consequently, WLRank depends not only on the network topology but also, as expected, on the natural frequencies and the coupling between the nodal degrees of freedom. As mentioned above, in the strong coupling limit, angle differences between coupled nodes remain small and L(θ(0))≃L(0). In that limit, one therefore expects nodal ranking to be given by resistance distances corresponding to the network Laplacian L(0). How long this remains true is of central interest, and to answer this question, we define further rankings LRank1,2 as the rankings using resistance centralities C1,2(0) obtained from the network Laplacian L(0). As long as angle differences between network-coupled nodes are not too large, the ranking LRank based on the network Laplacian matrix is almost the same as the true ranking WLRank based on the weighted Laplacian. This is shown in Fig. 4 for three electric power grid models and one random network of coupled oscillators. For the electric power grid models, injections/natural frequencies are limited by the standard operational constraint that the thermal limit of each power line is, at most, only weakly exceeded. This corresponds approximately to a maximal angle difference of max(Δθ) ≃ 30° between any pair of coupled nodes. Accordingly, we find that even in relatively strongly loaded power grids (corresponding, for instance, to the exceptional situation of the fall of 2016 when 20 French nuclear reactors were simultaneously offline; see red points in Fig. 4C), there is not much of a difference between LRank and WLRank. The two rankings start to differ from one another only when at least some natural frequencies become comparable with the corresponding nodal index, Pi ≲ ∑jbij, and angle differences become very large. This case has been investigated for an inertialess coupled oscillator system on a random rewired network with constant couplings (see Materials and Methods and the Supplementary Materials) (50). It is shown in green in Fig. 4D and corresponds to max(Δθ) = 106°.

Fig. 4 Comparison between LRank and WLRank corresponding to 𝒫1 for noisy disturbances with large correlation time τ0.

(A to C) Electric power grid models for normally (blue) and more heavily loaded (red) operating states governed by Eq. 1. (A) IEEE 57 bus test case where the more loaded case has injections six times larger than the moderately loaded, tabulated case (52). (B) MATPOWER Pegase 2869 test where the more loaded case has injections 30% larger than the moderately loaded, tabulated case (53). (C) European electric power grid model sketched in Fig. 2A (see Materials and Methods and the Supplementary Materials) where the moderately loaded case corresponds to a standard winter weekday and the more heavily loaded case to the November 2016 situation with 20 French nuclear reactors offline. For both cases, the operational state is obtained from an optimal power flow including physical, technological, and economic constraints (see Materials and Methods and the Supplementary Materials). (D) Inertialess coupled oscillators governed by Eq. 1 with mi = 0, ∀i, on a random network with 1000 nodes obtained by rewiring a cyclic graph with constant nearest and next-to-nearest neighbor coupling with a probability of 0.5 (see Materials and Methods and the Supplementary Materials) (50). Natural frequencies are randomly distributed as Pi ∈ [ − 1.8,1.63] (blue), Pi ∈ [ − 2.16,1.95] (red), and Pi ∈ [ − 2.7,2.45] (green), corresponding to maximal angle differences max(Δθ) = 31o, 70o, and 106o, respectively.

In Fig. 5, we investigate more closely when the approximate ranking LRank starts to differ from the true ranking WLRank. To that end, we used the randomly rewired model of inertialess coupled oscillators of Fig. 4D and calculated the percentage of nodes with highest LRank necessary to give the top 15% ranked nodes with WLRank2. The results are plotted as a function of the maximal angle difference between directly coupled nodes. Each of the 12,000 red crosses in Fig. 5 corresponds to one of 1000 natural frequency vectors P(0), with components randomly distributed in (−0.5,0.5) and summing to zero, multiplied by a prefactor β = 0.4,0.6,…2.4,2.6. The blue crosses correspond to running averages more than 500 red crosses with consecutive values of max(Δθ). One sees that, up to almost max(Δθ) ≃ 40°, the set of the 18% of nodes with highest LRank2 always includes the top 15% ranked nodes with WLRank2. Similar results for obtaining the top 10 and 20% ranked nodes with WLRank2 and for rankings using C1 instead of C2 are shown in Materials and Methods and the Supplementary Materials.

Fig. 5 Percentage of the nodes with highest LRank2 necessary to give the top 15% ranked nodes with WLRank2 for a random network of inertialess coupled oscillators with 1000 nodes obtained by rewiring with a probability of 0.5 of a cyclic network with constant nearest and next-to-nearest neighbor coupling (see Materials and Methods and the Supplementary Materials) (50).

Each of the 12,000 red crosses corresponds to one of 1000 random natural frequency vector P(0) with components randomly distributed in (−0.5,0.5) and summing to zero, multiplied by a prefactor β = 0.4,0.6,…,2.4,2.6. The blue crosses correspond to running averages more than 500 red crosses with consecutive values of max(Δθ). Inset: Running averages of the Frobenius distance between the matrices L(θ(0)) and L(0). The steps in the curve reflect discrete increments of β.

It is quite unexpected that nodal ranking remains almost the same up to angle differences of about 40°, since coupling nonlinearities are already well developed there. This is illustrated in the inset of Fig. 5, which plots the Frobenius distance ∑ij(Lij(θ(0))−Lij(0))2 between the network Laplacian L(0) and the weighted Laplacian L(θ(0)). When max(Δθ) ≃ 40°, the Frobenius distance has already reached about 27% of its maximal observed value, indicating that coupling nonlinearities are already strong. Yet, obtaining a desired set of the ns most critical nodes for any configuration with max(Δθ) ≲ 40°, including cases with nonegligible nonlinearities, is achieved with a single matrix inversion of the network Laplacian L(0), while considering a slightly extended set of ns + δns nodes with highest LRank, δns/ns ≪ 1. This is a moderate price to pay, compared to the price of calculating WLRank for each configuration, which each time requires inverting the weighted Laplacian matrix L(θ(0)). That latter procedure would be too time consuming for real-time assessment of large networks.

So far, we have assumed constant inertia and damping parameters, which led us to the analytical expressions given in Eq. 6 for the performance measures. Analytical results can further be obtained for inertialess systems with mi = 0 as well as in the case of homogeneous damping to inertia ratio, di/mi ≡ γ. In this latter case, the ranking is again given by a resistance centrality, but this time, related to the inertia-weighted matrix M−1/2LM−1/2 with M being the diagonal matrix whose ith diagonal entry is given by mi (see Materials and Methods and the Supplementary Materials) but not in the case of independently varying mi and di. We therefore lastly address this more general case using a purely numerical approach. This question is especially important for electric power grids where only nodes connected to rotating machines (such as conventional power plants) have inertia and consumer nodes have significantly smaller damping parameters (38). Time scales in electric power grids have typical values mi/di ∈ [1,3]s and di/λα ≲ 1s, and accordingly, we focus on the regime of large noise correlation time τ0 ≫ mi/di, di/λα, which is appropriate for persisting power fluctuations such as those arising from renewable energy sources. Figure 6 shows results corresponding to inertia and damping parameters fluctuating randomly from node to node by up to 40%. The ranking obtained from a full numerical calculation is compared to the ranking obtained from a direct calculation of the centrality of the weighted Laplacian Lij(θ(0)), corresponding to the long correlation time asymptotic limit of Eq. 6. One sees that the centrality-based ranking is close to the true, numerically obtained ranking, even in this case of strongly fluctuating inertia and damping parameters. This extends the validity of Eq. 6 for large τ0 in a much wider range of parameters than their derivation would suggest.

Fig. 6 Comparison between WLRank and numerical ranking for systems with inhomogeneous inertia and damping parameters.

(Left) Numerically obtained ranking based on the performance measure 𝒫1 plotted against the ranking WLRank2 based on the centrality C2 and (right) numerically obtained ranking based on the performance measure 𝒫2 plotted against the ranking WLRank1 based on the centrality C1. Each point is an average *over* more than 40 different noisy disturbances on a single node of the European electric power grid sketched in Fig. 2A, with independently fluctuating damping and inertia coefficients, di = d0 + δdi and mi = m0 + δmi with δmi/m0, δdi/d0 ∈ [−0.4,0.4] and γ = d0/m0 = 0.4 s−1. The noise correlation time is given by γτ0 = 4.

CONCLUSION

We have formulated a key player problem in deterministic, network-coupled dynamical systems. The formulation is based on the dynamical response to a nodal additive disturbance of the initial problem, and the most critical nodes—the key players—are defined as those where the response to the disturbance is largest. While this manuscript focused on (i) noisy Ornstein-Uhlenbeck disturbances, (ii) network-coupled systems on undirected graphs, particularly with symmetric couplings bij = bji in Eq. 1, and (iii) performance measures of the transient response that are quadratic forms in the system’s degrees of freedom, the method is not restricted to these cases. First, it can be used to deal with different disturbances, and in Materials and Methods and the Supplementary Materials, we calculate 𝒫1,2 for a box disturbance δPi(t) = δikδP0Θ(t)Θ(τ0 − t) with the Heaviside function Θ(t). This disturbance gives the same ranking as the Ornstein-Uhlenbeck noise disturbance considered above. Second, asymmetric couplings occurring, e.g., in directed graphs (51), in Kuramoto models with frustration (19), or in electric power grids with Ohmic dissipation (16) can also be considered. In this case, the internodal coupling is given by asymmetric real matrices instead of symmetric Laplacian matrices. However, the definition of the resistance distance (Eq. 4) remains valid even if L is replaced by an asymmetric matrix A, in that it still gives Ωii = 0, Ωij ≥ 0, and Ωij ≤ Ωik + Ωki, ∀i, j, k as long as the synchronous fixed point considered remains stable. Third, nonquadratic performance measures can, in principle, be considered within the spectral decomposition used in this article. One may think of average frequency nadir and rate of change of frequency, which are linear performance measures (30, 31). It is, at present, unclear whether these quantities can be analytically related to the location of disturbances via resistance or other centralities.

We gave an elegant answer to this key player problem: Ranking nodes from most to least critical is tantamount to ranking nodes from least to most central in the sense of resistance centralities. Depending on how the problem is formulated—mostly on details of the disturbance and on how the magnitude of the transient response is measured—different centralities have to be considered, giving different rankings. The key player problem in deterministic systems is therefore not uniquely defined, and its formulation must be tailored to reflect the most relevant dynamical properties one wants to evaluate. Averaged rankings, reflecting several such properties, simultaneously could also be considered. Last, we found numerically that resistance centralities are still accurate to identify the most critical nodes even when nodal dynamical parameters (damping and inertia) are not homogeneous.

The results shown in Fig. 6 are rather unexpected, and further inspection of our analytical results (Eq. 6 and eq. S14B) suggests that an inertia dependence could emerge in the opposite limit of short correlation time τ0 ≪ mi/di, di/λα. This point deserves further investigations. It would be furthermore interesting to extend our investigations to cases of distributions of inertia and damping parameters corresponding to realistic electric power grids. Work along those lines is in progress.

MATERIALS AND METHODS

Four different networks were considered. Data for implementing the IEEE 57 bus test case and the MATPOWER Pegase 2869 test case have been obtained from (52, 53). The random networks were obtained through the rewiring procedure of initially regular networks discussed in (50). We constructed the European power grid model from publicly available data on power plants and bus locations. More details on the procedure were given in the Supplementary Materials and in (54).

Several operational states of the European power grid were obtained from an optimal power flow. The latter minimizes a cost function including production costs for all available power plants, under constraints that power flows do not exceed the thermal line capacity of each power line and that productions do not exceed rated powers for each power plant.

For each network considered, Eq. 1 was numerically integrated using a fourth-order Runge-Kutta algorithm. Values for the performance measures were then calculated by numerical integration of the obtained time-dependent voltage angle and frequencies. Resistance distances, centralities, and Kirchhoff indices were calculated through a direct inversion of the Laplacian matrix.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/11/eaaw8359/DC1

Section S1. Calculation of the performance measures

Section S2. Resistance distances, centralities, and Kirchhoff indices

Section S3. Numerical models

Section S4. Numerical comparison of LRank with WLRank.

Fig. S1. Comparison between theoretical predictions and numerical results for both performance measures 𝒫1 and 𝒫2.

Fig. S2. Comparison of the performance measures 𝒫1, 𝒫2 obtained numerically and in eq. S14.

Fig. S3. Percentage of the nodes with highest LRank necessary to include the nodes with 10% and 20% highest WLRank.

References (55, 56)

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

REFERENCES AND NOTES

  1. ↵
    1. D. Kramer
    , Models poised to boost grid efficiency. Phys. Today 9, 25 (2016).
    OpenUrlCrossRefWeb of Science
  2. ↵
    National Academies of Sciences, Engineering, and Medicine, Analytic Research Foundations for the Next-Generation Electric Grid (The National Academies Press, 2016).
  3. ↵
    1. S. Fliscounakis,
    2. P. Panciatici,
    3. F. Capitanescu,
    4. L. Wehenkel
    , Contingency ranking with respect to overloads in very large power systems taking into account uncertainty, preventive, and corrective actions. IEEE Trans. Power Syst. 28, 4909–4917 (2013).
    OpenUrl
  4. ↵
    1. C. Ballester,
    2. A. Calvó-Armengol,
    3. Y. Zenou
    , Who’s who in networks. wanted: The key player. Econometrica 74, 1403–1417 (2006).
    OpenUrl
    1. S. P. Borgatti
    , Identifying sets of key players in a social network. Comput. Math. Organ. Theor. 12, 21–34 (2006).
    OpenUrl
    1. R. V. Solé,
    2. J. M. Montoya
    , Complexity and fragility in ecological networks. Proc. R. Soc. Lond. B Biol. Sci. 268, 2039 (2001).
    OpenUrlCrossRefPubMedWeb of Science
    1. J. M. Montoya,
    2. S. L. Pimm,
    3. R. V. Solé
    , Ecological networks and their fragility. Nature 442, 259 (2006).
    OpenUrlCrossRefPubMedWeb of Science
  8. ↵
    1. S. P. Borgatti
    , Centrality and network flow. Soc. Netw. 27, 55–71 (2005).
    OpenUrlCrossRefPubMedWeb of Science
  9. ↵
    1. P. Boldi,
    2. S. Vigna
    , Axioms for centrality. Internet Math. 10, 222 (2014).
    OpenUrlCrossRef
  10. ↵
    1. P. Bonacich
    , Power and centrality: A family of measures. Amer. J. Sociol. 92, 1170 (1987).
    OpenUrlCrossRef
  11. ↵
    1. S. Brin,
    2. L. Page
    , The anatomy of a large-scale hypertextual web search engine. Comput. Netw. ISDN Syst. 30, 107 (1998).
    OpenUrlCrossRef
  12. ↵
    1. M. Kitsak,
    2. L. K. Gallos,
    3. S. Havlin,
    4. F. Liljeros,
    5. L. Muchnik,
    6. H. E. Stanley,
    7. H. A. Makse
    , Identification of influential spreaders in complex networks. Nat. Phys. 6, 888 (2010).
    OpenUrl
  13. ↵
    1. J. Borge-Holthoefer,
    2. Y. Moreno
    , Absence of influential spreaders in rumor dynamics. Phys. Rev. E 85, 026116 (2012).
    OpenUrlAbstract/FREE Full Text
  14. ↵
    1. M. Girvan,
    2. M. E. J. Newman
    , Community structure in social and biological networks. Proc. Natl. Acad. Sci. U.S.A. 99, 7821 (2002).
    OpenUrl
  15. ↵
    1. M. K. Banda,
    2. M. Herty,
    3. A. Klar
    , Gas flow in pipeline networks. Netw. Heterog. Media 1, 41 (2006).
    OpenUrl
  16. ↵
    J. Machowski, J. W. Bialek, J. R. Bumby, Power System Dynamics: Stability and Control (Wiley, 2008), ed. 2.
  17. ↵
    N. A. Lynch, Distributed Algorithms (Morgan Kaufmann Publishers, 1997).
  18. ↵
    Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics (Springer, 1975), pp. 420–422.
  19. ↵
    1. J. A. Acebrón,
    2. L. L. Bonilla,
    3. C. J. Pérez Vicente,
    4. F. Ritort,
    5. R. Spigler
    , The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77, 137 (2005).
    OpenUrlCrossRefWeb of Science
  20. ↵
    1. E. Bompard,
    2. E. Pons,
    3. D. Wu
    , Analysis of the structural vulnerability of the interconnected power grid of continental Europe with the integrated power system and unified power system based on extended topological approach. Int. T. Electr. Energy 23, 620–627 (2013).
    OpenUrlCrossRefWeb of Science
  21. J. L. Sanchez Torres, “Vulnérabilité, Interdépendance et Analyse des Risques des Postes Sources et des Modes d’Exploitation Décentralisés des Réseaux Electriques,” PhD Thesis, University of Grenoble (2013).
  22. ↵
    1. K. N. Hasan,
    2. R. Preece,
    3. J. V. Milanović
    , Priority ranking of critical uncertainties affecting small-disturbance stability using sensitivity analysis techniques. IEEE Trans. Power Syst. 32, 2629 (2017).
    OpenUrl
  23. ↵
    1. S. V. Buldyrev,
    2. R. Parshani,
    3. G. Paul,
    4. H. E. Stanley,
    5. S. Havlin
    , Catastrophic cascade of failures in interdependent networks. Nature 464, 1025 (2010).
    OpenUrl
  24. ↵
    1. A. Basgan,
    2. Y. Berezin,
    3. S. V. Buldyrev,
    4. S. Havlin
    , The extreme vulnerability of interdependent spatially embedded networks. Nat. Phys. 9, 667–672 (2013).
    OpenUrlCrossRefPubMedWeb of Science
  25. ↵
    1. P. Hines,
    2. E. Cotilla-Sanchez,
    3. S. Blumsack
    , Do topological models provide good information about electricity infrastructure vulnerability? Chaos 20, 033122 (2010).
    OpenUrlCrossRef
  26. ↵
    1. M. Korkali,
    2. J. G. Veneman,
    3. B. F. Tivnan,
    4. P. D. H. Hines
    , Reducing cascading failure risk by increasing infrastructure network interdependency. Sci. Rep. 7, 44499 (2017).
    OpenUrlCrossRefPubMed
  27. ↵
    1. B. Bamieh,
    2. M. R. Jovanović,
    3. P. Mitra,
    4. S. Patterson
    , Coherence in large-scale networks: Dimension-dependent limitations of local feedback. IEEE Trans. Automat. Control 57, 2235 (2012).
    OpenUrl
  28. ↵
    1. F. Dörfler,
    2. M. Chertkov,
    3. F. Bullo
    , Synchronization in complex oscillator networks and smart grids. Proc. Natl. Acad. Sci. U.S.A. 110, 2005 (2013).
    OpenUrlCrossRefWeb of Science
  29. ↵
    P. Kundur, Power System Stability and Control (McGraw-Hill, 1994).
  30. ↵
    F. Paganini, E. Mallada, Global performance metrics for synchronization of heterogeneously rated power systems: The role of machine models and inertia, in Proceedings of the 55th Allerton Conference on Communication, Control, and Computing (IEEE, 2017), pp. 324–331.
  31. ↵
    L. Guo, G. Zhao, S. H. Low, Graph Laplacian spectrum and primary frequency regulation. arXiv:1803.03905 [cs.SY] (11 March 2018).
  32. ↵
    1. D. Manik,
    2. M. Rohden,
    3. H. Ronellenfitsch,
    4. X. Zhang,
    5. S. Hallerberg,
    6. D. Witthaut,
    7. M. Timme
    , Network susceptibilities: Theory and applications. Phys. Rev. E 95, 012319 (2017).
    OpenUrl
  33. ↵
    1. S. Tamrakar,
    2. M. Conrath,
    3. S. Kettemann
    , Propagation of disturbances in AC electricity grids. Sci. Rep. 8, 6459 (2018).
    OpenUrl
  34. ↵
    1. K. A. Stephenson,
    2. M. Zelen
    , Rethinking centrality: Methods and examples. Soc. Networks 11, 1 (1989).
    OpenUrlCrossRefWeb of Science
  35. ↵
    U. Brandes, D. Fleischer, Centrality measures based on current flow, Annual Symposium on Theoretical Aspects of Computer Science (Springer, 2005), pp. 533–544.
  36. ↵
    1. D. J. Klein,
    2. M. Randić
    , Resistance distance. J. Math. Chem. 12, 81 (1993).
    OpenUrlCrossRefWeb of Science
  37. ↵
    We allow ourselves a small abuse of language since, strictly speaking, the natural frequency of the ith oscillator would be Pi/di.
  38. ↵
    A. R. Bergen, D. J. Hill, A structure preserving model for power system stability analysis, in IEEE Transactions on Power Apparatus and Systems (IEEE, 1981), pp. 25–35.
  39. ↵
    A. R. Bergen, V. Vittal, Power Systems Analysis (Pearson/Prentice Hall, 2000), ed. 2.
  40. ↵
    1. M. Tyloo,
    2. T. Coletta,
    3. P. Jacquod
    , Robustness of synchrony in complex networks and generalized kirchhoff indices. Phys. Rev. Lett. 120, 084101 (2018).
    OpenUrl
    1. E. Tegling,
    2. B. Bamieh,
    3. D. F. Gayme
    , The price of synchrony: Evaluating the resistive losses in synchronizing power networks. IEEE Trans. Control Netw. Syst. 2, 254–266 (2015).
    OpenUrl
  42. ↵
    M. Siami, N. Motee, Systemic measures for performance and robustness of large-scale interconnected dynamical networks, in Proceedings of the 53rd IEEE Conference on Decision and Control (IEEE, 2014), pp. 5119–5124.
  43. ↵
    1. M. Siami,
    2. N. Motee
    , Fundamental limits and tradeoffs on disturbance propagation in linear dynamical networks. IEEE Trans. Automat. Contr. 61, 4055 (2016).
    OpenUrl
    1. B. K. Poolla,
    2. S. Bolognani,
    3. F. Dörfler
    , Optimal placement of virtual inertia in power grids. IEEE Trans. Automat. Contr. 62, 6209 (2017).
    OpenUrl
    1. T. Coletta,
    2. P. Jacquod
    , Performance measures in electric power networks under line contingencies. IEEE Trans. Control Netw. Syst. , (2019).
  46. ↵
    T. Coletta, B. Bamieh, P. Jacquod, Transient performance of electric power networks under colored noise, in IEEE Conference on Decision and Control (IEEE, 2018), pp. 6163–6167.
  47. ↵
    1. L. V. Gambuzza,
    2. A. Buscarino,
    3. L. Fortuna,
    4. M. Porfiri,
    5. M. Frasca
    , Analysis of dynamical robustness to noise in power grids. IEEE J. Emerg. Sel. Top. Circuits Syst. 7, 413 (2017).
    OpenUrl
  48. ↵
    1. P. J. Menck,
    2. J. Heitzig,
    3. J. Kurths,
    4. H. J. Schellnhuber
    , How dead ends undermine power grid stability. Nat. Commun. 5, 3969 (2014).
    OpenUrlPubMed
  49. ↵
    When used, subscripts LRanki and WLRanki indicate that these are rankings obtained from the centralities Ci, i = 1,2.
  50. ↵
    1. D. J. Watts,
    2. S. H. Strogatz
    , Collective dynamics of ‘small world’ networks. Nature 393, 440 (1998).
    OpenUrl
  51. ↵
    1. R. Delabays,
    2. P. Jacquod,
    3. F. Dörfler
    , The Kuramoto model on directed and signed graphs. SIAM J. Appl. Dyn. Syst. 18, 458–480 (2019).
    OpenUrlCrossRefPubMedWeb of Science
  52. ↵
    University of Washington, Electrical Engineering, Power Systems Test Case Archive (1993); www2.ee.washington.edu/research/pstca.
  53. ↵
    1. R. D. Zimmerman,
    2. C. E. Murillo-Sánchez,
    3. R. J. Thomas
    , MATPOWER: Steady-state operations, planning and analysis tools for power systems research and education. IEEE Trans. on Power Syst. 26, 12 (2011).
    OpenUrl
  54. ↵
    1. L. Pagnier,
    2. P. Jacquod
    , Inertia location and slow network modes determine disturbance propagation in large-scale power grids. PLOS ONE 14, e0213550 (2019).
    OpenUrl
  55. ↵
    B. Wiegmans, GridKit Extract of ENTSO-E Interactive Map (Zenodo, 2016); http://doi.org/10.5281/zenodo.55853.
  56. ↵
Acknowledgments: We thank R. Delabays and T. Coletta for interesting discussions. Funding: This work has been supported by the Swiss National Science Foundation under grants PYAPP2_154275 and 200020_182050. Author contributions: M.T. and L.P. wrote the numerical codes. L.P. constructed the model of the European power grid. M.T. performed the analytical and numerical calculations. P.J. formulated the ideas and wrote the initial draft. All authors participated in discussions of the results and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View this article with LENS

Article Tools

My saved folders

Stay Connected to Science Advances

Related Content

Similar Articles in:

Citing Articles in:

Navigate This Article