Introduction↩︎

Networks of coupled oscillators are indispensable for modeling natural and human-made systems. In fields as diverse as biology, neuroscience, ecology, physics, and engineering, important systems including the heart, the brain, food webs, coupled lasers, chemical reactions, power grids, and even firefly populations, are described as oscillators on complex networks [1]–[6]. The function of these systems is shaped in considerable parts by their connectivity, which is described by the network’s topology.

The paradigmatic model used to understand networks of oscillatory systems is the Kuramoto model [7], [8] and its variants. These consist of linear oscillators with a non-linear coupling, and feature extremely rich collective dynamical behavior, such as chimera states, frequency clusters, isolated desynchronization, and spatial chaos [1], [2], [4]–[6]. When including amplitude dynamics, one obtains Stuart-Landau oscillators. When combined with specific non-linear couplings, these have been shown to provide realistic models for the dynamics of renewable power grids [9], [10]. A complementary class of models are given by diffusively coupled non-linear oscillators, with van der Pol oscillators being a prominent example. Whereas the rich behavior and multistability in the Kuramoto setting are driven by the non-linear coupling [11] van der Pol oscillators also feature complex, multistable local dynamics.

The collective phenomena in all these systems are rooted in synchronization [3]. In some contexts, such as the brain, complete synchronization indicates dysfunction such as epilepsy. In others, such as power grids, it is a fundamental prerequisite for the system to work. Thus, a central question in the fields of networked complex systems is how robust the synchronous state is, either with the aim to design controls to disrupt it, e.g. deep brain stimulation [12], or to create a system that has favorable synchronization properties [13]–[16] (for power grids or neuronal networks).

In practice, the structure vs. function relationship is typically studied by considering correlations between network measures, which quantify network structure, and dynamical properties. Two prominent examples of this are i) explosive synchronization [17], in which the degree determines the synchronization properties, and ii) the probability that large localized perturbations desynchronize the network or lead to further failures [15], [18]. The latter is of particular importance in the context of power grids. Localized single component failures occur frequently. If these induce desynchronization of the whole network, large-scale blackouts are the consequence [19]. Understanding the probability that such a desynchronization event occurs turns out to be a particularly challenging task for network and complex systems’ science [19]. Despite a decade of effort since Menck et al. [15] found that dead ends have unfavorable properties, no clear picture has emerged [20]–[22].

Recently, an alternative approach to capture the relationship of structure and function has been developed: graph neural networks (GNNs), which are machine learning (ML) architectures specifically adapted to working with data on graphs. GNNs and their generalizations work well in different domains, such as epidemic spreading via social networks [23], and molecular properties [24]. Lately, they also succeeded in predicting the desynchronization probability in oscillator networks [25]–[27]. The performance reached by these GNNs approaches the level that would make them useful in real-world scenarios. However, it is challenging if not impossible to interpret GNNs and translate the results back to an understanding of the mechanisms in the system. The field of explainability methods for GNNs is still nascent [28].

Here, we curate a wide set of network measures, and systematically explore network science and ML-based approaches to the structure-function relationship. To study non-linear combination of network measures, we make use of node-wise machine learning, that is, we train conventional ML models to predict stability at a node solely based on the network measures associated to that node. This approach to featurizing the nodes using network science will be called NetSciML in the following. In contrast, GNNs consider the whole graph at once, see 1.

We study both transient and asymptotic stability of synchrony in the context of power grid topologies and Watts-Strogatz ensembles, featuring Kuramoto dynamics, Stuart-Landau like power grid models and van der Pol oscillators.

The central question we address is how robust the relationships identified are. Our central findings are:

1. In almost all cases we found that the relationship between network measures and stability properties can vary strongly as parameters of the network ensemble are varied.

2. Both GNNs and NetSciML can accurately predict a variety of stability metrics for a range of dynamical networks from the same ensemble. Individual network measures often cannot. Especially in situations featuring heterogeneous dynamical parameters and dynamics, hand-crafted features capturing the dynamical heterogeneity are required to achieve good predictions.

3. Even if GNNs and NetSciML achieve good predictive power for networks of the same size, drawn from the same ensemble, the predictions of GNNs and NetSciML can fail to generalize to other networks. We observe cases where NetSciML succeeds while GNNs fail and vice versa. We interpret the failure to generalize as a failure to accurately capture the underlying causes of instability.

Results↩︎

We now show how network measures, and non-linear combinations of network measures with high predictive power, relate to stability of synchrony for very different dynamical systems. We investigate three paradigmatic oscillator models, using appropriate notions of transient and asymptotic stability, on different classes of network structures using a wide range of curated network measures. Our main focus will be how the relationship between network measures and stability metrics changes as the network structure varies.

Oscillator models: First, we employ second-order Kuramoto-Sakaguchi oscillators [7], [8], [29]. These are linear oscillators with non-linear coupling. We use randomly sampled heterogeneous parameters and a more homogeneous power-grid-inspired parameter set. The latter is motivated by the fact that second-order Kuramoto oscillators are a paradigmatic model of the synchronization mechanism in power grids [30].

Second, we analyze van der Pol oscillators [31]. These are non-linear oscillators with linear, diffusive coupling. As opposed to the linear oscillators, here multistability and complex behavior is driven by the amplitude multistability of the individual oscillators.

Third, we examine realistic power grid dynamics, focusing on grid-forming inverters relevant for renewable energy integration. We make use of a Stuart-Landau like complex oscillator formulation of grid-forming inverters that has been experimentally validated [9], [10]. Dynamically they are higher-order linear amplitude-phase oscillators with non-linear coupling. (see ¿sec:sec:supplement95inverter95powergrids? for details).

Stability metrics: We assess both transient and asymptotic stability using probabilistic measures appropriate to the dynamical system [14], [18]. We investigate the stability of the synchronous state to random perturbations occurring at a single node. Thus we have two stability measures per node in each case.

Network structure: We examine two main classes of network topologies. First, we generate Watts-Strogatz networks [32] and vary average degree and rewiring probabilities. This way we can systematically explore both, the impact of sparsity, and the transition between regular, small-world, and random network regimes. Second, we analyze power grid topologies spanning a broad range of sizes, primarily generated using a random-growth algorithm that replicates key structural features of real electrical grids [33], supplemented by several authentic power grid networks.

Experimental Setups We study four setups:

• WS-Kura: Watts–Strogatz networks with inertial Kuramoto–Sakaguchi dynamics and heterogeneous node and edge parameters.

• WS-VDP: Watts–Strogatz networks with van der Pol oscillator dynamics and fully homogeneous parameters.

• PG-Kura: Two ensembles of synthetic power-grid like topologies (20 and 100 nodes), the power-grid topologies of four countries (DE, FR, GB, ES) and a large synthetic network of Texas generated using a different synthetic model.

• PG-Real: One ensemble of realistic synthetic power-grid topologies (between 70 and 80 nodes) and four larger girds generated the same way (128, 256, 512, 1024 nodes) with high-fidelity inverter-based dynamic models [10].

Network measures: We curated a set of established network measures, some of which have been previously used for studying synchronization properties. To account for the various dynamic heterogeneities of the dynamics, we augmented these by network measures derived from the systems parameters as well as properties of the synchronous state: In power grids and Kuramoto oscillators the synchronous state features non-zero phase differences on the edges \(\Delta \phi_{ij}\), giving rise to a non-zero "power flow" \(\text{Flow}_{ij}\). The system parameters used are the real valued coupling \(K_{ij}\), the adjacency matrix \(A_{ij}\), the power injection \(P^d_i\) and damping \(D_i\) for Kuramoto dynamics, and the complex valued \(Y_{ij}\) admittance, active and reactive power \(P^d\) and \(Q^d\) and an abstract "Dynamic type of node" variable that distinguishes different inverters for realistic power grids. The only setup with no extra dynamic network measures is WS-VDP, where all nodes are dynamically identical and the synchronous state has no spatial structure. The full set of network measures used is given in 2.

Machine Learning Predictors To investigate how informative combinations of network measures can be, we used the node-wise network measures as input for a gradient boosted tree machine learning model, trained to predict the stability metric under investigation (NetSciML). As a further point of comparison, we also trained GNNs to directly predict the stability metrics from dynamic parameters and graph structure (see 1). For WS we chose the average degree 6 and rewiring probability 0.5 as the training ensemble, for PG-Kura we chose the 100 node ensemble and for PG-Real we used the 70-80 node ensemble. We evaluated the predictors in-distribution, on networks drawn from the training ensemble but not used for training, as well as out of distribution on networks structurally different than the training ensemble.

The full set of results is presented in the Supplementary Material (see 1).

Pearson correlation of individual network measures and stability varies strongly↩︎

In almost all cases we found that the Pearson correlation between network measures and stability properties varies strongly, and can even flip sign, as parameters of the network ensemble are varied. In 2 we present the Pearson correlation between selected network measures and transient and asymptotic stability metrics in all four experiments. To confirm that the Pearson correlation does not miss non-linear relationships we also evaluated Theil’s U coefficient. Results for the full set of network measures with both Pearson and Theil’s U are provided in the supplemental material (1). Summary tables of the key results are provided in 3; these tables contain most of the detailed numerical values reported below.

In WS-VDP we observe a clear distinction between small and high degree with the boundary between degree 6 and 8, and low and high disorder with the boundary between 0.1 and 0.3 rewiring probability. The correlation of resistance distance centrality and transient stability is \(-80\%\) at rewiring probability \(0.5\) and mean degree \(8\), but \(+8\%\) at rewiring probability \(0.3\) and mean degree \(6\) (see 2). In the highly heterogeneous WS-Kura experiment we observe a smoother transition from low to high degrees. The universal Kuramoto order parameter depends on \(\Delta \phi_{ij}\), and is correlated to the transient and asymptotic stability at low degrees, but completely uncorrelated at higher degrees. Conversely the dynamical parameters \(P^d\) and \(D\) show the inverse behavior, they become more informative the higher the degree rises, with \(P^d\) showing no Pearson correlation but high Theil’s U at higher degree. In PG-Real various measures derived from the synchronous state, such as \(\max_j \cos (\Delta \phi_{ij})\) show weakly positive correlation for smaller grids, but become moderately negatively correlated for the larger grids, but there are no clear patterns to the observed changes. In PG-Kura and asymptotic stability we see average neighbor degree going from +0.45 for the Spanish grid to \(-0.52\) for the Texan grid (15). Resistance distance centrality and second order centrality show high correlation across most grids studied, but still vary substantially: from 0.76/-0.79 in the Spanish grid to 0.33/-0.32 for Great Britain (14).

Across all experiments, only very few network measures show a significantly non-zero correlation that is stable across topology variation. Those are the dynamic type of the node and the total inductance, which is a weighted degree measure, for PG-Real (6) and current-flow betweenness centrality for PG-Kura asymptotic (4).

Notably, among the strongly correlated network measures we find both, network measures specifically designed for synchronization dynamics, such as resistance distance centrality [34], as well as network measures that have not previously been considered in this context at all, most importantly second-order centrality[21], [35].

NetSciML and GNNs can predict stability well↩︎

3 shows the \(R^2\) score of linear regression using the best individual network measure, NetSciML and GNNs. The \(R^2\) score is evaluated on networks drawn from the same ensemble as the training networks. NetSciML achieves strong predictive performance for all experiments, demonstrating that a non-linear combination of multiple network measures can accurately predict diverse stability metrics across dynamical networks from the same ensemble, almost always significantly improving on the best individual network measure. We also see that GNNs typically outperform NetSciML with the sole exception of predicting asymptotic stability for WS-VDP. There NetSciML achieves 0.72, whereas GNNs only achieve 0.55 (3). Especially in WS-Kura, which features heterogeneous dynamical parameters and dynamics, hand-crafted features capturing the dynamical heterogeneity are required to achieve good predictions. Using only topological features, transient predictions reach an \(R^2\) of only 0.43 compared with 0.74 when only dynamical features are used. Combining all features yields only a marginal further gain, increasing \(R^2\) to 0.75 ([tab:WSKura32summary]).

As 4 shows for predicting asymptotic stability for PG-Kura, NetSciML can achieve its performance with considerably less training data than required by GNNs. When comparing the performances on both ensembles consisting of grids of size 20 or 100, the results suggest that the achievable performance depends more on the total number of nodes in the training set, than on the total number of grids or the grid size. This demonstrates clearly the utility of combining theory and data driven analysis in situations where data is comparatively sparse.

GNNs and NetSciML can fail to capture the root causes of instability↩︎

Even though NetSciML and GNNs achieve good performance when evaluated on networks drawn from the same ensemble, they can fail if evaluated on networks from other ensembles, indicating that they learn spurious correlations rather than root causes of instability. We evaluated the performance of the trained predictors for asymptotic stability for all experiments, and for transient stability for the PG experiments on all variations of the networks. Results are given in 5.

For power grid tasks (PG-Kura and PG-Real), GNNs achieve higher predictive performance than NetSciML; both tend to lose performance when predicting on other ensembles or real-world topologies, but GNNs appear considerably more robust. For PG-Kura asymptotic stability, within ensemble, GNNs and NetSciML achieve 0.88 and 0.82 \(R^2\). When evaluated on synthetic grids of size 20, these drop to 0.72 for and GNN but 0.38 for NetSciML. For PG-Kura transient stability, we see only a minor drop from 0.98 to 0.96 for GNNs when evaluating on the large Texas power grid, but NetSciML drops from 0.97 to 0.78.

For PG-Real, transient performance is flat and generalizes perfectly, indicating that in this case the predictors are actually capturing the underlying causes of instability. For asymptotic stability, GNNs achieve positive \(R^2\) across all grid sizes, while NetSciML does so only to size 128 and fails at larger scales.

For WS-VDP, neither NetSciML nor GNNs generalizes across ensembles with different average degrees, suggesting that neither method captures causal relations that transfer across regimes. Nonetheless, NetSciML retains its in-distribution advantage over GNNs for the same average degree. The shifts in network-measure correlations across ensembles already discussed, indicate substantial changes in dynamical behavior, which likely explain the difficulty both models face in learning transferable patterns.

For WS-Kura, NetSciML generalizes much better to high degree networks than GNNs, despite GNNs achieving much stronger in-distribution performance. At average degree 10, GNNs fail to predict dynamical stability, whereas NetSciML maintains strong predictive accuracy. Notably, NetSciML even performs better in some cross-ensemble settings than within ensemble. Across average degrees, several network measures exhibit clear structure, with strong correlations concentrated in either low- or high-degree regimes. The predictive results suggest that NetSciML correctly identifies dynamical patterns associated with larger average degrees; as these patterns become more dominant with increasing degree, predictive performance correspondingly improves.

It is also notable that NetSciML predictions of asymptotic stability using only topological features exhibit a similar pattern as GNNs here ([tab:WSKura32summary]), with a sharp reduction of predictive performance when going from degree 8 to 10. This indicates that GNN accurately identify topological causes for instability but miss the dynamic patterns that become dominant at larger degrees. This also aligns with the fact that GNNs generalize better to degree 4 networks.

As expected, the ability of NetSciML to generalize robustly corresponds to stable correlation patterns of the underlying network measures. For example for transient stability in power grids, node degree predicts stability across ensembles for PG-Kura ([tab:PGKura summary asymptotic,tab:PGKura summary transient]) and, and node type for PG-Real: (6). If the explanatory power of individual network measures varies dramatically, the network measure based non-linear predictors struggle to generalize.

Discussion↩︎

Much of the network theoretic analysis of asymptotic stability of power grid synchronization of the past decade implicitly assumed that network measures could capture the underlying causes for high or low basin stability, e.g., that dead ends undermine stability [15]. However, if network measures were causally determinative for stability, one would expect this relationship to be robust.

Instead we see that the relationship of stability metrics and network structure as captured by network measures is highly sensitive to the specific network ensemble the network comes from. Even relatively modest changes in the network ensemble, such as going from mean degree 6 to 8 to 10 can render observed correlations irrelevant or even completely misleading.

Similar observations hold for stability predictions using NetSciML or GNNs. Both can fail to predict stability in networks that are not dramatically different from the ones they are trained on. Correlations between network measures and stability, as well as GNN and NetSciML predictions should not be interpreted causally without careful validation. This poses a significant challenge for learning about dynamic phenomena from the analysis of such predictors using explainability methods.

For all these challenges, it is also notable that it always was possible to obtain good predictions of dynamic stability from structural information within the ensemble by drawing on a broad range of network measures. Using either method, we obtained \(R^2 > 0.5\) everywhere except for PG-Real Asymptotic. As expected, NetSciML was considerably more data efficient. We also observed that handcrafted features capturing heterogeneity of the synchronous state and the nodal dynamics are crucial to obtaining good performance in the NetSciML approach.

It remains unclear what determines the sensitivity of the network measure - stability relationship. WS-VDP is the most homogeneous experiment we ran, and showed the strongest sensitivity. It is plausible that dynamic heterogeneity attenuates the sensitivity. On the other hand we observed catastrophic failure to generalize in GNN predictions for WS-Kura, our most heterogeneous experiment.

An open question out of scope for this work is whether these sensitivities are typical more broadly for structure - function type relationships. It would be highly interesting to study dynamical properties other than stability, as well as dynamic networks that arise in other real-world contexts, especially neuronal dynamics and molecular dynamics. In the latter context, GNNs are already widely used. Hand-crafting features for nodes has also been explored, but to our knowledge not by using the systematic network-science based approach we introduce here. A central challenge remains generating appropriate datasets.

A further challenge lies in hand-crafting appropriate network measures to capture dynamic heterogeneity. No principled way to distill heterogeneous dynamic properties into network measures, or to select topological measures particularly pertinent to certain dynamic settings, exist. The fact that one of the most informative network measures for the much studied PG-Kura setting, second order centrality, was only serendipitously included here highlights how challenging this task is.

For the PG-Kura setting, the possibility to combine network measures and GNNs was explored in [36]. There it was shown that naively using the network measures as additional GNN input features is not a useful strategy unless it is known that there are no major distribution shifts between network measures. Instead a more robust approach was to use the GNN to predict both, stability and network measures, forcing the GNN to more robustly characterize the underlying graph.

From a practical perspective, the data efficiency of the NetSciML approach is promising for enabling predictions in data-sparse applications such as power grids. Central challenges remain the need for hand-crafted features and the lack of robust generalization. It remains an open question whether combinations of GNNs and NetSciML could also improve the data efficiency of GNNs.

Our results do not allow drawing strong conclusions on why NetSciML and GNNs sometimes fail to generalize for modest structural variations. Two essential possibilities are that they learn spurious correlations that change as the ensemble changes, or that the mechanism of instability itself actually changes. In the former case, no (direct) causal conclusions can be drawn from studying the predictors. In the latter, a careful analysis of the within ensemble predictions would still be causally meaningful. The results for WS-Kura indicate that in this case GNNs suffer from a failure of the first kind, but for the other experiments our results do not provide evidence either way. Thus, the central theoretical question raised by our findings is why we see such dramatic failures and sensitivity.

Methods↩︎

Networks and Network Ensembles↩︎

We consider two main classes of networks: (i) Watts-Strogatz networks and (ii) power grids with both synthetic and real-world topologies.

Watts-Strogatz networks↩︎

For the Watts-Strogatz networks [32], we systematically vary the mean degree and rewiring probability to generate ensembles with diverse structural properties. We consider four average degrees (4, 6, 8, and 10) and five rewiring probabilities (0.1, 0.3, 0.5, 0.7, 0.9).

For the Kuramoto oscillator experiments, we generate 2,048 networks per ensemble, except for rewiring probability of 0.5 and degree 6, where we generate 20,480 networks to provide larger training sets.

For the van-der-Pol oscillator experiments, 1,024 networks per ensemble are generated, except for rewiring probability of 0.5 and degree 6, where 21,504 networks are generated.

Ensembles of synthetic and real world power grid models↩︎

To systematically investigate the relationship between network measures and dynamic stability in power grid models, we utilize large datasets from [26], [27]. These datasets, specifically designed for machine learning applications, comprise 10,000 networks each for ensembles with 20 and 100 nodes, as well as a synthetic model of the Texan power grid with 1,910 nodes (TX).

To assess the generalizability of our ML models from synthetic to real-world topologies, we further analyze four high-voltage power grid models based on the actual network structures of France (FR), Great Britain (GB), Germany (DE), and Spain (ES). The French (146 nodes), British (120 nodes), and Spanish (98 nodes) topologies are sourced from [37], while the German grid (438 nodes) is provided by the German Institute for Economic Research (DIW) [38]. For all real-world grids, we parameterize nodes and edges following the procedure established for the Texan grid in [26], [27].

Realistic power grids↩︎

To further validate the robustness of our methods, we apply them to a more realistic power grid dataset introduced in [39]. This dataset features advanced node modeling, with loads represented by algebraic constraints and producers parameterized as inverters in three configurations. Line modeling incorporates losses, and admittances are based on the real-world German power grid. The synthetic grid generation framework is described in [40], which includes validation steps to ensure realistic topological and dynamical properties. The original dataset contains 1,000 grid samples for training, each with 70–80 nodes. To test generalization across grid sizes, we also conduct simulations on larger grids: G128 (128 nodes), G256 (256 nodes), G512 (512 nodes), and G1024 (1,024 nodes).

Dynamical models↩︎

We investigate three types of oscillator dynamics: (i) the second-order Kuramoto-Sakaguchi model, (ii) van der Pol oscillators, and (iii) inverter dynamics relevant for realistic power grids with high shares of renewable energy sources.

Second-order Kuramoto-Sakaguchi model: The paradigmatic second-order Kuramoto model [8] is given by \[\ddot{\phi}_i = P^d_i - D_i \dot{\phi}_i - \sum_{j=1}^n K_{ij} A_{ij} \sin(\phi_i - \phi_j - \alpha_{ij}), \label{eqKuramoto}\tag{1}\] where \(\phi_i\) is the phase angle at node \(i\), and \(\dot{\phi}_i\), \(\ddot{\phi}_i\) are its first and second time derivatives, respectively. The network topology is encoded in the adjacency matrix \(A_{ij}\).

In PG-Kura we use a homogeneous bimodal parametrization that is meant to capture conceptual properties of future transmission grids with distributed resources: \(P^d_i \in \{-1, 1\}\) indicates net-consumer or producer and is randomly drawn, damping \(D_i = 0.1\), and a relatively large, homogeneous overall coupling strength \(K_{ij} = 9\) and \(\alpha_{ij} = 0\). A synchronous state typically exists for this system; however, it is often not globally stable. For WS-Kura all parameters are drawn uniformly at random for every line/node, with \(\alpha_{ij} \in [0, 1]\), \(K_{ij} \in [0, 4]\), \(D_i \in [-0.15, -0.05]\) and \(P^d_i\) from a normal distribution of width \(1\) and centered at \(0\).

Van der Pol oscillator: The van der Pol oscillator is a classical nonlinear oscillator that exhibits both amplitude and phase dynamics. The equations for a network of coupled van der Pol oscillators are: \[\begin{align} \ddot{x}_i = - x_i + \mu (1 - x_i^2 + \alpha x_i^4 - \beta x_i^6) \dot{x}_i - \sum_{j=1}^n A_{ij} (x_i - x_j) \;, \label{eqVdP} \end{align}\tag{2}\] where \(x_i\) is the state at node \(i\). We considered a homogeneous population with \(\mu = 0.01\), \(\alpha = 0.1\) and \(\beta = 0.002\), a dynamically rich regime studied in [41].

Inverter dynamics: For realistic power grid modeling, we use inverter-based dynamics, which are crucial for grids with high shares of renewable energy. The modeling is based on the normal form, introduced by [9] and further validated by [10]. Importantly this model does not only consider frequency, but also voltage dynamics and are generalized Stuart-Landau equations. We have an oscillator \(v_i \in \mathbb{C}\) and scalar \(x_i\) at each node, the coupling is given in terms of a complex admittance matrix \(Y_{ij}\).

\[\begin{align} S_{i} &= v_i \sum_j \overline{Y}_{ij} \overline{v}_j - P^d_i - j Q^d_i\\ \frac{\dot{v}}{v} &= \mu_1 + \mu_2 |v|^2 v + \mu_3 S_i + \mu_4 \overline{S}_i + \mu_5 x\\ \dot{x} &= \mu'_1 + \mu'_2 |v|^2 v + \mu'_3 S_i + \mu'_4 \overline{S}_i + \mu'_5 x\; . \end{align}\]

\(P^d\) and \(Q^d\) are the real and reactive power injected at a node, \(S_i\) is the complex power imbalance. Passive load nodes are modeled by a constraint \(S_i = 0\) enforcing a constant active and reactive power draw. Following [39] we have three types of inverters corresponding to three different sets of coefficients \(\mu\).

Stability metrics↩︎

To ensure robust and meaningful results, we employ probabilistic stability metrics, considering both transient and asymptotic stability. The specific metric used depends on the task and system under investigation.

Asymptotic stability↩︎

The asymptotic stability measure we use is the single-node basin stability (SNBS) metric [15]. SNBS quantifies the probability that, after a large random perturbation applied to a single node, the system returns to synchrony. Mathematically, SNBS corresponds to the volume (with respect to the perturbation probability measure) of the basin of attraction of the synchronous state, restricted to perturbations at a single node. An example of such a phase space slice for PG-Kura is shown in 6 for a node with relatively low stability.

For inverter-based power grids, we extend the SNBS concept to account for both frequency and voltage stability. A configuration is only considered stable if it satisfies both criteria following a perturbation.

Transient stability↩︎

For Kuramoto power grids, we use the maximum frequency deviation metric from [27], defined as the largest frequency deviation during the transient period.

In the context of realistic power grids, we use the probabilistic ride-through probability (\(p_\mathrm{frt}\)) as introduced in [39], a variation of the survivability [18]. This metric measures the likelihood that a node or system does nto violate operational boundaries (i.e., “rides through”) after a disturbance such as a fault or sudden operational change. The stability limits are time-dependent, becoming stricter as time progresses to ensure the system returns to a stable state within a predefined window. For these models, both frequency and voltage stability are evaluated simultaneously.

For van der Pol oscillator networks, we use the Synchronization Norm considered in [34], [42], the time averaged variance of the state \(x_i(t)\).

Informativity measures↩︎

To quantify the relationship between network measures and dynamic stability, and the performance of predictors, we employ several informativity metrics: Pearson correlation, Theil’s U, and the coefficient of determination (\(R^2\)).

Pearson correlation (\(r\)) measures the linear correlation between two variables \(X\) and \(Y\). It is defined as \[r = \frac{\mathrm{cov}(X, Y)}{\sigma_X \sigma_Y},\] where \(\mathrm{cov}(X, Y)\) is the covariance and \(\sigma_X\), \(\sigma_Y\) are the standard deviations of \(X\) and \(Y\), respectively. \(r\) ranges from \(-1\) (perfect negative linear correlation) to \(1\) (perfect positive linear correlation), with \(0\) indicating no linear correlation.

Theil’s U is a normalized measure of mutual information between random variables. It is defined as the fraction of the total information describing of \(X\) we can predict using \(Y\): \[U(Y|X) = \frac{H(Y) - H(Y|X)}{H(Y)},\] where \(H(Y)\) is the entropy of \(Y\) and \(H(Y|X)\) is the conditional entropy of \(Y\) given \(X\). Theil’s U ranges from \(0\) (no information) to \(1\) (perfect information). We estimated U by binning stability and network measures into 25 bins.

Coefficient of determination (\(R^2\)) is used to assess the performance of NetSciML, GNNs and to compare the predictive power of network measures. It is given in terms of the mean square error \(\mathrm{mse}\): \[R^2 := 1 - \frac{\mathrm{mse}(\widehat y,y)}{\mathrm{mse}(y_\mathrm{mean},y)}, \label{eq:R2}\tag{3}\] where \(\hat{y}\) are the predictions, \(y\) are the target values, and \(y_{\mathrm{mean}}\) is the mean of the target values in the test set. \(R^2\) measures the proportion of variance in the data explained by the model, with \(R^2 = 0\) corresponding to a model that predicts the mean. Negative \(R^2\) indicates a performance worse than predicting the mean, and thus a complete misprediction.

Together, these informativity measures provide a comprehensive assessment of the relationships between network measures, predictors and dynamic stability, capturing both linear and nonlinear dependencies as well as model performance.

Network measures↩︎

A table listing the complete set of network measures and the models to which they apply is provided in 2. The table also reports each measure’s origin, and which experimental setting they apply to. Because many measures appear in multiple sources, we list at most two references per measure.

To capture the structural properties of the networks in this study, we assembled a comprehensive set of network measures informed by an extensive literature review. This collection includes measures previously proposed as relevant for oscillator stability, as well as additional features developed in the broader context of network science. We further expanded this set with custom measures and additional metrics available in the NetworkX software package [43].

One common application of oscillator stability analysis is power grids, so we place strong emphasis on identifying all relevant network measures from this domain. In one of the earliest applications of probabilistic dynamic stability assessment of power grids [15] found that leaf nodes correlate with low stability. Subsequently, [35] predicted single-node basin stability (SNBS) of power grids with logistic regression. They classified nodes as belonging to certain network motifs, such as tree-like subgraphs, and used these as inputs for the regression model, alongside the injected power \(P^d\) and the common network measures, degree, average neighbor degree, clustering coefficient, current-flow betweenness and closeness centrality.

[21] expand upon this work by studying survivability as well as SNBS, and by refining the topological classification scheme. To compute the proposed node categories [21], the open source software from GitHub [44] was used. Furthermore, we consider connections between nodes of walks of lengths higher than 1, by computing the row-sum of \(A^k\), where \(k=3\) would represent walks of length 3.

To identify local vulnerabilities in power grids with intermittent fluctuations, [34], [42] use centrality measures based on resistance distance and Kirchhoff indices. Feld et al. [45] construct very stable and very unstable power grids and study the correlation of basin stability with various power grid specific network measures, such as flow backup capacity, power sign ratio [46] and universal Kuramoto order parameter [47]. [48] predict re-synchronization after individual perturbations in synthetic grid topologies generated with the same model that we use [33]. As input features for their ML models, they use a few common network measures, such as eigenvector centrality and \(k\)-core index. Furthermore, they use the specific perturbation vector, and the so-called community inconsistency. The community inconsistency depends on a free parameter for which the authors do not provide a heuristic choice [49]. Since, furthermore, [48] find no large correlation of this measure with dynamic stability, we do not include it in our study.

The related problem of predicting critical links in power grids has been studied by Witthaut et al. [50]. For this problem, Titz et al. [51] used Gradient Boosted Trees to predict de-synchronization event after line failures at high accuracy. They include many network measures as features, however most of them are only applicable to lines.

To capture the heterogeneity of the dynamics and the synchronous state we further added hand-crafted measures where appropriate. We considered the dynamic parameters \(P^d\), \(Q^d\), \(D\), where they apply, and a categorical node type feature for PG-Real. For PG-Real, PG-Kura and WS-Kura, the synchronous state is determined by phase angles \(\phi\). To capture the local properties of the synchronous state we calculated edge wise \(\Delta \phi_{ij}\), \(\cos(\Delta \phi_{ij})\) and the total flow on the edge \(\text{Flow}_{ij} = K_{ij} \sin(\Delta \phi_{ij} - \alpha_{ij})\). We then constructed nodal measures by taking the sum, minimum, maximum, and median of these edge wise values for each edge incident on a node. For PG-Real, the flow is calculated as the magnitude of the average apparent power on the line, and we further include various measures that characterize \(Y_{ij}\), including the total and shunt inducance.

Dynamic stability from structure - NetSciML↩︎

We employ two supervised machine learning approaches to predict dynamic stability labels: graph neural networks (GNNs) and NetSciML. The GNN-based setup is described in the following section; here, we focus on the NetSciML approach.

To systematically assess the predictive power of network measures for capturing the structure–function relationship, we use gradient boosted trees to predict stability from the full set of network measures. Histogram Gradient Boosted Regressors were used for all predictions using the library scikit-learn [52].

Dynamic stability from structure - GNNs↩︎

Graph neural networks (GNNs) are artificial neural networks designed for learning relationships within data structured as graphs. Their adjustable internal weights can be optimized to suit the specific task at hand. They take the graph structure itself as input, along with optional features associated to the nodes and edges. GNNs can generate outputs in the form of global attributes that describe the entire graph, attributes specific to sub-graphs, or properties associated with individual nodes or edges.

To capture long-range dependencies—an important property for analyzing our problem—we focus on two GNN architectures: (a) Topology Adaptive Graph Convolution (TAGNet) [[53]], which incorporates information from small subgraphs and indirect neighbors within a single layer, and (b) Dirac–Bianconi graph neural networks (DBGNNs) [[54]]. In contrast to conventional GNNs, which propagate information diffusively (analogous to the heat equation), DBGNNs enable coherent long-range feature propagation. This mechanism has shown clear benefits for prediction tasks on power-grid networks.

The GNN layers used in this study are described in detail in ¿sec:appendix95GNN? and we rely on PyTorch Geometric [55].

Acknowledgments↩︎

All authors gratefully acknowledge Land Brandenburg for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research. Michael Lindner greatly acknowledges support by the Berlin International Graduate School in Model and Simulation based Resarch (BIMoS) and Christian Nauck would like to thank the German Federal Environmental Foundation (DBU) for funding his PhD scholarship. Further, the work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, HE 6698/4-1 and KU 837/39-2) and the German Federal Ministry for Economy and Energy (03EI1092A). Special thanks to Anna Reckwitz for sharing her code for computing the semi-analytic bound for survivability. AI Tools such as LanguageTool and ChatGPT are used to improve the grammar and wording on a sub-sentence level. We also want to thank Konstantin Schürholt for carefully reading the manuscript and his valuable feedback.

Data Availability Statement↩︎

The used datasets of the ensembles and the synthetic Texan power grid are available at: https://zenodo.org/record/6572973 and we will also publish all code for the training of the ML models (NetSciML and GNNs), the semi-analytic boundary, the results of the real power grid topologies, as well as code to generate the figures upon publication on Zenodo: https://doi.org/10.5281/zenodo.10686691 and Github: https://github.com/PIK-ICoNe/NetworkScienceML_paper-companion. It will be licensed under the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

1 Supplementary material↩︎

The supplementary material is structured as follows. 2 shows the full set of network measures. 3 contains a summary of the experimental results, followed by detailed results for each experiment ([sec:appendix_experiments_WSVDP,sec:appendix_experiments_WSKura,sec:appendix_experiments_PGKura,sec:appendix_experiments_PGReal]). 8 contains supplementary material on the ML methods.

2 The full set of network measures↩︎

1 lists all network measures used, along with details on their origin and application in this study.

3 Supplementary material: Summary of experimental results↩︎

This section provides a concise overview of the experimental outcomes, highlighting selected results. Additional results are presented in the dedicated sections for each experiment. Specifically, [tab:WSKura32summary] summarizes the WS-Kura results; [tab:VDP pearson selected nm,tab:VDP NetSciML] summarize the WS-VDP results; [tab:PGKura summary asymptotic,tab:PGKura summary transient] summarize the PG-Kura results; and 6 summarizes the PG-Real.

4 Supplementary material on WS-Kura experiments↩︎

We analyze Watts–Strogatz ensembles with mean degree \(d \in \{4,6,8,10\}\) and rewiring probability \(p_r \in \{0.1,0.3,0.5,0.7,0.9\}\), governed by the inertial Kuramoto–Sakaguchi model with strongly heterogeneous parameters. The varying parameters are the coupling strengths \(K_{ij}\), phase shift \(\alpha\), damping \(D\), and injected power \(P^d\). Asymptotic stability is quantified via single-node basin stability, defined as the fraction of perturbations for which the system returns to synchrony with frequencies within narrow bounds, while transient stability is assessed using a synchronization norm that measures the magnitude of desynchronization during recovery.

7 reports Pearson correlations for the network measures in WS-Kura, while 8 reports Theil’s \(U\). For a selected subset of network measures, 9 provides the corresponding Pearson correlations.

In [tab:WSKura32summary], correlations vary systematically across two regimes: a sparse/regular regime and a denser more random regime. A nonlinear combination of network measures (NetSciML), trained at \(p_r=0.5\) and \(d=6\), generalizes well within the dense regime but degrades markedly in the sparse regime. NetSciML attains strong performance using only dynamical features: Theil’s \(U\) for the injected power \(P^d\) indicates high informativeness in the dense regime but substantially reduced informativeness in the sparse regime. Using only three dynamical inputs, \(P^d\), \(D\), and \(\sum_j K_{ij}\), NetSciML already matches—or sometimes exceeds—the performance obtained with the full set of network descriptors. By contrast, in the sparse regime the universal Kuramoto order parameter correlates strongly with stability, indicating that properties of the synchronous state become increasingly predictive there.

Graph neural networks (GNNs) trained at \((d=6,\, p_r=0.5)\) achieve moderate extrapolation to the sparse/regular regime, but their loss in accuracy mirrors that of NetSciML ([tab:WSKura32NetSciML32GNN]). GNNs outperform NetSciML in the dense regime, and part of this advantage carries over when extrapolating; however, models trained on dense data still fail to capture additional mechanisms driving (in)stability that emerge in sparse networks.

Training and model details↩︎

The GNN used is a DBGNN with 2 layers, each comprimising 20 internal steps. The hidden channel dimension for both nodes and edges is set to 500. No regression head is applied. The model has The model is trained for 2000 epochs with a batch size of 50. OneCycle learning rate scheduling is used with cosine anneal strategy, the final div factor is set to \(9694991443.141111\), the initial div factor to 9363750, and a maximum learning rate of 180.978942.

11 reports the mean and standard errors over multiple GNN initializations.

5 Supplementary material on WS-VDP experiments↩︎

We consider Watts–Strogatz ensembles with mean degree \(d \in \{4,6,8,10\}\) and rewiring probability \(p_r \in \{0.1,0.3,0.5,0.7,0.9\}\). Node dynamics are van der Pol oscillators with homogeneous parameters chosen to yield a bistable regime. Asymptotic stability is assessed via single-node basin stability (return of frequency to narrow bounds), and transient stability via a synchronization norm.

11 reports Pearson correlations for the network measures in WS-Kura, while 12 reports Theil’s \(U\). For a selected subset of network measures, 2 provides the corresponding Pearson correlations.

For the VDP dynamics, correlations between network structure and stability vary markedly across ensembles 2. For instance, the correlation of resistance distance centrality with transient stability is \(+0.08\) for \(d=6\), \(p_r=0.3\) but \(-0.8\) for \(d=8\), \(p_r=0.5\), indicating pronounced sign changes and limited consistency across regimes.

Models trained at \((d=6,\, p_r=0.5)\) generalize unreliably: neither GNNs (3) nor NetSciML transfer well beyond increases in rewiring probability. Overall, the GNN underperforms NetSciML even in-distribution, and both approaches fail to capture stability mechanisms that shift with topology in sparser or more regular networks.

Training and model details↩︎

The GNN used is a DBGNN with 2 layers, each comprising 10 internal steps. The hidden channel dimension for both nodes and edges is set to 100. A linear layer of dimension 100 follows the GNN convolutions. The model consists of 111,001 parameters. Dropout is set to 0. The batch size is 32, the learning rate is \(1 \times 10^{-5}\), and the model is trained for 2000 epochs. Training on 5 seeds consecutively with an H100 GPU takes approximately 15 hours.

15 reports the mean and standard errors over different initializations of the ML model.

6 Supplementary material on PG-Kura experiments↩︎

We analyze synthetic grids with 20 and 100 nodes, a large synthetic Texas system (approximately 2000 nodes) generated via a distinct procedure, and several real-world topologies. Node dynamics follow an inertial Kuramoto model with nearly homogeneous parameters; the only node-specific variation is the injected power \(P^d\). Asymptotic stability is assessed via single-node basin stability (return of frequency to narrow bounds), and transient stability via the maximum frequency deviation.

Correlations between structural descriptors and stability vary substantially across systems. For example, the correlation between average neighbor degree and asymptotic stability is \(0.45\) in Spain but \(-0.52\) in Texas, whereas second-order centrality and resistance-distance centrality remain comparatively stable (15). 14 reports the most strongly correlated network measures, and 15 the network measures with the highest variance.

Transient stability is predicted well by NetSciML 16. For asymptotic stability, GNNs trained on the 100-node ensemble generalize robustly to real-world grids (16). A NetSciML model without synchronous-state (here also the power-flow) features attains comparable accuracy within the 100-node ensemble but fails to generalize to real-world topologies; NetSciML augmented with the full feature set improves generalization yet still lags behind GNNs.

These results indicate that the power-flow structure is a key determinant of stability. Network measures combined with the single dynamical feature \(P^d\) do not recover this information. Even with power-flow features, NetSciML generalizes less robustly than GNNs: performance degrades on both smaller (20-node) and much larger (Texas) systems, whereas GNN accuracy remains largely stable across sizes, with the notable exception of GB.

17 reports the performance with fewer training samples. This table includes OMP3, which selects the subset of three network measures that collectively delivers the best predictive performance (see ¿sec:sec:app95best95subset95selection? for details).

Details on the study of power grids with Kuramoto oscillators↩︎

As a key challenge, we analyze the capability of predicting the dynamic stability of larger grids generated by the same random process, an even larger grid generated by a different random process and real-world topologies.

This can be interpreted as measuring how well the prediction captures underlying causal relationships, rather than simply fitting statistical relationships specific to the ensemble. In power grid applications, this ability is of crucial practical importance. The computational cost of the dynamic simulations grows at least quadratically with the size of the grid. Using ML methods for predicting dynamic stability becomes useful for real-world application, if they can be trained on datasets of small synthetic networks, which are easy to simulate, while still performing well on large, complex real-world grids.

Ensembles of synthetic and real world power grid topologies modeled with Kuramoto oscillators↩︎

To compare statistical differences regarding the dynamic stability of the real topologies in comparison to the synthetic datasets, we visualize the histograms of SNBS in 7. The differences in distribution are caused by the differences regarding their generation. The grids of the 20 and 100 node ensembles are generated using a spatial random growth model [33] tuned to infrastructure and power grid properties [33]. The Texan grid is generated by a process that aims to ensure good electrical properties of the generated grid [59], [60]. Lastly, the real topologies have their own unique characteristics. The German grid, which is the largest among these four, has a small third peak for low SNBS values which is similar to the synthetic Texan grid. This might indicate that the three modes are not a unique property of the Texan grid, but a property that only emerges for grids of certain sizes. The GB grid is less stable than the other topologies. Compared to the other real and synthetic grids, there are considerably fewer nodes that are stable against all perturbations. This suggests systematic differences in the structure of the GB grid and the other grids. Therefore, we expect the GB grid to be particularly challenging for extrapolation tasks.

The datasets also contain alternative prediction targets related to dynamic stability, notably the maximum frequency deviation throughout the transient. In contrast to SNBS, NetSciML predicts this measure at high accuracy and generalizes well to other ensembles.

Complexity and training details of the used ML models↩︎

GNN models require considerable training time, especially when conducting hyperparameter studies. Training of NetSciML is significantly faster. The NetSciML approaches based on gradient-boosted trees and linear regression have negligible training times of less than 2 CPU minutes. By far the most significant computational effort is required for the hyperparameter studies of the GNNs. For this analysis, it is common to analyze hundreds or thousands of possible configurations. In contrast, hyperparameter studies considerably less resource-intensive for gradient boosting trees (GBT). This underscores the significant advantage of employing the NetSciML approach, which incurs substantially lower computational costs. Once the models are trained, ML approaches have the huge advantage of negligible evaluation times. This is especially helpful when evaluating unknown grids of the size of the Texan grid. For such a large grid, the required Monte-Carlo simulations take roughly 15,000 CPU hours [27] when aiming for an accuracy comparable to current GNNs, whereas the evaluation takes only a few seconds for ML methods.

6.0.0.1 Training and model details

We use the trained The TAGNet model from [27]. The TAG model has has 13 layers and consists of 415,320 parameters. The number of channels is optimized per layer and varies between 30 and 141. No dropout is used.

[tab:PG-Kura-Asymptotic_gnn_std,tab:PG-Kura-Transient_gnn_std] report the performances across different model initializations.

7 Supplementary material on PG-Real experiments↩︎

The PG-Real dataset consists of synthetic topologies with validated dynamical properties, comprising an ensemble of networks with 70 to 80 nodes and individual networks of 128, 256, 512, and 1024 nodes. The dynamics are modeled using universal inverter-based resource models that accurately capture real power grid behavior. Stability metrics include transient stability measured by fault ride-through probability and asymptotic stability measured by the probability of returning to narrow frequency and voltage bounds.

The correlation between network measures and dynamical properties shows considerable variation ([tab:NF all correlations,tab:NF abs correlations,tab:NF var correlations]. Overall correlation values are low, with the notable exception of dynamic node types and transient stability. The heterogeneous dynamical properties appear more directly influential than topological characteristics. For transient stability prediction, GNNs demonstrate excellent performance (23), while NetSciML lagg significantly. Predicting asymptotic stability proves much more challenging, achieving only modest \(R^2\) values of \(0.33\) and \(0.16\) on the 70 to 80 node ensemble (5). However, GNNs maintained relatively consistent performance when scaling to larger grids, whereas NetSciML fail on substantially on larger networks. This pattern aligns with observations from the PG-Kura experiments. 24 reports the performance depending on the number of training samples.

Details on analyzing realistic power grids with inverters↩︎

Building on the work of [39], we perform simulations on larger power grids to evaluate the generalization capabilities of our models across varying grid sizes. We follow the same algorithms as outlined in the original paper, with one key modification during the grid generation process. Specifically, for larger grids, it was challenging to identify a valid operating point. To address this, we initially treated the lines as static to obtain a feasible operating point before reintroducing the dynamic aspects of the actual power grids.

7.0.0.1 Training and model details

The GNN used is a two-layer DBGNN, with each layer comprising 30 internal steps. The number of hidden channels is set to 120 for both nodes and edges. In total, the model contains 147,961 parameters. For the transient task, node and edge dropout are set to 0.014877 and 0.002454, respectively. For the asymptotic task, dropout is set to 0.1 at both node and edge level.

One-cycle learning-rate schedulers are used with the following settings: initial division factor \(413{,}759\) (transient) and \(4{,}365{,}873\) (asymptotic); final division factor \(1{,}241{,}798{,}872{,}312.98\) (transient) and \(3{,}680{,}955{,}167{,}359.520508\) (asymptotic); and maximum learning rate \(863.4449718452441\) (transient) and \(851.067199\) (asymptotic).

[tab:PG-Real-Asymptotic_gnn_std,tab:PG-Real-Transient_gnn_std] show mean performance with standard deviation over model initializations.

8 Supplementary material on ML methods↩︎

Graph Neural Network architectures↩︎

Different types of GNN have been introduced, some of which are detailed below. For the composition of the GNN models, we closely followed [25]. Many GNN layers build on the Graph Convolution Network (GCN) introduced by [61]: \[H = \sigma(\overline{A} X \Theta),\] where \(H\) is the output of a layer, \(\sigma\) denotes the activation function, using the input features \(X\), the matrix \(\Theta\) containing the learnable weights and a slightly modified and re-normalized adjacency matrix \(\overline{A}\). To increase the considered region and to consider neighbors at further distance, multiple GCN layers can be applied consecutively.

To increase the potential region per layer, [53] use multiple exponents \(i\) of \(\tilde{A}\) within one layer according to the following scheme: \[H = \sum_{z=0}^Z D^{-\frac{1}{2}} A^z D^{-\frac{1}{2}} X \Theta_z.\]

In case of DBGNN [54], both node and edge features are updated simultaneously using the following idea for the updates: \[\begin{align} \mathbf{h}_v^{(l+1)} = \sigma \left( \sum_{e \in \mathcal{E}(v)} \mathbf{W}_1 \mathbf{h}_e^{(l)} + \mathbf{W}_2 \mathbf{h}_v^{(l)} \right), \\ \mathbf{h}_e^{(l+1)} = \sigma \left( \sum_{v \in \mathcal{V}(e)} \mathbf{W}_3 \mathbf{h}_v^{(l)} + \mathbf{W}_4 \mathbf{h}_e^{(l)} \right), \end{align}\] where \(\mathbf{h}_v^{(l)}\) and \(\mathbf{h}_e^{(l)}\) are the node and edge features at layer \(l\), \(\mathcal{E}(v)\) is the set of edges incident to node \(v\), \(\mathcal{V}(e)\) is the set of nodes incident to edge \(e\), and \(\mathbf{W}_i\) are learnable weight matrices.

Best subset selection using orthogonal matching pursuit↩︎

To analyze the performance of small models with just a few independent variables, we perform best subset selections for the linear models. Given input features \(X_1,\ldots, X_p \in \mathbb{R}^n\) and target variables \(y \in \mathbb{R}^n\), the best subset selection problem is to find parameters \(\beta \in \mathbb{R}^n\), s.t., \[\min_\beta \|y - (X\beta)\|^2_2 - \| \beta \|_0, \label{secBestSubsetSelection}\tag{4}\] with \(\|.\|_0\) denoting the number of non-zero coefficients of the argument (sometimes called \(\ell_0\)-‘norm’). This is just the least squares formulation of a linear model, with an added regularization term. As an algorithm for the best subset selection problem we use orthogonal matching pursuit (OMP) [62], which approximates the subset of the network measures , that best fits the variance in the data. If the reduced model captures the main information and achieves good performance, it allows for a more intuitive understanding of the model predictions. We approximate the subsets 3 predictors and denote the corresponding models as OMP3.

References↩︎