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.
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 form
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
When the natural frequencies Pi are not too large, synchronous solutions exist that satisfy Eq. 1 with
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,
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
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
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
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
The generalized resistance centralities and Kirchhoff indices appearing in Eq. 6 depend on the operational state via the weighted Laplacian
(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
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
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.
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
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.
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
(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.
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
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
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
(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
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.
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
- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).