Network-induced multistability through lossy coupling and exotic solitary states

The stability of synchronised networked systems is a multi-faceted challenge for many natural and technological fields, from cardiac and neuronal tissue pacemakers to power grids. For these, the ongoing transition to distributed renewable energy sources leads to a proliferation of dynamical actors. The desynchronisation of a few or even one of those would likely result in a substantial blackout. Thus the dynamical stability of the synchronous state has become a leading topic in power grid research. Here we uncover that, when taking into account physical losses in the network, the back-reaction of the network induces new exotic solitary states in the individual actors and the stability characteristics of the synchronous state are dramatically altered. These effects will have to be explicitly taken into account in the design of future power grids. We expect the results presented here to transfer to other systems of coupled heterogeneous Newtonian oscillators.

In the Methods section, we have discussed the numerical procedure of the average single-node basin stability estimation and motivated the selection of different parameter regimes for Eqn. 1. For the Scandinavian and the circle topology we have performed additional simulations for all regimes and report them here in an overview.
Supplementary Figure 3 accompanies the results of Fig. 4 in the main text with additional simulations of (1:3, K+)-parametrisations. Note that here we plot the ASBS results on a linear scale that is cut-off at 0.6 to obtain an alternative visualisation also of the regimes with low ASBS. Furthermore, Supplementary Figure 3j   The columns refer to a-j standard (D) b-k strongly (D+) and c-l very strongly damped (D++) regimes. The first two rows refer to the status-quo scenario 1:3 with a-c standard (K) and d-f strongly coupled regime. The last two rows refer to the prosumer scenario 1:1 with g-i standard (K) and j-l strongly coupled regime. The right axis in each subfigure shows the expected number of desynchronised nodes N d . For each value of α, the shares sum up to 1. Note the linear scale of the ordinate.    Phase at t = 0 Phase at t = 0 Phase at t = 1000 Phase at t = 1000 Phase at t = 1000

Homoclinic bifurcation of solitary oscillations and a continuation study
As it can be seen for instance in Fig. 3, there are regions in the α − K plane for which solitary oscillations do not exist. Our numerical studies give a strong indication that the underlying mechanism is a homoclinic bifurcation of the solitary limit cycle towards synchronisation. Furthermore, this is also suggested by our analogy with the single machine infinite bus model. The numerical procedure is as follows: the period of the solitary oscillation is estimated as the inverse time average of the phase velocity ω s at the solitary node with index s: for some large enough T such that we can be sure to measure several periods. Varying the phase lag α as the control parameter, we have used the period τ estimated at one parameter value to set T = 10τ for the next (updated) parameter value. The numerical integration is performed using SciPy's CVode solver with step size 10 −3 .
In the case of a homoclinic bifurcation, it is expected that period τ scales logarithmically with the parameter in proximity to the bifurcation point α c , i.e. τ − log |α c − α| (see for instance [1]).
Consider the phase portrait of the infinite bus model (P = 1, D = 0.1, K = 2, H = 1) close to the bifurcation point at α c ≈ 0.263 pictured in Supplementary Figure 6. The red circles correspond to 100 initial conditions equally spaced on the limit cycle. Firstly, the locus of the limit cycle is close to the stable manifold of the saddle point (determined by backward integration) before they merge into a homoclinic orbit at α c . Secondly, the distribution of points clusters near the saddle, indicating an increase of τ close to α c due to a slowing down of trajectories. The numerically determined values of τ are given in Supplementary Figure 5a. As a guidance to the eye, the solid yellow line depicts a logarithmic scaling as expected for the homoclinic bifurcation in the infinite bus model. For comparison, the dashed yellow line gives a square root scaling in the case of an infinite-period bifurcation. It appears that close to α c the logarithmic function is a better candidate for the scaling τ (α), the onset of the divergence is much steeper than expected for the square root case. Analogously, τ (α) is estimated for the circle model, here for the standard parametrisation and control. From Fig. 3 at K = 6, one expects a bifurcation at slightly negative α c ≈ −0.009 from normal solitary oscillations to synchronisation. Repeating our analysis from the infinite bus model for the circle indeed strongly suggests a homoclinic bifurcation (Supplementary Figure 5b). In summary, our numerical investigations and the analogy to the infinite bus approximation support the hypothesis that the global bifurcation curves in the parameter plane of the circle model are homoclinic.
Finally we give the results of a continuation study of the solitary states in Supplementary Figure 7. Note in particular that the region in which solitaries exist seems to continue to arbitrarily high coupling strength.

Global basin stability
Formally, the basin stability [2][3][4] and [5] S B of attractor A with basin of attraction B (A) is defined as where I B(A) (x) is the indicator function of the basin giving 1 if x ∈ B (A) and 0 otherwise. In this definition, S B depends on the choice of a compact subset R ⊆ X of the d-dimensional phase space X (except X is bounded, e.g. a Torus) and the selection of a probability measure µ on R, typically in terms of probability density ρ(x). ρ can be suitably chosen to incorporate the knowledge about the statistical nature of perturbations, for instance short circuit durations in power grids.
Set R can then be parametrised by a maximum frequency deviation ∆ω around the synchronous state φ * : We refer to S B as global basin stability because it can be interpreted as applying a perturbation to every node. If the perturbation at each node has an independent fixed probability to destabilize the system, this means that S B becomes exponentially small as one increases the number of nodes. Therefore, care needs to be taken when interpreting the results. In the Scandinavian topology it means that once a small fraction of nodes becomes vulnerable to perturbations, the global basin stability will be almost zero and lost in the sampling error. On the other hand the average single node basin stability that we consider in the paper uses perturbations that are local on the network. The likelihood to destabilize the system is thus proportional to the number of vulnerable nodes. Tracking how additional localized asymptotic states come in and cause more vulnerability in the network is thus more robustly done by using ASBS.
For the circle topology, where the system is homogeneous, all nodes become vulnerable at the same time. We show the dramatic reduction in the global basin stability (∆ω = 1) that results in Figure 8. Notably, even here solitaries play a role.
As mentioned above and discussed in the Methods section, the main tool of our investigation is based on networklocal perturbations and termed average single-node basin stability, ASBS. For completeness, we node that this version of basin stability can be written as above by constructing the spaces R k which consist of the states which satisfy: where φ * and ω global give the synchronous state. The perturbation space is then given by And the measure we chose is homogeneous within the region |ω| < ∆ω.

Analytic analysis of the circle model of the power grid
In this further discussion we consider an abstract power grid organisation, the so-called circle model, which demonstrates the main features we find responsible for prompt desynchronisation. In this model, the desired synchronised state exists and is stable. However, the model itself appears to be genuinely multistable. There arise indeed many other, the so-called solitary states in which one or a few power grid units start to behave differently with respect to all others perfectly synchronised such that not only the phase but also the average frequency deviates from the synchronous behaviour. Remarkably, these solitary states are Lyapunov stable too and hence, they cannot be excluded from the global system dynamics even by an essential increase of the inner network connectivity (as it would be naturally expected). Furthermore, it is found that this kind of network multistability arises under realistic operating parameters applied to power grids. If so, there is always a danger that sudden non-small disturbances associated e.g. with fluctuations in production or demand may trigger desynchronisation of one or several power generators/consumers even if the synchronous state is perfectly operating in a stable regime. While our model is topologically simplified and does not account for complication of the real power grid structure, we show through extensive numerical experiments that our predictions remain robust when including these effects for the Scandinavian power grid.
We define the circle power grid as follows (see also Supplementary Figure 1b). Let it contains N units and let a half of them, N/2, be identical generators (representing power plants) and the other half identical motors (representing consumers). Place all of them on a ring with the "alternative" disposition and assume that each unit is connected with R of its nearest neighbours to the left and to the right, and with an equal coupling strength K. Supplementary  Figure 1b shows the case of R = 2 which we have used for our numerical explorations. In the case that all units obey the same type of equation of motion, the Kuramoto model with inertia, of the form (i = 1, . . . , N ), where the term K sin α shifts the oscillator eigenfrequency Ω up or down depending on K and the parameter α.
Here, we use the notation In the lossless case α = 0 this term disappears and one comes to the classical Kuramoto model with bi-modal frequency distributionφ which dynamics is comparably simpler. Introducing non-zero α, however, makes the global network behaviour more sensitive to desynchronisation due to multiple appearance of solitary states, not avoidable even with large increase of the coupling strength K.
a. Synchronous state of the circle model To identify a synchronous state of the circle power grid model (Eqn. 5), re-write it in a two-group form denoting the generators by θ i and the consumers ψ i , i = 1, . . . , N/2: Assume θ 1 = · · · = θ N/2 ≡ θ and ψ 1 = · · · = ψ N/2 ≡ ψ, and subtract the second equation from the first one. The following equation is obtained for the phase difference η = θ − ψ: By simple trigonometry, it can be re-written as where a new coupling parameter µ = 2RK is introduced. This is an equation of motion in the two-dimensional invariant manifold indicated by the conditions θ 1 = · · · = θ N/2 ≡ θ, ψ 1 = · · · = ψ N/2 ≡ ψ as above. Synchronous states of the original N -dim Eqn. 5 are given by the roots of its right-hand side, which are corresponding to stable node O = (η O ; 0) and saddle S = (η S ; 0) of the in-manifold scalar Eqn. 9. The states are born in a saddle-node bifurcation at the bifurcation curve and they exist for all µ > µ sync , where node O transforms into a stable focus soon after the bifurcation at µ f oc = 2Ω The synchronous state we are looking for is given by the stable equilibrium O. The in-manifold stability, however, does not guarantee its transverse stability. Our numerics indicate, however, that the equilibrium O is also stable in the whole N -dimensional space of the original model Eqn. 5.
We conclude that the synchronous state of the circle model (Eqn. 5) is born in a saddle-node bifurcation as the coupling strength µ exceeds the bifurcation value µ sync and it preserves the stability with further increase of µ. In the state, generators and consumers create two phase-synchronised groups with the phase shift between them equal η O . At the saddle-node bifurcation curve µ = µ sync , the phase shift is η O = π/2 and it decreases to zero as µ → ∞.
Frequency Ω O of the synchronous state O is obtained from Eqn. 5 as: It follows therefore that synchronous frequency Ω O is zero in the lossless case α = 0, negative for α > 0 and positive for α < 0.
b. Solitary states in mean-field coupled circle model Solitary states can be derived analytically in the case of mean-field approximation of the circle model, i.e. when R = N/2 in Eqn. 5. For convenience, we write the mean-field model in the standard normalised form Here parameter α can be positive or negative. Due to the assigned symmetry, these two cases are reducing explicitly to each other by simple transformation α → −α, Ω → −Ω, φ i → −φ i , at which the generators and consumers exchange their network disposition. Therefore, the system dynamics appears to be equivalent for α < 0 and α > 0 and if so, only one of the cases should be derived. The symmetry property is also valid for the general, non mean-field form of the circle power grid model (Eqn. 5), which is clearly demonstrated by our simulations. However, it is not the case for non-symmetric generator/consumer disposition as it is normally the case in real power grid networks, e.g. in the Scandinavian power grid.
The synchronous state O of the mean field-model Eqn. 12 is obtained by the same procedure as in the previous section applied to the two-group representation Eqn. 7 including generators θ i and consumers ψ i : So, it has the form with a constant phase shift η O = θ − ψ between generators and consumers, where η O = arcsin 2Ω µ cos α . The synchronous state exists and is stable for all µ above the bifurcation curve as in 11.
Beyond the existence of synchronous state O which is phase and frequency synchronised, numerous solitary states of different configurations are typically arising in the mean-field model (Eqn. 12). They consist normally of several phase clusters, each rotating with its own frequency, manifesting in such a way a phenomenon of frequency clustering. This peculiar phenomenon differs from usual phase clustering well-known in (pure) phase oscillatory networks, i.e. the standard Kuramoto model. De facto, due to the presence of inertia, the system dynamics becomes much more involved. The main peculiarity is an enormous multistability, i.e., when the number of co-existing stable solitary and other states can grow exponentially with the system size N , moreover, in an essential domain of the system parameters. If so, spatial chaos [6][7][8] becomes an inherent characteristic of the model. We are coming apparently to a new reality, still not comprehended, in the field of coupled oscillators with inertia.
The simplest but not trivial solitary-type behaviour is caused by the appearance of a 1-solitary state. Here just one generator or consumer -let for definiteness it be a generator θ 1 .-splits off from the others, which are frequency synchronised consisting of two phase clusters: one of remaining generators and the other of all consumers, i.e., The examples of four different 1-solitary states for the mean-field model (Eqn. 12) are illustrated in Supplementary  Figure 9, for solitary generator (b,c) and consumer (a,d). Due to the model symmetry, the case with solitary consumers can be reduced to the generator one by simple parameter/variable change α → −α, Ω → −Ω, φ i → −φ i . We skip it in the section and concentrate further on the states with the solitary generator as in Supplementary Figure 9.
As it is illustrated in Supplementary Figure 9, there are two types of 1-solitary states in which a generator is detached. Actually, what happens in both cases is that the generator continues to rotate close to its eigenfrequency Ω = 10 in the situation when all others units synchronise with smaller (b) or larger (c) mean frequency. The first situation (b) may be considered as a "normal", as the desirable synchronous frequency should naturally be less then these of the generators and it equals zero if α = 0 and deviates to negative values as α > 0. Another 1-solitary state shown in (c) obeys a counter-intuitive property: the mean frequency of the synchronous cluster appears to be larger than 10. We call this state an "exotic" and it exists for negative α only.
The parameter region for both, "normal" and "exotic" c. Reduction to solitary manifold In the new variables (γ, θ, ψ) for 1-solitary state (Eqn. 14) the mean-field model (Eqn. 12) consists of only three equations: This is a reduced model of three coupled virtual "pendula" governing the dynamics in the respective 6-Dim solitary manifold. The pendula are linear but the coupling terms, however, are strongly non-linear (contain sinusoidal functions). This results in complex dynamical behaviour characterised by multiple solitary states and chaos. As coupling terms depend only on the phase differences, the system dimension can be further reduced by introducing phase difference variables ν = γ − ψ, η = θ − ψ. In the new variables, only two equations are left: This is a system of only two coupled "pendula", maintaining on a 4-Dim invariant manifold of the whole Ndimensional phase space. It is complicated for analytic study, as each pendulum as well as all coupling terms are nonlinear. However, the model is low-dimensional so it can be easily derived in a standard way with the use of numerical simulations to identify 1-solitary states for the N -dimensional mean-field model (Eqn. 12), up to the thermodynamic limit N → ∞. Note, that all in-manifold regimes which develop within the framework of the reduced system (Eqn. 15) should be tested on the transverse stability with respect to out-of-manifold perturbations, to guarantee the existence in the whole N -dimensional model.
Reduced low-dimensional systems for the k-solitary states can be obtained in a similar way as in the k = 1 case above. They are slightly different with the system parameters however, their dimensions are the same. Note finally, that the behaviour of solitary states becomes even more puzzled if the synchronised elements, variables θ i and ψ i in (Eqn. 16), split off by more than two clusters. The above reduced procedure still works, however, the dimension of the "solitary" manifold becomes larger then 4. Nevertheless, the study in this case is yet much simpler than the original N -dimensional model.