Enhancing power grid synchronization and stability through time delayed feedback control

TL;DR

This paper investigates power grid synchronization and stability using a Kuramoto model with inertia, introducing time-delayed feedback control to enhance stability, demonstrated on a German grid model.

Contribution

It proposes a novel control method based on time-delayed feedback to improve power grid stability and synchronization, tested on a realistic grid model.

Findings

Time-delayed feedback control effectively stabilizes power grid synchronization.

Fewer controlled nodes are needed for stability with optimized control strategies.

The method is validated on a model of the German power transmission grid.

Abstract

We study the synchronization and stability of power grids within the Kuramoto phase oscillator model with inertia with a bimodal frequency distribution representing the generators and the loads. We identify critical nodes through solitary frequency deviations and Lyapunov vectors corresponding to unstable Lyapunov exponents. To cure dangerous deviations from synchronization we propose time-delayed feedback control, which is an efficient control concept in nonlinear dynamic systems. Different control strategies are tested and compared with respect to the minimum number of controlled nodes required to achieve synchronization and Lyapunov stability. As a proof of principle, this fast-acting control method is demonstrated using a model of the German power transmission grid.

Figures

Equations

$${\ddot{\theta }}_i+\alpha {\dot{\theta }}_i=\frac{P_i}{I_i{\omega }_{\text{G}}}+\frac{K}{I_i{\omega }_{\text{G}}}\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$${\ddot{\theta }}_i+\alpha {\dot{\theta }}_i$$

$$m_v{\dot{E}}_i$$

$$p(P)=\frac{1}{2\sigma \sqrt{2\pi }}\left(e^{-\frac{(P-P_0{)}^2}{2{\sigma }^2}}+e^{-\frac{(P+P_0{)}^2}{2{\sigma }^2}}\right).$$

$$P_{\text{Total}}=\sum _{i=1}^N{p}_i{P}_{\text{Total}}=36\text{ GW}$$

$$C_i=\{\begin{array}{ll}{\sum }_{k=1}^{n_i}{C}_i^k& n_i>0\\ 0& n_i=0\end{array}$$

$$C_{\text{Total}}=\sum _{i=1}^N{C}_i=88.354\text{GW}.$$

$$P_i=0.41\cdot C_i-p_i{P}_{\text{Total}}$$

$$\Delta \omega (t)\equiv \frac{1}{N}\sqrt{\sum _{i=1}^N{\left({\omega }_i(t)-\bar{\omega }(t)\right)}^2},$$

$$\delta {\ddot{\theta }}_i+\alpha \delta {\dot{\theta }}_i=\frac{K}{I{\omega }_{\text{G}}}\sum _{j=1}^N{A}_{ij}\cos ({\theta }_j-{\theta }_i)(\delta {\theta }_j-\delta {\theta }_i).$$

$$\delta {\ddot{\theta }}_i+\alpha \delta {\dot{\theta }}_i$$

$$m_v\delta {\dot{E}}_i$$

$${\ddot{\theta }}_i+\alpha {\dot{\theta }}_i$$

$$-\frac{g_i\alpha }{\tau }\left[{\theta }_i(t)-{\theta }_i(t-\tau )\right],$$

$$m{\ddot{\theta }}_i+{\dot{\theta }}_i={\Omega }_i+k\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i).$$

Full text

nonlinear complex networks, power grids, synchronization, stability analysis, time-delayed feedback control

I Introduction

Synchronization phenomena in nonlinear dynamical networks are of major interest to a wide field of applications in natural and technological systems Pikovsky et al. (2001); Boccaletti et al. (2018), e.g., neural networks in the human brain, or supply and communication networks and power grids, which naturally have a strong link to economy. Research in these fields has revealed diverse phenomena related to synchronization, ranging from partial synchronization patterns to asynchronous states van Vreeswijk (1996, 2000); Strogatz (2001). In particular, scenarios leading from full synchronization to asynchronicity via solitary states, i.e., single nodes which are desynchronized from the rest, play an important role for complex dynamical systems Maistrenko et al. (2014); Jaros et al. (2018), and in this work we will show that they are fundamental also for power grids.

Infrastructure, e.g., public transportation, medical care and a vast number of other everyday life applications, rely on electrical power supply. Given the fact that modern power transmission grids, notably if they include renewable energy sources, differ significantly from conventional power grids with regards to topology and local dynamics Milan et al. (2013); Heide et al. (2010, 2011), it is necessary to identify, understand, and cure the arising challenges and problems. In particular, malfunctioning grids can be the result of power outages, which occur for various reasons, including line overload or voltage collapse. Here we will focus on the loss of synchrony. In normal operation, a power grid runs in the synchronous state in which all frequencies equal the nominal frequency (50 or 60 Hz) and in which steady power flows balance supply and demand at all nodes. When parts of a power grid desynchronize, destructive power oscillations emerge. To avoid damage, affected components must then be switched off. However, such switchings can in turn desynchronize other grid components, possibly provoking a cascade of further shut-downs and ending in a large-scale blackout Union for the Coordination of Transmission of Electricity () (UCTE); Motter and Lai (2002); Buldyrev et al. (2010).

The failure of a transmission line during a blackout can be determined not only by the network topology and the static distribution of electric flow but also by the collective transient dynamics of the entire system where the time scale of system instabilities is of seconds Schäfer et al. (2018); Simonsen et al. (2008). In general, grids are designed such that the synchronous state is locally stable, implying that a cascade-triggering desynchronization cannot be caused by a small perturbation. However, even if the synchronous state is stable against small perturbations, the state space of power grids is also populated by numerous stable non-synchronous states to which the grid might be driven by short circuits, fluctuations in renewable energy generation or other large perturbations Chiang (2010); Anvari et al. (2016); Schäfer et al. (2018). Therefore it is of fundamental interest to explore the relation between network properties and grid stability against large perturbations Schäfer et al. (2015, 2018); Menck et al. (2013). Yet many intriguing questions on the relation between grid topology and local stability are still not understood. Decentralized grids tend to be less robust with respect to dynamical perturbations, but more robust against structural perturbations to the grid topology Rohden et al. (2012). However, adding new links may not only promote but also destroy synchrony, thus inducing power outages when geometric frustration occurs Witthaut and Timme (2012); Tchuisseu et al. (2018). The local stability can be improved by relating the specifics of the dynamical units and the network structure Motter et al. (2013); Dörfler et al. (2013); Menck et al. (2014), or predicting a priori which links are critical via the link’s redundant capacity and a renormalized response theory Witthaut et al. (2016).

In this paper we will demonstrate the role played by the solitary nodes in driving the populations out of synchrony and the necessity to control these nodes when restoring both stability and synchronization. Solitary nodes can be related to local instabilities via the application of a standard stability toolbox (i.e., Lyapunov exponents and Lyapunov vectors), and to topological properties of the network, like dead ends, thus complementing the analysis reported in Menck et al. (2014). Once we have identified the critical power grid nodes which undermine stability and synchronization, we will apply time-delayed feedback control to a small subset of these nodes, in order to cure a desynchronized and unstable power grid. Time-delayed feedback is an efficient mechanism known in nonlinear dynamics and often used to control unstable systems Pyragas (1992); Schöll and Schuster (2008). Generator and consumer dynamics will be described in terms of: i) Kuramoto oscillators with inertia Filatrella et al. (2008); ii) extended model of Kuramoto rotators with non trivial voltage dynamics or synchronous machines Schmietendorf et al. (2014). As a specific example, we consider the topology of the German ultra-high voltage power transmission grid (220 kV and 380 kV).

II Model and methods

II.1 The Kuramoto model with inertia

The Kuramoto model with inertia describes the phase and frequency dynamics of $N$ coupled synchronous machines, i.e., generators or consumers within the power grid, where mechanical and electrical phase and frequency are assumed to be identical:

[TABLE]

with the phase $\theta_{i}(t)$ and frequency $\dot{\theta}_{i}(t)=\frac{d\theta_{i}}{dt}$ of node $i=1,...,N$. Both dynamic variables $\theta_{i}(t)$, $\dot{\theta}_{i}(t)$ are defined relative to a frame rotating with the reference power line frequency $\omega_{\text{G}}$, e.g., 50 Hz for the European transmission grid. The distribution of net power generation ($P_{i}>0$) and consumption ($P_{i}<0$) is bimodal; it corresponds to the inherent frequency distribution in the Kuramoto model with rescaled parameters (see Appendix A for a detailed discussion on the parameter selection). The power balance requires $\sum_{i}P_{i}=0$.We assume homogeneously distributed transmission capacities $K$. The adjacency matrix $A_{ij}$ takes values 1 if node $i$ has a transmission line connected to node $j$, and 0 otherwise. Moreover $\alpha$ is the dissipation parameter and takes typical values of 0.1-1 s*-1* Menck et al. (2014); Machowski et al. (2008). Finally, the moment of inertia $I_{i}$ of turbine $i$ is $I_{i}=I=40\times 10^{3}\text{kg}\text{ m}^{2}$, corresponding to generation capacities of a single power plant equal to 400 MW Menck et al. (2014); Horowitz and Phadke (2008). With the above definitions, the frequency synchronization criterion reads $\dot{\theta}_{i}(t)=0,\quad\forall i=$1,…N, i.e., deviations from the reference frequency are zero.

II.2 Synchronous machine

Eq. (1) has been derived in Filatrella et al. (2008) from the swing equation governing the rotor’s mechanical dynamics Machowski et al. (2008), by assuming constant voltage amplitude and constant mechanical power $P_{i}$. The former assumptions make the model incapable of modeling voltage dynamics or the interplay of amplitude and phase. However it is possible to extend the model straightforwardly by including the voltage dynamics, thus taking into account the machine’s electrodynamical behavior. In the following we consider a lossless network of synchronous machines, whose dynamics is described by the extended model derived in Schmietendorf et al. (2014). The coupled dynamics of the phases $\left\{\theta_{i}\right\}$ and magnitudes $\left\{E_{i}\right\}$ of the complex nodal voltages $\left\{\textbf{E}_{i}=E_{i}e^{i\theta_{i}}\right\}_{i\in\{1,...,N\}}$ are given by

[TABLE]

where $\dot{\theta}_{i}$ is the individual frequency of the $i-$th oscillator. $P_{i}$ denotes the mechanical input or output power and $KA_{ij}E_{i}E_{j}\sin(\theta_{j}-\theta_{i})$ is the electrical real power transferred between machines $i$ and $j$. The susceptance matrix coefficients $A_{ij}$ allow for variations concerning the network topology; as for the previous model, $A_{ij}=1$ if node $i$ has a transmission line connected to node $j$, 0 otherwise. In particular the diagonal entries of $A_{ij}$ are chosen such that the matrix has zero row sum $\sum_{i=1}^{N}A_{ij}=0$. $m_{v}$, $E_{f,i}$, $X_{i}$ take into account machine and line parameters. In particular these parameters are set to be homogeneous and of the same order of magnitude as in Schmietendorf et al. (2014): $m_{v}=1$, $E_{f,i}=1$, $X_{i}=1$, while the remaining quantities, already discussed in the original model Eq. (1), are chosen as $\alpha=5/6$ s*-1*, $I_{i}=I=40\cdot 10^{3}\text{kg}\text{ m}^{2}$, $\omega_{\text{G}}=2\pi\cdot 50~{}$ Hz.

II.3 German power grid and power distributions

In our numerical example we extract the topology $A_{ij}$ from the Open Source Electricity Model for Germany (elmod-de) Egerer (2016), which describes the German ultra-high voltage transmission grid using $N=438$ nodes connected by 662 transmission lines (see Figure 1a).

In many previous studies using the Kuramoto model with inertia to model power grid networks, the distribution of net power generation and consumption $P_{i}$ is set to be a bimodal $\delta$-distribution Rohden et al. (2012, 2014); Witthaut and Timme (2012); Lozano and Buzna (2012); Menck et al. (2014); Olmi et al. (2014). Here we consider more complex distributions: first of all, an artificial bimodal Gaussian distribution $P^{\mathrm{G}}$ Olmi and Torcini (2016); Tumash et al. (2018) is generated, whose probability density function $p(P)$ is given by the superposition of two Gaussians centered at $\pm P_{0}$ with standard deviation $\sigma$

[TABLE]

Figure 1b shows a histogram of the realization $P^{\mathrm{G}}$ used in the numerical simulations of this study. The second distribution $P^{\mathrm{R}}$ shown in Figure 1c is calculated based on data provided by elmod-de Egerer (2016) and will be referred to as real-world distribution.

According to the data documentation Egerer (2016), elmod-de is an open source nodal DC load flow model, minimizing generation costs, for the German electric power transmission grid. In the following, we point out how the information in elmod-de is translated into realistic values for the parameters used in our network of Kuramoto oscillators with inertia. As anticipated above, the data set contains nodal information on $N=438$ network nodes within the 220 kV and 380 kV ultra-high voltage transmission grid, of which 393 are substations. The remaining nodes are used to model interactions with neighboring countries (22) and auxiliary nodes (23), e.g., points in the grid without a transformer station. The nodes are connected with 697 transmission lines, 35 of them appearing twice in the data set, which will be neglected, such that 662 unique transmission lines remain. We will furthermore assume identical power transmission capacities for all transmission lines, resulting in a generic coupling strength for the network, thus reducing the values of the coupling matrix to 0 or 1. Besides geographical locations of all nodes, local power demand values $p_{i}$ are given in parts of the total power demand of Germany at off-peak times:

[TABLE]

Following these definitions the absolute power demand at node $i=1,...,N=438$ is given by $p_{i}P_{\text{Total}}$. The spatial distribution of $p_{i}P_{\text{Total}}$ is illustrated in Fig. 2.

Furthermore 562 conventional power plants, e.g., coal or atomic plants, are listed. Information on the topological location of plants, i.e., to which node $i$ they belong, and their maximum power generation capacities is provided. Let $n_{i}$ be the number of power plants associated with node $i$. The maximum capacity of plant $k=1,...,n_{i}$ located at node $i$ will be denoted by $C_{i}^{k}$. In order to obtain node-wise generation capacities $C_{i}$, $C_{i}^{k}$ will be aggregated for each node:

[TABLE]

The spatial distribution of $C_{i}$ is illustrated in Fig. 3. The total generation capacity $C_{\text{Total}}$ reads:

[TABLE]

Due to the fact that plants being operated at 100% of their maximum generation capacities would cause a large oversupply of power generation and break power balance, we will assume each plant to be operated at 41% of its maximum capacity, since $P_{\text{Total}}/C_{\text{Total}}\approx 0.41$. With this intermediate level of power generation, the power balance is fulfilled and the net generation/consumption $P_{i}$ at node $i$ is given by:

[TABLE]

II.4 Macroscopic indicators and Lyapunov analysis

We consider a scenario where, due to an arbitrary dynamical perturbation, some critical nodes have become desynchronized, where we define as critical those nodes withstanding self-organized resynchronization. Synchronization is first gained by performing an adiabatic transition from the asynchronous to the synchronized state for increasing coupling constant: starting with random initial conditions $\theta_{i}(0)\in[-2\pi,2\pi)$, $\dot{\theta}_{i}(0)\in[-1,1)$ at $K=0$, the coupling strength $K$ is increased adiabatically up to $K_{\mathrm{Max}}$ where the system shows synchronized behavior. For each investigated value of $K$, the system is initialized with the final conditions found for the previous coupling value, then the system evolves for a transient time $T_{A}$, such that it can reach a steady state. After the transient time $T_{A}$, characteristic measures are calculated in order to assess the quality of synchronization and the stability of the underlying state $\{\theta_{i}(T_{A}),\dot{\theta}_{i}(T_{A})\}$. In particular the time-averaged phase velocity profile $\left<\omega_{i}\right>_{t}\equiv\left<\dot{\theta}_{i}\right>_{t}$ provides information on frequency synchronization of individual nodes $i$, whereas the standard deviation of frequencies

[TABLE]

is used to estimate the deviation from complete frequency synchronization ($\bar{\omega}(t)$ indicates the instantaneous average grid frequency).

Once a desired synchronized state is reached, a perturbation can occur leading the state out of synchrony. In this situation the overall stability of the power grid might be lost, therefore it is necessary to analyze the time-evolution of small dynamic perturbations $\delta\theta_{i}:=\theta_{i}^{*}-\theta_{i}$ around the steady state $\theta_{i}^{*}$, whose dynamics is ruled by the linearization of Eq. (1) as follows

[TABLE]

For the extended model, the linearization of Eqs. (2, 3) reads as

[TABLE]

The exponential growth rates of the infinitesimal perturbations are measured in term of the associated Lyapunov spectrum $\{\lambda_{k}\}$, with $k=1,...,2N$, numerically estimated by employing the method developed by Benettin et al. Benettin et al. (1980). In particular one should consider for each Lyapunov exponent $\lambda_{k}$ the corresponding 2N-dimensional tangent vector $\mathcal{T}^{(k)}=(\delta\dot{\theta}_{1},...,\delta\dot{\theta}_{N},\delta\theta_{1},...\delta\theta_{N})$ whose time evolution is given by Eq. (10) (resp. Eqs. (11, 12) for the extended model). Important information about the sources of instability and, in particular, about the oscillators that are more actively contributing to the chaotic dynamics, can be gained by calculating the time averaged evolution of the tangent vector $\mathcal{T}^{(1)}$, here referred to as maximum Lyapunov vector. The Euclidean norm of each $\{\theta_{i},\dot{\theta}_{i}\}$ pair in $\mathcal{T}^{(1)}$, averaged in time, is measured for each oscillator as $\xi_{i}:=\left<\sqrt{[\delta\theta_{i}(t)]^{2}+[\delta\dot{\theta}_{i}(t)]^{2}}\right>_{t}$, once the tangent vector is orthonormalized, i.e. $||\mathcal{T}^{(1)}||=1$.

III Results for a network of Kuramoto oscillators with inertia

III.1 Emergence of solitary states

In general we have performed sequences of simulations by varying adiabatically the coupling parameter $K$ with two different protocols. Namely, for the upsweep protocol, as described in the previous section, the series of simulations is initialized for the decoupled system by considering random initial conditions both for phases and frequencies. Afterwards the coupling is increased in steps of $\Delta K$ until a maximum coupling strength $K_{\mathrm{Max}}$ is reached. For the downsweep protocol, starting from the maximum coupling strength $K_{\mathrm{Max}}$ achieved by employing the upsweep protocol simulation, the coupling is reduced in steps of $\Delta K$ until $K=0$ is recovered. At each step the system is simulated for a transient time $T_{A}$ followed by a time interval $T_{W}$ during which the average frequencies $\left\langle\omega_{i}\right\rangle_{t}$, as well as the components of the Lyapunov vector $\xi_{i}$ and the maximum Lyapunov exponent $\lambda_{1}$ are calculated. An example of the results obtained by performing the sequence of simulations of upsweep followed by downsweep is shown in Figs. 4 and 5 for the bimodal Gaussian distribution $P^{\mathrm{G}}$ and the real-world distribution $P^{\mathrm{R}}$, respectively.

In both cases, at low coupling, a large fraction of the network is unsynchronized (panel (a)) and the system is chaotic, i.e., $\lambda_{1}>0$ (b). A considerable part of the oscillators rotates with average frequency $\left|\left<\omega_{i}\right>\right|>0$, while relatively few oscillators are locked at average zero frequency (c). Other clusters at $\left|\left<\omega_{i}\right>\right|\neq 0$ may emerge. The solitary nodes, which are desynchronized from the rest of the network, and oscillate with high frequency, are those mostly responsible for the lack of synchronization. This is revealed by the analysis of the components of the maximum Lyapunov vector $\xi_{i}$, which assume large values for those nodes which are solitary, thus indicating that the directions identified by solitary nodes are the most unstable in the network (as shown in panel (f)).

For intermediate $K$ values, the majority of nodes is synchronized on average, with a small set of nodes being solitary, for instance, 9 for the Gaussian and 11 for the real-world distribution (see panels (a) and (d)). The system is still chaotic (panel (b)) and the components of the Lyapunov vector $\xi_{i}$ are still localized around solitary nodes (panel (g)). The number of solitary nodes diminishes for increasing coupling values, since more and more nodes join the main synchronized cluster at zero average frequency. Just before full synchronization (see panel (e)), one solitary node is left and no instability emerges in the system (panel (h)). The full synchronized state ($\left\langle\Delta\omega\right\rangle_{t}=0$) is stable and it is characterized by a single cluster with no solitary nodes. In particular complete frequency synchronization with $\left<\Delta\omega\right>_{t}=0,\lambda_{1}=0$ is achieved at $K\geq 1320$ MW ($K\geq 4200$ MW) for $P^{\mathrm{G}}$ ($P^{\mathrm{R}}$).

When $K$ is decreased starting from the synchronized states, the systems remains synchronized for a larger $K$ interval, due to the hysteretic nature of the transition, and the synchronized state loses stability (i.e., $\left\langle\Delta\omega\right\rangle_{t}>0$) for a coupling value smaller than the one found during the upsweep protocol (see panel (a)). The system is multistable and partially synchronized states (as those shown in panel (d)) coexist with the synchronized one. Depending on the initial state of the system, the dynamics can approach either the synchronized state or one of the upper branch states. This also means, that, starting from the synchronized states, large perturbations can kick the system out of synchrony. The goal of this paper it to give a proof of principle that once such a partially synchronized state is approached, our control method is capable of synchronizing and stabilizing the system. Thus in the following we consider the unstable states present in panels (d), (g) of Figures 4 and 5, which we aim to control.

III.2 Application of time-delayed feedback control

To facilitate understanding we report in Fig. 6 the main features of the unstable states, briefly introduced in the previous section, that we aim to control. In particular Fig. 6 shows the time-averaged standard frequency deviation $\left<\Delta\omega\right>_{t}$ and the maximum Lyapunov exponent $\lambda_{1}$ for each value $K$ of the adiabatic increase for the bimodal Gaussian (panel a) and for the real-world distribution (panel b) and highlights the considered operating points via dashed black lines.

If a perturbation pushes the system out of synchrony at an intermediate state with finite values of $\left<\Delta\omega\right>_{t}$, in a chaotic regime characterized by $\lambda_{1}>0$, would it be possible to enhance synchronization and stability by controlling a small subset of nodes? In the following we will give a positive answer to this question, by exploring the dynamics of the system at $K\approx 816$ MW ($K\approx 729$ MW) for $P^{\mathrm{G}}$ ($P^{\mathrm{R}}$), where deterministic chaos is present, i.e., $\lambda_{1}=0.0187$ ($\lambda_{1}=0.096$), and the system is not perfectly frequency synchronized: $\left<\Delta\omega\right>_{t}\approx 0.34$ Hz ($\left<\Delta\omega\right>_{t}\approx 0.91$ Hz), modeling a strongly perturbed power grid pert . Even though we are considering a partially synchronized regime with an intermediate transmission capacity value, as a resulting regime in case of strongly perturbed grid, we made sure not to artificially drive the system to an unrealistic range of capacity values. Indeed the operating point at which we are working is in a realistic regime when considering the average transmission capacity ($\approx 1500$ MW) at which the German ultra-high voltage transmission grid works, according to elmod-de data set.

From the average frequency profile shown in Fig. 6, panel c (panel d) for $P^{\mathrm{G}}$ ($P^{\mathrm{R}}$), we can see that a major part of the power grid is frequency synchronized while few nodes have a significant frequency deviation and are identified as solitary states: 9 nodes for $P^{\mathrm{G}}$, 11 nodes for $P^{\mathrm{R}}$. (Note that the three solitary nodes $i=1,2,3$ can only be resolved in the blown-up inset.) Solitary nodes oscillate with their own average frequency and do not resynchronize in a self-organized way at a given coupling strength, being thus critical for desynchronization. Note that the solitary nodes include those with the largest $\xi_{i}$, but not only those.

In order to enhance frequency synchronization and stability at the intermediate coupling strength discussed above, the Kuramoto model with inertia is now extended by time-delayed feedback control which is an efficient control concept, well known in nonlinear dynamic systems Pyragas (1992); Schöll and Schuster (2008), but also commonly employed in power grid engineering Machowski et al. (2008); Kundur et al. (1994):

[TABLE]

where $g_{i}$ is the control gain of node $i$ and $\tau$ is the delay time. The control method turns out to be robust against changes in the parameters $\tau,g_{i}$. We propose in the following to apply the control term only to a small subset of nodes selected according to their dynamical properties.

In order to find such a set, different control strategies are proposed in the following: (i) the first strategy takes into consideration all solitary nodes, sorted in descending order of $\xi_{i}$; (ii) the second strategy orders the solitary nodes by their absolute average frequency $\left|\left<\omega\right>_{t}\right|$; (iii) the third strategy consider all nodes, not only solitary ones, randomly picked. The outcome of the different strategies is shown in Fig. 7(a)-(c) and (d)-(f) for the bimodal Gaussian distribution $P^{\mathrm{G}}$ and the real-world distribution $P^{\mathrm{R}}$, respectively. First of all, strategy (i) is able to achieve stability if just one node is controlled, and frequency synchronization if the number of controlled solitary nodes is sufficiently large: 8 controlled nodes for both $P^{\mathrm{G}}$ and $P^{\mathrm{R}}$. Strategy (ii) requires 4 controlled nodes for stabilization and 8 for synchronization in case of $P^{\mathrm{G}}$, and one controlled node for stabilization and 9 nodes for synchronization in case of $P^{\mathrm{R}}$. The third strategy is not able to frequency-synchronize and stabilize, it can at most mitigate to some extent the desynchronization and the instability. For the given setup, strategy (i) is the best choice: it is particularly efficient since the Lyapunov vector is re-calculated every time when an additional solitary node is controlled, thus taking into account the interplay between solitary states and emerging instabilities. However, both strategies (i) and (ii) highlight the role played by solitary nodes, role that will be clarified in more details in the next section.

III.3 Lyapunov analysis

The presence of solitary nodes deeply influences the dynamics emerging in the system, since they behave almost independently, adding complexity and conveying the instability. In particular the role played by the solitary nodes can be understood by the change in the Lyapunov spectrum when the control strategy (i) is applied, i.e., when solitary nodes are controlled, ordered according to their Lyapunov vector component (for the definition of the other strategies see previous section).

If we first consider the bimodal Gaussian frequency distribution, the uncontrolled state is characterized by a cluster of synchronized oscillators plus 9 solitary nodes. The system is chaotic and the maximum Lyapunov exponent is positive (see Fig. 8a): the interplay between solitary nodes and cluster state gives rise to low-dimensional chaos in the system. When the first solitary node is controlled (Fig. 8b), the dynamics becomes quasiperiodic and the collective behavior is a high-dimensional torus, as can be deduced by the consistent number of 8 Lyapunov exponents that are exactly zero. Each solitary node, at the microscopic level, moves with an average velocity which is different from the velocity of the cluster and from the velocity of the other solitary states: the self-emergent dynamics, at the macroscopic level, is a quasiperiodic motion characterized by multiple incommensurable frequencies. When solitary nodes are controlled and frequency synchronized to the cluster, they do no longer contribute to the collective dynamics with their own frequency, thus decreasing the dimensionality of the macroscopic behavior. Thus, the further control of more solitary nodes has the effect of stabilizing the system: negative exponents becomes more and more negative while the zero ones become negative. When 5 solitary states are controlled, the macroscopic dynamics evolves on a 2-dimensional torus (see Fig. 8f). This can be explained considering that in the system under investigation one might expect two Lyapunov exponents to be zero due to the symmetries of the system: one is always present for a system with continuous time, while the second zero exponent is related to the invariance of the model under uniform phase shift. Therefore when 5 solitary nodes are controlled, 2 exponents are zero due to symmetries, while the other 2 zero exponents identify the emergent quasiperiodicity. Finally, when the system is synchronized, thanks to the control of 8 solitary nodes, the typical spectrum of a stable periodic synchronized state appears, with a negative plateau at $\lambda_{n}=-\alpha/2$ (for $1<n<2N-1$) and $\lambda_{1}=0$ (see Fig. 8i, j). The synchronized state is degenerate and the phase shift of all the phases corresponds to a perturbation along the orbit of the fully synchronized state, which explains why the two invariances, and thus the Lyapunov exponents, coincide, as already shown in Olmi S. (2015) for a globally coupled network.

A similar behavior can be observed for the real-world frequency distribution case, where the initial uncontrolled state is chaotic ($\lambda_{1}>0$) and 11 solitary nodes emerge from the synchronized cluster state (see Fig. 9 a). When the solitary state with largest Lyapunov component is controlled and synchronized to the cluster, the system is no longer unstable, which indicates that the instability was conveyed by the selected solitary node (see Fig. 9b). Due to the interaction of the remaining solitary states, characterized by different average frequencies, the collective dynamics of the system turns out to be quasiperiodic and high-dimensional. The dimensionality of the quasiperiodic motion is reduced by controlling more and more nodes and results in a 2-dimensional torus when 5 solitary nodes are controlled (see Fig. 9f). Finally the system is synchronized when 8 solitary states are controlled (see Fig. 9i), while the additional control of further nodes does not alter nor enhance the synchronization.

III.4 Topological features vs Extreme events

In Menck et al. (2014) numerical evidence was given that dead ends and dead trees undermine basin stability of nodes in Kuramoto power grid networks, which means that the basin of attraction of the frequency synchronized solution for single nodes tends to be small if a node is placed at a dead end, thus making such nodes hard to synchronize. Indeed, in the case of the bimodal Gaussian distribution $P^{\mathrm{G}}$, all the identified solitary nodes belong to a dead tree (see Fig. 10a). However, this trend cannot be observed for the real-world distribution $P^{\mathrm{R}}$, where just 3 of the 11 solitary nodes belong to a dead tree (see Fig. 10b) and dead trees do not correspond to the most unstable nodes. In general we have observed that the most unstable solitary nodes, for $P^{\mathrm{G}}$, are dead ends adjacent to well connected nodes, whereas for $P^{\mathrm{R}}$ they are nodes with $P_{i}>4\Delta P$, where $\Delta P$ is the standard deviation of the distribution. The discrepancy between the two cases can be explained if, starting from $P^{\mathrm{G}}$, we arbitrarily add $4\Delta P$ to the net power ($\widehat{=}$ inherent frequency) of a non-solitary node $k$. This altered node then becomes solitary and causes other adjacent nodes to become solitary, some of them belonging to dead trees. If we control all the newly emerged solitary dead trees, the system does not synchronize and the dynamics of node $k$ is almost unchanged (Fig. 10c), whereas we can achieve synchronization via controlling node $k$ only (Fig. 10d). This means that dead trees are fundamental in determining the power grid stability whenever the power distribution does not contain fat tails or extreme events, which is the case for $P^{\mathrm{G}}$; for the real-world distribution $P^{\mathrm{R}}$, however, nodes with significant power difference are common and the stability is undermined by these nodes rather than by dead trees.

IV Results for a lossless network of synchronous machines

Applying the same procedure as previously done for the standard Kuramoto model with inertia with different frequency distributions, we perform an adiabatic parameter scan in $K$, thus identifying the synchronization transition of the system during the upsweep and downsweep protocols. The system is initialized at $K=0$ with uniformely distributed initial conditions not only for phases and frequencies $\{\theta_{i},\dot{\theta}_{i}\}$, but also for the voltage amplitudes $\{E_{i}\}$, that are set uniformly random: $E_{i}(0)\in[0.5,1.5)$.

As for the previously investigated setups, the system undergoes a hysteretic transition to synchronization (see Fig. 11a). It shows an asynchronous state for low coupling values $K$, and partially synchronized states for intermediate $K$ values (panels d, e). In particular the number of whirling nodes diminishes with increasing $K$ and it is possible to identify a state, in proximity of the synchronization transition, where almost all nodes are synchronized, while few of them are solitary nodes still oscillating with average frequency different from zero (panel e). Similarly to the previous setups, the Lyapunov vector is (mostly) localized around solitary nodes (see Fig. 11, panels (c)-(h) corresponding to different stages of the adiabatic upsweep), thus indicating that solitary nodes are leading the synchronization transition even when considering voltage dynamics. Finally, the system is chaotic for a larger $K$ interval (see panel b) as compared to the original Kuramoto model with inertia.

Strategies (i) and (ii) to synchronize and stabilize the system are applied to the partially synchronized state at $K\approx 1307$ MW (see Figure 11 panels e, h), where 13 solitary nodes are present: a comparison of the strategies is shown in Fig. 12.

The first strategy requires to control one node in order to stabilize the system and 11 to synchronize, whereas the second strategy performs worse when stabilizing the system (4 nodes required) but performs better when synchronizing (10 nodes). However both control schemes require not all solitary nodes to be controlled in order to achieve synchronization and stability. All in all our approach is not only applicable to the example systems presented in Sec. 3, but works for different models. In Appendix B, the generality of the approach will be further explained considering different topologies and different operating points. Even though it is not possible to provide an analytical proof of the efficiency and generality of our control approach, our results indicate how powerful and robust time-delayed feedback control is, and that it can be applied to a diversity of topologies and power grid models. The hysteretic nature of the transition to synchronization, the bistability of the system, and the emergence of solitary states driving the dynamics, are fundamental ingredients for enhancing the stability of power grids, which have not been recognized until now.

V Conclusions

In conclusion, we have proposed a time-delayed feedback control scheme to restore frequency synchronization and stability of the power grid after perturbations. To this purpose we have firstly studied the Kuramoto model with inertia in the presence of two different bimodal distributions of generator and load power (an artificial distribution, and one adapted from the real German high-voltage transmission grid), which both lead to a fully frequency synchronized, stable network for large transmission capacities $K$. We have focussed on the operating regime of intermediate $K$ characterized by a number of solitary nodes whose mean frequency deviates from that of all other nodes.

We have shown that stability and synchronization can be enhanced by time-delayed feedback control in this $K$ regime by applying delayed feedback to a small subset of nodes: frequency synchronization and stability can be restored in a short time and persist even if control is turned off. Different control strategies were tested. For the shown setup the best strategy is to control the most unstable solitary nodes, characterized by the largest Lyapunov vector components. However, both strategies (i) and (ii) are efficient, being based on the solitary nodes that turn out to be fundamental in regulating the dynamics of the system. Solitary nodes exhibit independent dynamics, giving rise to low-dimensional chaos that turns into high-dimensional quasi-periodic motion when the most unstable node is controlled, until synchronization is achieved. Therefore, due to their independence, the set of controlled nodes cannot be much smaller than the number of solitary nodes.

The proposed fast-acting control method might offer an interesting approach to cure disturbances in real-world power grids, due to its general applicability and validity, as shown in Sec. IV, where we have applied our control strategy to a more sophisticated model including the voltage dynamics Schmietendorf et al. (2014) and, more in general, as shown in Appendix B, where we have extended our analysis to a different network (i.e., the Italian grid) and to different operating points, keeping the German grid topology.

Acknowledgements.

We acknowledge A. Torcini and S. Lepri for valuable discussions. Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 163436311 - SFB 910.

Appendix A: Parameter choice

As already detailed in Sec. II A, the Kuramoto model with inertia describes the phase and frequency dynamics of $N$ coupled synchronous machines, i.e., generators or consumers within the power grid, where mechanical and electrical phase and frequency are assumed to be identical. The $N$ dynamic equations describing the time evolution of the phase $\theta_{i}(t)$ and frequency $\dot{\theta}_{i}(t)=\frac{d\theta_{i}}{dt}$ of node $i=1,...,N$ are given by Eq. 1. In particular $\alpha$ represents the dissipation parameter and takes typical values of 0.1-1 s*-1* Menck et al. (2014); Machowski et al. (2008). However, in a realistic power grid there are additional sources of dissipation, especially Ohmic losses, and losses caused by damper windings Machowski et al. (2008), which are not taken into account directly in the coupled oscillator model. Therefore, for this parameter we have chosen slightly higher values: $\alpha=5/6$ s*-1* when a bimodal Gaussian distribution is considered and $\alpha=2$ s*-1* when the real-world distribution is taken into account to describe the distribution of the net power $P_{i}$. Different dissipation values are necessary for the different distributions in order to obtain comparable setups , i.e., unstable, partially synchronized states at comparable coupling strengths, $K=819$ MW for the bimodal Gaussian distribution and $K=729$ MW for the real-world one.

For both net power distributions, the coupling strength $K$, which represents the maximum power transmission capacity of transmission lines, was set homogeneously throughout the grid. A more realistic approach would have been to use a coupling matrix $K_{ij}$, containing not only the topology, but also individual transmission capacities to schematize different transmission line lengths. However the goal of the present paper is to gain insight into the principal behavior of large power grids depending on the network topology, and their capability to synchronize by controlling a minimal set of nodes and, for a proof of principle of our control approach, the choice of identical transmission lines suffices. The choice of using simplified homogeneous transmission line capacities (coupling constants) turned out to be a good compromise when using heterogeneous power distributions, whose realistic values were available in the open data source as opposed to the power distribution data.

Eq. (1) can be simplified by rescaling the parameters $m:=\frac{1}{\alpha}$, $\Omega_{i}:=\frac{P_{i}}{I_{i}\omega_{\text{G}}\alpha}$, $k:=\frac{K}{I_{i}\omega_{\text{G}}\alpha}$, thus giving

[TABLE]

In comparison with Eq.(1), the inertial mass $m$ now represents the inverse of the dissipation $\alpha$ in the grid, and the coupling constant $k$ now represents the maximum power which can be transmitted between two connected nodes. Moreover each node $i$, when uncoupled, oscillates with an angular frequency $\Omega_{i}$, referred to as natural frequency or inherent frequency. Therefore the distribution of natural frequencies and the distribution of net power $P_{i}$ are equivalent, up to a constant factor.

Finally, adiabatic simulations (upsweep of $k$) are performed to measure the level of synchronization in the network starting from the asynchronous state towards the partially synchronized state. In particular the rescaled coupling strength $k$ is increased from $k=0$ to $k=60$ in steps of $\Delta k=0.4$ (from $k=0$ to $k=60$ in steps of $\Delta k=0.2$) for the bimodal Gaussian distribution (real-world distribution, respectively). Specifically, for the bimodal Gaussian distribution with $\alpha=5/6$ s*-1*, $I_{i}=40\times 10^{3}\text{kg}\text{ m}^{2}$, $\omega_{\text{G}}=2\pi\cdot 50$ Hz and $\Delta k=0.4/\Delta t$ one obtains $\Delta K=\Delta kI_{i}\omega_{\text{G}}\alpha\approx 21$ MW, if a time unit $\Delta t=0.2$ s is considered.

Appendix B: Generality of the results

In order to show that the efficiency of our proposed control strategies is not restricted to the setups shown in the main text, we will present additional results: (a) keeping the setups shown in the main text, but analyzing different operating points and different configurations by considering different coupling strengths; (b) taking into consideration a different topology.

Different operating points

In this section we present the results for a different operating point, thus giving rise to a different configuration of solitary nodes. In particular, keeping the same setups presented in the main text, we show a comparison between the strategies (i) and (ii) obtained when the system is evaluated at different coupling strengths, thus investigating different working points with respect to the results shown in the main text. For the bimodal Gaussian distribution $P^{G}$, we investigate the state at $K\approx 565$ MW, which is a partially synchronized state found during the upsweep protocol, characterized by 19 solitary nodes. This configuration is unstable, with $\lambda_{1}=0.0144\pm 0.0005$. Regarding the real-world distribution $P^{R}$, the different working point that we have investigated is characterized by $K\approx 578$ MW, 19 solitary nodes and $\lambda_{1}=0.116\pm 0.005$.

The outcome of the control schemes is shown in Fig. 13. For the $P^{G}$ distribution strategy (i) requires the control of 2 solitary nodes to stabilize the system and 15 to synchronize, while strategy (ii) requires the control of 8 nodes to stabilize and 19 to synchronize the system. For the realworld distribution both strategies require one controlled node to stabilize. Synchronization is reached with 13 and 14 nodes using strategy (i) and (ii) respectively.

Italian grid

In this section we apply our control strategy to a different grid topology. The dynamics of the single node is still described by Eq. (1), but we now consider the Italian high-voltage (380 kV) power grid (Sardinia excluded), which is composed of N = 127 nodes, divided into 34 generators (hydroelectric and thermal power plants) and 93 consumers, connected by 171 transmission lines ita . This network is characterized by a quite low average connectivity $\left\langle N_{c}\right\rangle=2.865$, due to the geographical distributions of the nodes along Italy (see Fig. 14a). Since we have no access to a distribution of generator powers and nodal power consumption, we restrict the application of our method to the artificial distribution, using a bimodal Gaussian distribution (shown in Fig. 14b) with the same probability density function as the one used for the German grid (see Eq. 4 of the main text).

Like for the German grid, the synchronization transition is hysteretic (see Fig. 15a), but the formation of frequency clusters at different stages of the upsweep protocol is more pronounced since the local architecture favours a splitting based on the proximity of the oscillators. At $K\approx 461$ MW (middle black point of Fig. 15a) the system is partially synchronized and unstable ($\lambda_{1}>0$): it represents a big cluster of locked oscillators with zero average frequency and 20 unsynchronized whirling oscillators (see panel d). Besides the main frequency-synchronized cluster, two other clusters can be found: one with positive and one with negative average frequency, consisting of eight and five nodes, respectively. The remaining seven nodes are solitary. As before, we will take this state as an example to be controlled using our proposed strategies. For smaller coupling the system is unstable, but completely asynchronous (see panels b, c), while for larger coupling the system is (almost) completely synchronized (see panel e): one solitary state corresponding to the last node in Sicily hardly synchronizes due to the peripheric position in the network.

In Fig. 16 a comparison of strategies (i) and (ii) is presented. First of all, as for the German grid, the delayed feedback control is able to synchronize and stabilize the grid when enough nodes are controlled. Strategy (i), which controls preferably the most unstable nodes, sorted according to their Lyapunov vector component $\xi_{i}$, needs two nodes to stabilize and three controlled nodes to synchronize the system (see panel a). On the other hand, by employing strategy (ii), which orders the controlled nodes with respect to their frequency deviation $\left|\left<\omega_{i}\right>_{t}\right|$, the control of one node is required to stabilize, and two controlled nodes to synchronize the system (panel b). In both cases a remarkably small fraction of the 20 whirling nodes has to be controlled to gain the suitable conditions for operating power grids, thus highlighting the role played by solitary nodes in driving the network dynamics.