Deciphering the imprint of topology on nonlinear dynamical network stability

Many biological, chemical and technical systems of N coupled oscillators can effectively be described by the Kuramoto model [30, 34], which is given by the following temporal evolution law for an oscillator i's phase ${\phi }_{i}$:

$$\begin{eqnarray}&&\dot{{\phi }_{i}}=\displaystyle \frac{{P}_{i}}{{\alpha }_{i}}-\displaystyle \frac{1}{{\alpha }_{i}}\displaystyle \sum _{j=1}^{N}{K}_{{ij}}\sin ({\phi }_{i}-{\phi }_{j})\qquad \forall \ i=1...N,\end{eqnarray} \tag{ 2.1 }$$

where ${P}_{i}/{\alpha }_{i}$ denotes the oscillators natural frequency and K_ij reflects the coupling strength between nodes i and j, satisfying ${K}_{{ij}}={K}_{{ji}}\gt 0$ if nodes i and j are connected and ${K}_{{ij}}=0$ otherwise. Hence, the K_ij define the topology of a (weighted) network of oscillators. ${\alpha }_{i}$ denotes a damping coefficient of node i.

In many contexts the dynamics of coupled oscillators is described more accurately when inertia is accounted for. This is obtained by the second-order Kuramoto model which additionally describes the temporal evolution of the node's frequency ${\omega }_{i}$:

$$\begin{eqnarray}&&\dot{{\phi }_{i}}={\omega }_{i},\,\end{eqnarray} \tag{ 2.2 }$$

$$\begin{eqnarray}&&\dot{{\omega }_{i}}={{P}}_{i}-{\alpha }_{i}{\omega }_{i}-\displaystyle \sum _{j=1}^{N}{K}_{{ij}}\sin ({\phi }_{i}-{\phi }_{j}).\end{eqnarray} \tag{ 2.3 }$$

This model was shown to effectively describe the dynamics of synchronous machines within power grids [35–39] where the ${\phi }_{i}$ denote the phase angles and the ${\omega }_{i}$ the frequency deviations from the grid's rated frequency. In this context the _Pi correspond to the net power input at node i which is positive (negative) for generator (consumer) buses and the ${\alpha }_{i}$ describe the strength of the electro-mechanical damping and droop control. The coupling coefficients K_ij correspond to the capacity of the transmission lines of the power grid.

For typical parameter values, the dynamical system given by (2.2) and (2.3) features a stable steady state $({\phi }^{* },{\omega }^{* })=({\phi }_{1}^{* },\ldots \,{\phi }_{N}^{* },{\omega }_{1}^{* },\ldots ,{\omega }_{N}^{* })$ with constant phase angles ${\phi }_{i}^{* }-{\phi }_{j}^{* }$ and vanishing frequency deviations ${\omega }_{i}^{* }=0$ at every node i, given by a solution the following system of nonlinear equations:

$$\begin{eqnarray}&&{P}_{i}=\displaystyle \sum _{j=1}^{N}{K}_{{ij}}\sin ({\phi }_{i}^{* }-{\phi }_{j}^{* })\qquad \forall \ i=1...N.\end{eqnarray} \tag{ 2.4 }$$

This state corresponds to the desirable synchronous operating mode of the power grid. However, depending on the parameter choices, there might be also undesirable non-synchronous states which correspond to attracting limit cycles within the phase space of the dynamical system.

Menck et al showed that for a single-node, connected to an infinite power grid, the limit cycle solution can be approximated as

$$\begin{eqnarray}&&{\omega }_{\mathrm{LC}}(t)\approx \displaystyle \frac{P}{\alpha }+\displaystyle \frac{\alpha K}{P}\cos \left(\displaystyle \frac{P}{\alpha }\,t\right),\end{eqnarray} \tag{ 2.5 }$$

provided that $| P| /{\alpha }^{2}\gg 1$ and ${| P| }^{2}/{\alpha }^{2}\gg K$ [31].

In order to generate a representative ensemble of spatially embedded power grid topologies, we used a suitable random growth model [32]. It aims at generating synthetic network topologies that reproduce topological properties of real-world power grids and other spatially embedded infrastructure networks. The two-phase algorithm starts with a minimum spanning tree of size N₀ to which further nodes and (redundant) lines are added iteratively. The network growth is subject to a heuristic redundancy-versus-cost optimisation, which takes not only the line lengths but also additionally-created redundancy in the form of alternative routes into account. The growth model parameters have been set to (${N}_{0}=1$, $p=1/5$, $q=3/10$, $r=1/3$, $s=1/10$), where p, q and s are probabilities related to the creation and splitting of lines, and r specifies the redundancy-versus-cost trade-off. For a detailed explanation of all parameters, we refer to [32]. For this choice of parameter values the randomly generated networks match characteristics of real-world power grids, for instance the sparsity with a mean degree of about $\bar{d}\approx 2.7$. The ensemble consists of M = 50 random networks of size N = 100, an exemplary topology is shown in figure 1.

Zoom In Zoom Out Reset image size

We further make use of a topological classification scheme for nodes, which particularly distinguishes nodes located in or adjacent to tree-shaped parts of the networks. We start by giving concrete definitions for trees and tree-shaped parts and their roots:

Definition 2.1 A graph $G=(V,E)$ is called a tree if it is connected and has no cycles.

Definition 2.2 Let $G=(V,E)$ be an undirected graph that is connected but not a tree. A tree-shaped part is an induced subgraph ${T}^{\prime }=({V}^{\prime },{E}^{\prime })$ of $G$ which is a tree and is maximal with the property that there is exactly one node $r\in {V}^{\prime }$ that has at least one neighbour in $G-{T}^{\prime }$. $r$ is called the root of ${T}^{\prime }$ and has degree $d(r)\geqslant 3$.

The union of all nodes located in tree-shaped parts is subsequently denoted by $T={\bigcup }_{i}{T}_{i}$ and called the 'forest' part of G, where the T_i are the different tree-shaped parts within G. The remaining parts of G are referred to as the bulk, denoted by $B=G-T$. A finer partition of T is achieved by distinguishing the root nodes R from the non-root nodes N. The non-root nodes can be further subdivided into the leaves $L=\{l\in N\,| d(l)=1\}$ which have degree one, and the inner tree nodes $I=N-L$ which are located between the root and the leaves.

We will see that for stability assessment, an even finer partition is useful for whose definition we first need to introduce the following properties of nodes in tree-shaped parts:

Definition 2.3 Given a tree-shaped part ${T}^{\prime }=({V}^{\prime },{E}^{\prime })$ of a graph $G=(V,E)$, the depth $\delta (x)$ of a node $x\in {T}^{\prime }$ is the length of the shortest path from x to the root of ${T}^{\prime }$. The root $r\in {T}^{\prime }$ has depth $\delta (r)=0$.

The height $\eta (x)$ of a node $x\in {T}^{\prime }$ is the length of the longest outward path from x to a leaf of ${T}^{\prime }$. All leaves $l\in {T}^{\prime }$ have height $\eta (l)=0$.

Note that the presented definitions of height and depth of nodes might appear counter-intuitive when applied to the picture of trees growing upward from the root. However, this terminology originates from the data structure of a 'tree' in informatics, which is typically depicted as growing downwards, and became standard in graph theory.

The smallest possible type of tree-shaped part consists of a root and some adjacent leaves. Such leaves are subsequently termed sprouts and form the class $S=\{x\in N\,| \eta (x)=0\wedge \delta (x)=1\}=\{x\in L\,| \delta (x)=1\}$. The leaves of larger tree-shaped parts are called proper leaves and form the class $P=\{x\in N\,| \eta (x)=0\wedge \delta (x)\gt 1\}=\{x\in L\,| \delta (x)\gt 1\}$. Note that $G=B+R+N=B+R+I+L=B+R+I+P+S$.

Finally, the group of sprouts can be separated into those which are connected to high-degree roots, called dense sprouts ${S}_{d}=\{x\in S\,| {\bar{d}}_{{ \mathcal N }}\gt 5\}$, and those connected to rather low-degree roots, the sparse sprouts ${S}_{s}=\{x\in S\,| {\bar{d}}_{{ \mathcal N }}\lt 6\}$ where ${\bar{d}}_{{ \mathcal N }}(x)$ denotes the (average) degree of the neighbour (s) of x. For this last distinction, we chose the threshold of 5 so that the stability properties of the two groups are separated best.

For an efficient algorithm which provides both the partition of the nodes into the topological groups and their respective height and depth levels see appendix A.

The nodes in the exemplary network in figure 1 are coloured according to this classification and a representative node of each group is labelled accordingly. Definitions and total shares of the node categories in the ensemble of randomly generated network topologies are summarised in table 1.

Table 1.  Overview of names, symbols and definitions of the hierarchically ordered topological groups of nodes in tree-shaped parts of networks. The last column shows the shares of nodes of each category in the ensemble of the M = 50 randomly generated network topologies. More than half of the nodes belong to tree-shaped structures within the networks and about a quarter is given by leaf nodes. An exemplary network topology with the nodes coloured according to these groups is shown in figure 1. The simulation results shown in figure 5 are also coloured according to this scheme.

Group	Symbol	Definition	Share of all nodes
Bulk nodes	B	$\{x\in G\,| x\notin T\}$	48.0%
Roots	R	$\{x\in T\,| \exists b\in B\,:\{x,b\}\in E\}$	19.6%
Non-roots	N	$\{x\in T\,| x\notin R\}$	32.4%
Inner tree nodes	I	$\{x\in N\,| d(x)\gt 1\}$	7.2%
Leaves	L	$\{x\in N\,| d(x)=1\}$	25.2%
Proper leaves	P	$\{x\in L\,| \delta (x)\gt 1\}$	7.1%
Sprouts	S	$\{x\in L\,| \delta (x)=1\}$	18.1%
Sparse sprouts	Ss	$\{x\in S\,| {\bar{d}}_{{ \mathcal N }}(x)\lt 6\}$	12.9%
Dense sprouts	Sd	$\{x\in S\,| {\bar{d}}_{{ \mathcal N }}(x)\gt 5\}$	5.2%

The subsequently introduced measures assess the stability of a dynamical system with respect to large perturbations. They reflect global characteristics of the system's phase space which distinguishes them from the local perspective taken in conventional linear stability analysis in which only small (infinitesimal) perturbations are regarded [12, 40]. These measures are particularly suited for assessing the stability of power grids, since in this context large perturbations from the normal operating state are a common threat [41].

The first measure to be introduced, basin stability, is an indicator for the system's likelihood to asymptotically return to a desirable state following a large perturbation [17]. The second, survivability, in turn is sensitive to whether the transient behaviour after a larger perturbation remains within a desirable region of the system's phase space [18].

In complex networks, it is instructive to regard only localised perturbations originating at a single node which makes the stability measures node-wise quantities. However, it should be noted that for the stability assessment the response of the whole system with all nodes is relevant and hence the violation of the stability constraints might happen at other nodes than the perturbed one.

Both measures necessitate the specification of a probability distribution from which large (finite) perturbations are drawn. For the case of power grids modelled by (2.2) and (2.3) these are given by values chosen uniformly at random $(\delta \phi ,\delta \omega )\in [-\pi ,\pi ]\times [-{\rm{\Delta }}\omega ,{\rm{\Delta }}\omega ]$ which are added at ${t}_{0}=0$ to the state variables of a single node j while all unperturbed nodes are initialised to the desirable steady state:

$$\begin{eqnarray}&&{\phi }_{i}(0)={\phi }_{i}^{* }+{\delta }_{{ij}}\delta \phi ,\end{eqnarray} \tag{ 2.6 }$$

$$\begin{eqnarray}&&{\omega }_{i}(0)={\delta }_{{ij}}\delta \omega ,\,\end{eqnarray} \tag{ 2.7 }$$

where ${\delta }_{{ij}}$ denotes the Kronecker delta. Recall that ${\omega }_{i}^{* }=0$.

2.4.1. Basin stability

The concept of basin stability is applicable to multistable dynamical systems with state space X for which there exist attracting states distinct from the set of desirable attracting states, where the latter is denoted by ${X}^{\star }\subset X$ in the following. The basin of attraction ${{ \mathcal B }}^{\star }$ of ${X}^{\star }$ is given by all initial states from which the system asymptotically converges to the desirable attractor:

$$\begin{eqnarray}&&{{ \mathcal B }}^{\star }=\{x(0)\in X\,| x(t)\to {X}^{\star }\}.\end{eqnarray} \tag{ 2.8 }$$

When the perturbations are drawn from a certain subset ${X}^{0}\subseteq X$, it is instructive to define the basin of attraction restricted to this region: ${X}^{B}={{ \mathcal B }}^{\star }\cap {X}^{0}$ (see figure 2). This is especially the case if the phase space is not compact. In the case of perturbations drawn uniformly at random, the basin stability β is then simply the ratio of the volume of the (restricted) basin of attraction X^B to the overall region of perturbations X⁰ [17, 42]:

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\mathrm{Vol}({X}^{B})}{\mathrm{Vol}({X}^{0})}.\end{eqnarray} \tag{ 2.9 }$$

In other words, this quantity corresponds to the probability for the system to return to the desirable attractor after a perturbation from X⁰.

Zoom In Zoom Out Reset image size

For the synchronous machine power grid model described by (2.2) and (2.3) the desirable attractor ${X}^{\star }$ is identical to the set of synchronous states $({\phi }^{* },{\omega }^{* }=0)$, while there exist several non-synchronous attracting states which are undesirable. Since basin stability is determined for each individual node j, the region of perturbations is the subset

$$\begin{eqnarray}{X}_{j}^{0} & = & \{(\varphi ,\omega )\in X\phantom{\rule{thinmathspace}{0ex}}| -\pi \leqslant {\varphi }_{j}\leqslant \pi \wedge -{\rm{\Delta }}\omega \leqslant {\omega }_{j}\leqslant {\rm{\Delta }}\omega \\ & & \wedge \phantom{\rule{thinmathspace}{0ex}}{\rm{\forall }}\phantom{\rule{thinmathspace}{0ex}}i\ne j\phantom{\rule{thinmathspace}{0ex}}:({\varphi }_{i}={\varphi }_{i}^{\ast }\wedge {\omega }_{i}=0)\}.\end{eqnarray} \tag{ 2.10 }$$

Hence the single-node basin stability corresponds to the ratio of areas in the phase space cross-section spanned by the dimensions associated to the node j (see figure 2).

As the basin of attraction and its geometry are typically not known a priori and difficult to determine, especially in high-dimensional systems [43], β needs to be determined numerically via a Monte-Carlo method. For each node L = 200 independent perturbations have been chosen uniformly at random from X⁰_j and the corresponding trajectories simulated for t = 100 time units. The basin stability β of node j can be estimated from the number s of trajectories which asymptotically return to ${X}^{\star }$. More details are given in appendix B.

2.4.2. Survivability

In contrast to basin stability the survivability concept [18] presumes a desirable region ${X}^{+}\subseteq X$ which must not be left by a trajectory for a time t following a large perturbation in order to call the system 'survived'. The finite-time basin of survival X^S_t is given by the fraction of those initial states of the system which give rise to evolutions that stay within ${X}^{+}$ until time t:

$$\begin{eqnarray}&&{X}_{t}^{S}=\{x(0)\in {X}^{0}\,| x({t}^{\prime })\in {X}^{+}\ \forall \ 0\leqslant {t}^{\prime }\leqslant t\}.\end{eqnarray} \tag{ 2.11 }$$

The infinite-time basin of survival is obtained by taking the limit ${X}^{S}={\mathrm{lim}}_{t\to \infty }{X}_{t}^{S}$. The survivability σ of a system with respect to uniformly drawn random perturbations from the region X⁰ is then analogously given as the ratio (see figure 2)

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{\mathrm{Vol}({X}^{S})}{\mathrm{Vol}({X}^{0})}.\end{eqnarray} \tag{ 2.12 }$$

Hence it can be regarded as the probability for the system to remain within a desirable region ${X}^{+}$ after being hit by a perturbation from ${X}^{0}$.

For the operation of power grids it is necessary to keep the frequency deviations of all generators and consumers below a certain level [44]. Hence for the synchronous machine model the desirable region of the state space is given as

$$\begin{eqnarray}&&{X}^{+}=\{({\varphi }_{1},\ldots ,{\varphi }_{N},{\omega }_{1},\ldots ,{\omega }_{N})\in X\phantom{\rule{thinmathspace}{0ex}}| \,{\rm{\forall }}\phantom{\rule{thinmathspace}{0ex}}i\,:\phantom{\rule{thinmathspace}{0ex}}| {\omega }_{i}| \phantom{\rule{thinmathspace}{0ex}}\leqslant \phantom{\rule{thinmathspace}{0ex}}{\omega }^{+}\},\end{eqnarray} \tag{ 2.13 }$$

with ${\omega }^{+}\gt 0$ being the maximally tolerable frequency deviation. Again for each node j the perturbations are chosen at random from ${X}_{j}^{0}$ (see (2.10)) and the basin of survival is determined node-wise as ${X}_{j}^{S}\subseteq {X}_{j}^{0}$ (figure 2). Note that it is not of relevance at which node the frequency constraint $| \omega | \leqslant {\omega }^{+}$ is violated for the perturbation to count as not survived. Subsequently, the desirable region boundaries are chosen identical to the maximal perturbation level of the frequency deviations, ${\omega }^{+}\equiv {\rm{\Delta }}\omega$.

In order to estimate the single-node survivability ${\sigma }_{j}$ for each node, L = 200 trajectories with independent perturbations to the synchronous operating state, drawn uniformly from ${X}_{j}^{0}$ (see 2.10), have been simulated. The fraction of trajectories which did not leave the desirable region ${X}^{+}$ has been used as a statistical estimator for ${\sigma }_{j}$ (see appendix B for details).

Single-node basin stabilities and survivabilities have been estimated for all nodes in an ensemble of M = 50 randomly generated networks with N = 100 nodes each. Within each network half of the nodes act as net generators ($P=+1$), while the remaining nodes are net consumers ($P=-1$) such that there is an overall power balance, ${\sum }_{i=1}^{N}{P}_{i}=0$. Even though the synthetic power grid topologies generated by the random growth model are spatially embedded, the coupling strengths of the transmissions lines have been chosen uniformly to K = 6 for simplicity. The damping coefficients have been set uniformly to $\alpha =0.1$. The uniform choice of (i) the distribution of generators and consumers, (ii) the coupling strength of the transmission lines and (iii) the damping coefficients allows highlighting purely topological effects in the output data. Note, that for this choice of parameters the preconditions for the approximation of the limit cycle position, equation (2.5), are met.

Firstly, we investigated interrelations between the node-wise stability measures and topological properties of the perturbed node. We highlight topological effects in the example of a node's degree d and shortest-path betweenness b [3]. They turned out to reveal the most prominent insights for transient and asymptotic stability, respectively, and are basic established local/mesoscale network characteristics. Furthermore, this choice facilitates a comparison with previous findings [31, 32].

3.1.1. Basin stability

We regard different maximal perturbation levels ${\rm{\Delta }}\omega =2.5$ to 12.5, and observe that as expected the mean basin stability of all $N\cdot M=5000$ nodes decreases with growing perturbation levels ${\rm{\Delta }}\omega$.

In order to detect which topological node characteristics influence asymptotic stability, the basin stability scores are regarded in dependence of the degree d of the node which is defined as the number of neighbouring nodes (figure 3(a)). While there is no significant dependency for lower perturbations levels (${\rm{\Delta }}\omega \leqslant 5.0$) there is a slight increase of β with d for larger perturbation levels, for which there are, however, relatively large standard deviations in the stability estimates. Hence, degree alone is a weak predictor for basin stability which is in line with previous findings [31, 32].

Zoom In Zoom Out Reset image size

There is, however, a characteristic dependence of basin stability on the betweenness b of the node, which is defined as the number of shortest paths between all pairs of nodes within the network which are passing through the regarded node (figure 3(b)). While there is no general trend for any perturbation level, distinctive down-peaks of basin stability are observable at certain betweenness values (as illustrated by Menck et al [31]). These particular values of b correspond to nodes which lie within tree-shaped parts of the network (see the green-coloured nodes of the 'inner tree node' category shown in figure 1). A similar dependency was found in [31] for a maximal perturbation level of ${\rm{\Delta }}\omega =100$. Hence we were able to qualitatively reproduce these findings for considerably smaller perturbation levels which appear more realistic by comparison to real-world cases.

3.1.2. Survivability

Next, we studied the transient stability against large perturbations as measured by single-node survivability, again for different values of the maximal perturbation strength ${\rm{\Delta }}\omega$ (figure 4). The generally larger survivability scores achieved with increasing perturbations can be explained by the fact that the desirable region boundaries ${\omega }^{+}$ are increased simultaneously with ${\rm{\Delta }}\omega$ (see section 2). Another convention, which is not followed in this study, would be to hold the desirable region (${\omega }^{+}$) constant when increasing the perturbations (${\rm{\Delta }}\omega$). In this case the average survivability of all nodes would decrease since for ${\rm{\Delta }}\omega \gt {\omega }^{+}$ a certain fraction of trajectories would start outside ${X}^{+}$ and hence could not 'survive'.

Zoom In Zoom Out Reset image size

In contrast to basin stability, the single-node survivability features a significant negative correlation with the degree which is found for all levels of perturbations (figure 4(a)). This means that the frequency deviations within the network following a perturbation tend to be larger when nodes with a high degree are hit, thereby making the dynamics leave the desirable region. For example, for ${\rm{\Delta }}\omega =12.5$ the probability of the trajectory to survive a perturbation is close to 1 if a node with degree d = 1 is hit, while it lies below 0.5 for nodes with degree d = 13. For smaller perturbation levels the same trend is observable. Hence the degree may serve as a suitable predictor for the node's survivability.

For the betweenness measure the same but less striking trend is found (figure 4(b)). Survivability is generally decreasing with the node's betweenness. Inner tree nodes which are characterised by specific values of b do not stick out as dominantly as it is the case for basin stability. For lower perturbations (${\rm{\Delta }}\omega \leqslant 5.0$) we observe small up-peaks in σ at the specific b-values while for ${\rm{\Delta }}\omega =10.0$ the same down-peaks as with basin stability occur. As we show later, for ${\rm{\Delta }}\omega =10.0$ the survivability of inner tree nodes is strongly correlated with their basin stability which explains the occurring down-peaks. These observations reveal that nodes adjacent to tree-like structures are also crucial for predicting survivability, however, in a more subtle way compared to basin stability. Overall, betweenness alone is only a weak predictor for survivability, showing different features at different perturbation levels.

As for basin stability comparing survivability to other nodal network measures did not lead to more insights. Instead, a direct comparison of basin stability and survivability of the nodes was found to reveal non-trivial aspects of the system's dynamics and helped in identifying the node classification scheme introduced in 2.3.

While this section focuses on the general distribution patterns of basin stability and survivability, the stability characteristics of the different node classes are described in the next section (3.3). In order to reveal more insights about the interdependencies of asymptotic and transient dynamics as well as the relation to the network topology, we plotted the joint distribution of single-node survivability σ and single-node basin stability β for different perturbation levels ${\rm{\Delta }}\omega$ (figure 5). Each panel shows a total of $N\times M=5000$ individual estimates.

Zoom In Zoom Out Reset image size

Let us focus on the two marginal distributions first. For a rather low maximal perturbation level (${\rm{\Delta }}\omega =5.0$) the distribution of basin stability is extremely skewed with $73.0 \%$ of nodes featuring $\beta \geqslant 0.95$. The survivabilities in turn show a widely spread bimodal distribution. This fact shows that survivability is generally much more influenced by topology than basin stability at lower perturbation levels.

At larger perturbation levels (${\rm{\Delta }}\omega \geqslant 7.5$) the distributions of both basin stability and survivability are unimodal and skewed. For ${\rm{\Delta }}\omega =7.5$ still $62.4 \%$ of nodes feature a high asymptotic stability of $\beta \geqslant 0.9$. The distribution of transient stability is still wide but shifted towards larger values. While some groups of nodes show a strong correlation between β and σ, the overall Pearson correlation coefficient is close to zero ($\rho =0.07$).

This picture is reversed when looking at the highest perturbation level of ${\rm{\Delta }}\omega =12.5$. Here $62.9 \%$ of nodes have survivabilities of $\sigma \geqslant 0.9$ and the distribution of basin stability is rather widespread. It is known that at considerably larger perturbation levels (${\rm{\Delta }}\omega =100$) also basin stabilities show a widespread multimodal distribution [31, 45]. Hence basin stability might be a more useful measure for the system's stability if very large perturbations are expected.

The intermediate case with ${\rm{\Delta }}\omega =10.0$ is particularly interesting. Here the distributions of β and σ are very similarly shaped and numerous nodes show a strong correlation. However, there are also many nodes for which there is no such correlation, with the overall Pearson correlation coming out at $\rho =0.41$.

The findings presented so far are in line with the reasoning of Hellmann et al who state that basin stability and survivability are generally not correlated and hence the asymptotic behaviour does not allow conclusions about the transient one [18]. We now want to achieve a more precise statement by putting a stronger focus on the patterns observed in the joint distribution of the stability measures in figure 5.

In order to decipher how these patterns relate to topological characteristics, it turned out to be helpful to focus on nodes located inside or adjacent to tree-shaped parts of the network and to distinguish several types of nodes by how far inside the tree they are located (see section 2.3 and table 1). This is also suggested by the findings presented in section 3.1.

Each of the previously defined classes features typical characteristics regarding asymptotic and transient stability at different perturbation levels (figure 5). We begin with discussing ${\rm{\Delta }}\omega =5.0$. In this case, for all nodes the survivability is lower than the basin stability. This indicates that there are no undesirable attractors within the survival region ${X}^{+}$. Hence basin stability sets an upper limit to survivability. All trajectories which stay within the desirable region ${X}^{+}$ have to converge to the synchronous operating state of the grid, ${X}^{\star }$ [18]. This is expected as a desynchronization is expected to lead to oscillators rotating with their natural frequency, which here is ${\omega }_{\mathrm{LC}}=P/\alpha =10.0$ (see (2.5)).

At rather low perturbation levels (${\rm{\Delta }}\omega =5.0$) the nodes forming the upper mode of the survivability distribution are the leaves. Dense sprouts show lower survivabilities than sparse sprouts, indicating that also a low neighbour degree is beneficial for survivability. The lower mode of the distribution is formed by bulk, root and inner tree nodes. As expected from the degree dependence (see figure 4) inner nodes show higher survivabilities than root nodes which are most critical. Hence, the transient dynamics following a perturbation of a root node tend to exhibit large frequency deviations, leading to a transgression of the boundaries of the desirable region.

For ${\rm{\Delta }}\omega \geqslant 7.5$ we observe surprising new behaviour. While perturbations at most nodes still behave as if there was no undesirable attractor in the survival region, perturbations originating at a dense sprout almost all have $\sigma \gt \beta$. This means a novel asymptotic state or a very long transient, is reached, in which the system is not synchronised, but as the system does not leave the survival region ${X}^{+}$, no node is swinging at its natural frequency ${\omega }_{\mathrm{LC}}$ either. Exemplary trajectories for perturbations at a dense sprout are shown in figure C2 in appendix C and reveal that we indeed observe a novel asymptotic state, with a solitary desynchronized oscillator not swinging at its natural frequency. To our knowledge, this state has not previously been observed in these systems, for example in the bifurcation studies of [46, 47]. Indeed we expect that this state would be very hard to observe, if the initial conditions are drawn fully randomly, and not localised at individual nodes.

For sparse sprouts, proper leaves and inner tree nodes β and σ are strongly positively correlated. This means that for these nodes trajectories which leave the desired region tend to converge to non-synchronous states. In other words, the basin of attraction of the synchronised state is entirely contained in ${X}^{+}$. Root nodes show a pattern contrary to dense sprouts with a wide range of survivabilities and rather high basin stabilities.

The patterns yielded at ${\rm{\Delta }}\omega =10.0$ are very similar to those at ${\rm{\Delta }}\omega =7.5$. As the frequency of the limit cycle fluctuates around $\omega =10.0$, drawing the boundary exactly there means that still all nodes that fully desynchronise, and go to their natural frequency, will hit the survivability threshold. Now a few nodes besides the dense sprouts feature $\sigma \gt \beta$. A further anomaly that can be observed here is that the mean survivabilities of the leaves are smaller for ${\rm{\Delta }}\omega =10.0$ than for ${\rm{\Delta }}\omega =7.5$, opposed to the general observation of increasing survivabilities for larger values of ${\rm{\Delta }}\omega$.

Finally, at ${\rm{\Delta }}\omega =12.5$, the natural frequency of the oscillators is fully within the survival region, and a desynchronized node is characterised as having survived. Consequently, the majority of nodes has larger survivabilities than basin stabilities. There are numerous non-synchronous states whose trajectories lie completely within the desirable region ${X}^{+}$. The limit cycle trajectory is hence also prominent in the multi-node network model. While the leaves of the tree-like parts show a similar pattern in the σ–β-space, the inner tree nodes are clearly separated at slightly lower survivabilities and fairly lower basin stabilities. Nodes from the bulk feature the largest basin stabilities.

Independent of the perturbation level ${\rm{\Delta }}\omega$, dense sprouts tend to feature lower basin stabilities than sparse sprouts. Menck et al studied the dependence of basin stability on the (average) degree of the neighbour(s) for all leaves ('dead ends') but did not find a significant correlation [31]. It is our finer partition of leaves into proper leaves and sprouts, the neccesity of which was recognised by studying the joint plots, that allows more detailed insights.

Complementary to figure 5, the partition achieved using the defined topological classes of nodes is highlighted in three-dimensional representations of the joint distributions of basin stability and survivability at different perturbation levels (figure C1 in appendix C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The insights gained from this study are particularly relevant for both the stability assessment and the design of stable power grids. The employed nonlinear stability measures, basin stability and survivability, provide useful information on the asymptotic and transient stability of a grid against large perturbations, thereby surpassing the insights gained from local, linear stability assessments. The probabilistic measures have the benefits of being intuitive to understand and efficient to estimate for the overall system, irrespective of the complexity of the dynamics. They thus allow, for the first time, a systemic understanding of dynamical effects in the overall system.

Of particular relevance to the systemic stability of the power system are the novel asymptotic states we discovered, which primarily arise from perturbations at dense sprouts. The perturbed nodes are desynchronized and oscillate at a frequency much smaller than their natural. They are thus a pure network effect that can not be understood by studying individual machines. If such desynchronized states exist within sufficiently narrow frequency bounds, they would be a severe systemic risk to the power grid. They might be related to observed phenomena like Inter Area Oscillations, in which a deviation from perfect synchrony is observed. While an extensive analysis of these states is not part of this work, the topological characterisation of dynamic systems already suggests a mean to prevent them preemptively in the design of the power grid, by avoiding dense sprouts.

Beyond this novel dynamical phenomenon, our study revealed various interrelations between the pure topological properties of a node and its stability. Tree-like structures, which contain about half of the nodes in the network ensemble, were found to be characterised by stability properties which significantly distinguish them from the remaining bulk of nodes. We were able to identify a small set of topological groups whose nodes feature similar patterns of asymptotic and transient stability at various perturbation levels. The knowledge about these patterns can in turn be utilised to predict a node's relevance for the whole grid's stability knowing only its topological embedding within the network. Given a concrete real-world power grid topology, our classification scheme hence has the potential to facilitate the identification of particularly vulnerable nodes. While outside the scope of this work we expect that it will be possible to extend this identification to more nodes by means of statistic regressions taking network measures as input [45].

Menck et al found that the weak basin stability of inner nodes can be cured by reconnecting leaf nodes (there termed 'dead ends') to the grid [31]. This agrees well with our observation above to avoid sprouts and that 'detour' nodes are favoured by our stability measures [45]. Another strategy suggested by our findings would be to avoid high-degree nodes as these feature the worst survivabilities, independent of the perturbation level. By combining the above-mentioned rules, reconnection of leaves and avoidance of hubs with high centralities, a possible design principle for stable power grids would be to strive for rather homogeneous network topologies, characterised by narrow degree distributions. Using linear stability techniques such homogeneous topologies have also be found to be generally easier to synchronise [48], while tree-shaped structures apparently show rather bad synchronizability [49]. Interestingly, while tree-shaped structures feature both poor linear (synchronizability) and nonlinear (basin stability, survivability) stability properties, in small-world topologies synchronizability and basin stability were found to behave contrarily [17].

It should be pointed out that for our analysis we employed conceptual models of the transmission grid. A direct transfer of the findings to lower grid levels might not be valid. Another simplifying assumption which does not hold for real-world grids is the uniform distribution of both the power injections and loads. Whether rather heterogeneous power distributions would affect the findings has not been investigated, but might be particularly interesting for future research in the context of decentral renewable energy production. We are also aware that the maximal perturbation levels assumed for our simulations are rather high. Instead of assuming uniformly distributed perturbations it would be more realistic to assume a unimodal distribution in which small perturbations in both the phase and frequency deviations are more frequent and large perturbations are rare events. Realistic distributions of perturbations could, for instance, be derived from the research on intermittency in power fluctuations from wind or solar power systems [50, 51]. This work should thus be seen as an exploration of phase space properties like the structure of a state's basin of attraction, rather than a concrete study of realistic power grid models.

Beyond the direct implications for the resilience of power grids, our findings show the power of sampling-based dynamical methods, like survivability and basin stability. These are evaluated by sampling the system's dynamical reaction to non-small perturbations, here in phase space, but more generally also in parameter space. Combined with network measures and topological classifications, they provide general insights for complex dynamical networks. We expect that the crucial role of tree-shaped parts for the dynamics found here, is not specific to power grids but might rather be a general phenomenon in oscillator networks, infrastructure networks or even biological networks, e.g., in neural or genomic dynamics. The suggested topological classification scheme and the terminology for nodes in trees are independent of the nodal dynamics and are thus easily applicable to networks of different types.

As the employed synchronous machine model (or 'Swing equation') is equivalent to the Kuramoto model with inertia [30], the findings are particularly relevant for the general study of oscillator networks. The combined analysis of transient and asymptotic behaviour via basin stability and survivability allows indirect insights into the geometry of the system's phase space. The identification of nodes for which $\sigma \gt \beta$ allowed us to infer the existence and position of new types of limit cycles. Particularly, for the group of dense sprouts we found that non-synchronous states exist besides those given by the approximated solution to the single-node system (2.5).

In a first order approximation the coupling strength K in (2.5) is proportional to a node's degree d, while P and α are independent of topological characteristics. Hence the approximation (2.5) suggests that the amplitudes ${A}_{\mathrm{LC}}={PK}/\alpha$ of a limit cycle are proportional to the degree of the respective non-synchronously rotating nodes. This is a contributing mechanism to the low survivability of high degree nodes for very large perturbations. We suspect that it will be possible to gain further understanding of the relationship between topology and dynamics through more sophisticated approximations. Understanding the dependencies of the asymptotic spectrum of networked dynamical systems on the ambient topology in more detail will however require significant further work, both numerical and analytical.