Oscillatory spreading and inertia in power grids

New sources of dynamical instability can arise as more variable renewable generators (VRGs) are used in power systems. Unlike synchronous generators (SGs, e.g., coal, natural gas, nuclear), VRGs connect to a power system via a power electronic device called an inverter, which results in different dynamical interactions with the system. Specifically, SGs provide inertial support—an automatic and inherent SG response to changes in the power balance—which is a critical stability mechanism that VRGs often lack. As inertial support decreases, changes in the power balance become larger and faster, which can overwhelm control devices and decrease the time available for operator intervention.¹ These changes can lead to further instabilities in the power system, such as cascading failures or the loss of system synchronization.

To address these issues, a large body of work has emerged to devise new inverters that provide inertial support,² develop new techniques to optimize inertial support,^3–5 and formulate new operational policies to ensure sufficient inertial support.⁶ Much of this previous work has focused on larger perturbations, such as the loss of a node, but small perturbations, such as those due to resource variations of VRGs, can be equally damaging.⁷ VRG power fluctuations differ from those of SGs or loads^8,9—larger fluctuations tend to be more likely, and inertia influences how far away from the VRG these fluctuations can be observed.¹⁰ Specifically, it has been shown that perturbations decay more quickly when there is more inertia in the power system,^11–13 but the effects of heterogeneous inertia, which have been shown to significantly alter the spread of power perturbations,¹⁴ are less understood. In general, power perturbations on networks that are more tree-like and less clustered tend to decay more quickly than in networks that are more meshed,^15,16 but the impact of where these sub-structures are located in a network, and how this changes the spread of power fluctuations, has not been considered. Finally, much of this prior work has defined the spread of a perturbation by relating the response of individual nodes and network distance. However, the relationship between the collective response of a power system and the location of the perturbed node in the network has practical applications and is underexplored.

The goal of this work is to evaluate the impact of the size of an inertia value at a particular network location on deviations of a power system from its equilibrium. We focus on small oscillations that result from perturbations that are representative of VRG power fluctuations.⁹ Furthermore, the inertia in the power system can change the magnitude and location of these small oscillations. To study these effects, we perturb a node and analyze the impact of different inertia values at that location. To understand the impact of a network structure, we repeat this experiment at numerous locations across the network. We model two different perturbations: a frequency impulse that represents a sudden and short duration power loss (e.g., from a fault) and a stochastic perturbation that is representative of the second-to-second power fluctuations of a renewable resource.

Unlike other work that has focused on characterizing perturbation arrival times and distances,^11–13 we instead focus on the collective deviation of all nodes from their equilibrium states. To achieve this, we measure the power-system response to these perturbations using the concept of localization. Specifically, we extend an existing localization metric, the inverse participation ratio, the IPR, in a novel way to quantify the transition from localized (involving only a few nodes) to delocalized (involving many nodes) system responses and relate this transition to the perturbed node’s inertia parameter. We then study the relationship between this transition and the network location of the perturbed node by introducing a new measure, the nodal mode angle (NMA), that captures the participation of a node in the total system dynamics and compares this among the nodes in the network. We compare this metric across a set of small networks that differ by only one edge and show that it effectively captures differences in the dynamics of these networks. Finally, we illustrate the practical utility of this metric by studying the two perturbation types for a large canonical power-grid test case, showing that it can accurately describe the differences in the delocalization transition at different network locations.

Understanding the interplay between network structure and power system inertia aids in the implementation of solutions for power-system frequency control. The ability to quickly assess the likely oscillatory behavior of a power system given a set of parameters is useful and important, especially as sources of inertia become more diverse and more coordination among power-system stakeholders is required.

We use the structure-preserving (SP) model to simulate the dynamics of a power system.¹⁷ Our focus is on the impact of generator inertia on power oscillations, with the ultimate goal of helping with the implementation and placement of control devices, and we do not directly model VRG inverters; therefore, higher-order models that consider control devices and voltage dynamics are not necessary here. Equations (1a)–(1c) define the dynamics of the nodes in a power system as stated in the SP model, where $n$ is the number of nodes in the network, the collection of which we refer to as the system,

There are $g=1,2,…,ng<n$ generator nodes, where Eqs. (1a) and (1b) define the circular motion of the rotors, while load/intermediary node dynamics are defined by Eq. (1c). $θg$ is the angle of rotation for the generator rotor, and $δi$ is the voltage phase angle of a node. It is common to designate one generator’s rotational angle as the reference from which to measure all other generator angles. Here, we assume that this is the first generator (i.e., $θ1=0$⁠). When power production matches power consumption, $∑i=1nPi=0$⁠, the system reaches a steady state (i.e., without an external force, the state variables do not change over time), and all generators have the same angular frequency (60 or 50 Hz). $Δωg$ is a generator’s angular frequency deviation as measured from the steady-state value (i.e., $Δωg=0$ means that the generator is rotating at the steady-state frequency). Every node has a damping constant $Di$⁠, which defines how quickly oscillations at that node dissipate. The susceptance of an edge that connects node $i$ to node $j$ is $Bij$⁠, and physically, it describes how easy it is for electric current to flow across the line. Finally, every generator has an inertia constant $Mg$⁠, the parameter of interest in this work, that describes the rotor’s resistance to change in its angular frequency. Although VRG inverters implement inertia as a virtual parameter via control mechanisms,¹⁸ modeling inertia as a physical parameter in this manner is a valid assumption for power systems today because it is still the dominant source of inertial support.¹⁹ 

We perform our experiments on a representative power-system test case: the 500-node network from Birchfield et al.,²⁰ constructed such that it accurately reflects the transmission system located in South Carolina, and parameter values from Xu et al.²¹ as defaults, unless stated otherwise. Figure 1 depicts this network, where green squares are generators, pink triangles are loads, and gray circles are intermediary nodes.

We study the system response, determined by Eqs. (1a)–(1c), to two realistic power perturbation types, adding a term to the equation of motion of the perturbed generator $gΔ$⁠,

We set the magnitude of both perturbations such that other power system instabilities—e.g., cascading failures or generator disconnections—are unlikely to occur. The first perturbation is an impulse modeled by Eq. (2) and represents a sudden change in the power output of $gΔ$ (e.g., from a fault),

The magnitude of the impulse perturbation is $<1%$ of the total power generation/consumption of the power system, well within the range that is considered to be small.²² The second perturbation is stochastic, which models the second-to-second variation in the power output of a VRG, specifically when the renewable resource is wind. We use the model and parameters from Schmietendorf et al.⁸ to generate a representative time series for the stochastic perturbation, shown in Fig. 2. The model and relevant parameters are shown in Eq. (3), where $Γ(t)$ is Gaussian white noise with $μ=0$ and $σ=1$⁠,

The term $2$ reflects the strength of the randomness injected into the power fluctuations, where we have chosen a value that corresponds with strong intermittency.⁸ Additionally, the term $−β$ in Eq. (3b) is a damping term that impacts the variance of the variable $u$ and introduces a one second time correlation in the noise, which has been shown to occur in wind power fluctuations.²³ 

We perform a set of numerical experiments to evaluate the dynamical implications of these two perturbations delivered to a generator node with a specific inertia value. For a given simulation, one generator is perturbed, and we denote this generator by $gΔ$⁠. In order to explore the effects of network structures, we run simulations at each generator location defined in the power-system test case. To aid with clarity, we define the set of simulations associated with a particular generator location as a scenario. For a given simulation in a scenario, we measure the (de)localization of the system response (discussed further in Sec. III) due to each of the perturbations. We repeat this simulation using different inertia values for $gΔ$ in the range of $MgΔ∈$ [0.01–20.0] with a step of $0.01$⁠, where the minimum value is the floor of the default inertia values and the maximum value was chosen to be double the maximum inertia value typically observed in power systems.²¹ Therefore, a scenario consists of a chosen $gΔ$ and 4000 separate simulations, 2000 for each perturbation type. We use the Julia package DifferentialEquations.jl²⁴ to simulate Eqs. (1a)–(1c). Specifically, we use the adaptive time-step KenCarp4() solver with an absolute tolerance of $1×10−9$ for the impulse perturbation and the fixed time-step TangXiaoSROCK2() solver with $dt=1×10−3$ for the stochastic perturbation.²⁵ 

For each simulation described in Sec. II, we quantify the impact of a perturbation on the power system by measuring the deviations of the state variables. Let

be the state vector of the system described by Eqs. (1a)–(1c) at time $t$⁠. We define the system response to a power perturbation by the normalized vector $x^=[⟨Δω1⟩,⟨Δω2⟩,…,⟨δn−1⟩,⟨δn⟩]$⁠, where $⟨.⟩$ denotes the root-mean square of the variable as measured from its steady-state value over all time steps. Therefore, for example, $⟨Δω1⟩=1tf∑t=0tf(Δω1(t)−Δω1∗)2$⁠. $x^$ captures the magnitude of state variable deviations over time and the time each state variable spends away from its steady-state value. By normalizing $x$⁠, the deviations of the state variables are compared to that of the perturbed node (which will generally have the largest entry in $x$⁠). Previous work^2,3 has found utility in the metric $x^$⁠, especially for small perturbations where typical frequency control metrics such as the maximum frequency deviation or the maximum rate of change of frequency are incomplete descriptions of the collective impact of such power fluctuations.¹ 

We now discuss our metric for system (de)localization. The inverse participation ratio (IPR), which quantifies the density of a vector, was first introduced in Anderson²⁶ to quantify the density of electron states on lattices of various materials. Since then, it has been used increasingly in network-science applications. For example, Pastor-Satorras and Castellano^27,28 relate the (de)localization of the primary eigenvector of the adjacency matrix with epidemic spreading on a network. In related work, Torres-Sánchez et al.¹⁶ show that increased randomness in the small-world network model increases the IPR of the Fiedler vector, which indicates that small perturbations are likely to remain more localized.

The IPR of a normalized vector $z$ with $n$ entries is defined by Eq. (4),

When $z$ is totally localized (i.e., only one entry is non-zero), $IPR(z)=1$⁠; when the vector is totally delocalized, all entries are non-zero, have equal magnitude, and $IPR(z)=1/n$⁠. We define system (de)localization using $IPR(x^)$⁠. This differs from past work^11–13 that has used the IPR on network matrices or where localization was quantified by arrival times and network distances (i.e., localization was defined by the spread across the network). Here, $IPR(x^)$ quantifies the number of state variables that fluctuate and the magnitude of those fluctuations. This is effective for power-system applications because it is indicative of the oscillatory behavior of the system (i.e., local or wide-spread fluctuations) without requiring any knowledge of the network structure. Additionally, it has potential in other network applications where spreading phenomena are important to understand (e.g., brain networks).

As an example, Fig. 3(a) shows $IPR(x^)$ for both perturbation types as a function of $MgΔ$ where the node $gΔ$ is identified by an A in Fig. 1. The results illustrate some important and representative features in our experiments, beginning with a clear relationship between system (de)localization and $MgΔ$⁠. Specifically, as $MgΔ$ increases, the system response $x^$ transitions from mostly localized, i.e., large values of $IPR(x^)$⁠, to delocalized, indicated by small values of $IPR(x^)$⁠. Furthermore, there are obvious peaks and valleys in the plot, showing that at some inertia values, $x^$ is more localized (peaks) or delocalized (valleys) as compared to nearby inertia values. In other words, the changes in $IPR(x^)$ correspond to switches in which node(s) have a resonant response to the perturbed node. Clearly, $IPR(x^)$ is sensitive to $MgΔ$⁠, especially when it is small.

As shown in Fig. 3(a), the system response to the impulse perturbation is always more localized than to the stochastic, given the same inertia value. This is the case for every scenario, not just this example. This difference occurs because the frequency of the system oscillations due to the stochastic perturbation is smaller than for the impulse perturbation. A smaller frequency is usually associated with a more delocalized system response.¹⁶ Furthermore, the $IPR(x^)$ curve associated with the stochastic perturbation has more peaks and valleys than the impulse. These additional peaks are due to resonance phenomena induced by the frequency content of the stochastic perturbation that is absent from the impulse perturbation.^29,30 The sensitivity of $IPR(x^)$ to the inertia value of the perturbed generator clearly depends on the dynamics of the perturbation, which highlights the importance of studying both of these perturbation types.

We introduce a secondary metric on $IPR(x^)$ that accounts for the changes in $IPR(x^)$ due to changes in $MgΔ$⁠. Consider Fig. 3(b), which compares $IPR(x^)$ of two different scenarios, marked as B and C in Fig. 1, for the impulse perturbation. $IPR(x^)$ overlaps at several values for both scenarios when $MgΔ$ is small, but scenario B has significantly more peaks and valleys. Physically speaking, this means that there are more nodes that have significant and isolated frequency responses to perturbations at $gΔ$ for scenario B. As the inertia value of the perturbed node increases, scenario C exhibits a decrease in $IPR(x^)$ at a smaller inertia value than scenario B, but $IPR(x^)$ of scenario B eventually becomes smaller around $MgΔ≥5$⁠. This means that at larger inertia values, more nodes deviate from the synchronized frequency for perturbations at scenario B, and those responses are smaller in magnitude as compared to those of scenario C.

These differences in the (de)localization transition are also important to capture. Our approach to this is to integrate $IPR(x^)$ from $M=0.01$ to $M=20$ (double the largest default inertia value) as shown in Eq. (5), which contextualizes the change in system (de)localization in terms of the change in $MgΔ$⁠. We refer to this quantity as the (de)localization transition integral (DTI) throughout the rest of this paper,

Whether or not a large or small DTI is beneficial or detrimental to system stability is dependent on the particular system and its operating requirements. A larger DTI indicates that the system response is typically localized for most values of $MgΔ$⁠. This has important possible applications, such as for system planning of inertial support.⁶ Nodes with larger DTI could be a source of flexibility because their inertia value does not drastically change the system response to perturbations (i.e., guaranteeing stability is less dependent on that particular inertia value). Furthermore, a localized system response could be beneficial if it is necessary to isolate the oscillating system. In this case, less of the network would need to be disconnected in order to prevent a perturbation from spreading further because there are fewer nodes involved. However, localized responses could also cause devices to go outside their valid operating ranges, which presents a stability risk (e.g., causing cascading failures). To avoid this, these nodes may require extra damping, via a power system stabilizer (PSS), because the majority of the network nodes do not help absorb the perturbation.

On the other hand, a smaller DTI suggests that $gΔ$ is particularly influential on the spread of power fluctuations because many nodes tend to respond to perturbations at that location. The advantage in this situation is that the power system’s stability benefits from larger inertia values at this location because it helps prevent oscillations from spreading further through the system. The disadvantage is that malfunctions and failures at this location could impact more nodes, and this would have a larger impact on system stability (e.g., loss of synchronism or further outages). Therefore, we do not make a priori assertions about which values of the DTI are better for system stability. Alternatively, we seek to identify the dependence of the DTI on the network structure. A better understanding of this relationship can help with power system planning and operation, such as the placement and sizing of inertia sources or power-system stabilizers, especially as more VRGs are integrated.

As discussed in Sec. III, quantifying the dependence of system (de)localization on the inertia value of a node has important possible applications. However, computing the DTI requires numerous simulations, which is not desirable or feasible in many situations (e.g., hourly planning for very large systems). To address this, we develop a heuristic to relate the DTI with the network structure via the eigenvalues and eigenvectors of the network Laplacian.

A number of studies exist that have shown the utility of measures constructed from the spectral characteristics of this matrix in power system applications. Zhang et al.^30,31 showed that the entries of the Laplacian eigenvectors can be used to predict the maximum response (i.e., largest frequency deviation) of nodes in a power system due to stochastic fluctuations at a particular location. Similarly, Torres-Sánchez et al.¹⁶ use a node’s linear response to perturbations, defined as an average across its entries in the Laplacian eigenvectors, to devise a control metric for power system stabilizers. Our metric shares similarities with these studies, but rather than predict individual node responses or arrival times, we focus on comparing the metric among all of the nodes in the network to better understand the role of network structures in the system dynamics.

The Laplacian matrix is defined as $L=D−A$⁠, where $D$ contains the weighted degree of each node on the diagonal and $A$ is the weighted adjacency matrix. The Jacobian matrix associated with the system of Eqs. (1a)–(1c) takes the form of a Laplacian matrix.¹⁶ This means that the eigenvalues and eigenvectors of the network Laplacian are descriptors of independent patterns of motion of the system dynamics. Therefore, for small perturbations from the steady state, two nodes will have equivalent oscillations if they have equivalent entries in every eigenvector. We use this observation to assess the dynamical similarity of two nodes.

To quantify the above ideas, we introduce a metric we call the nodal mode angle (NMA). The NMA is computed in the following way. Define the sensitivity vector of a node $i$ to be

In other words, $vki$ is the $i$th entry of the $k$th eigenvector where eigenvalues maintain the relationship $λ2≤λ3≤⋯≤λn−1≤λn$ (excluding $λ1=0$⁠). We define the NMA between nodes $i$ and $j$ by Eq. (7),

where $||.||$ is the $L2$ norm. In other words, $αi,j$ is the angle between the sensitivity vectors of node $i$ and node $j$⁠. We use the notation $αi$ to indicate the set of all $αij$ for node $i$⁠.

We explore how $αij$ and $IPR(x^)$ change with network structures for five small graphs: motifs that occur commonly in power networks. The star graph in Fig. 4(a), for instance, occurs 76 512 times in the network of Fig. 1; the motifs in b, c, d, and e occur 776, 11 168, 658, and 516 times, respectively.

The graph topologies in Figs. 4(a) and 4(b) differ by one edge and similarly for those in Figs. 4(c)–4(e). The heatmaps in the second row indicate the value of the NMA for all pairs of nodes within a graph, where pink indicates $αij<π/2$ and green indicates $αij>π/2$⁠. Finally, we simulate the node dynamics defined by Eqs. (1a) and (1b) for an impulse and stochastic perturbation, separately applied to each node in each graph. We then measure $IPR(x^)$⁠, as shown in the scatterplots in the two bottom rows of Fig. 4, where the horizontal axis indicates the perturbed node from the associated graph. As discussed in Sec. III, heterogeneous parameters have a significant influence on $IPR(x^)$⁠. To isolate the effects of network structures, we use homogeneous parameters for these small networks, as listed in the caption of Fig. 4.

The NMAs of a network illustrate the tendency of nodes to behave similarly in the total system dynamics, which influences the spread of a perturbation. Consider the star graph of Fig. 4(a), which is representative of numerous generators or loads that connect to a power system at a common point. The NMAs of the pendant nodes suggest that they are most likely to oscillate in anti-phase. On the other hand, node 5 has a sensitivity vector that is nearly orthogonal to that of the pendant nodes, suggesting that perturbations at this location are likely to remain localized. This is confirmed by the $IPR(x^)$ calculations for both perturbation types. Adding an edge between two pendant nodes, as in Fig. 4(b), changes the NMAs and the spread of frequency fluctuations for those nodes. However, the change in $IPR(x^)$ is dependent on the perturbation type. For impulse perturbations at nodes 3 and 4, the system response becomes more delocalized with an edge addition. This is because the impulse perturbation causes nodes 1 and 2 to oscillate anti-phase with nodes 3 and 4, and the edge addition strengthened this trend, as indicated by the NMA. On the other hand, stochastic perturbations at nodes 3 and 4 induce a more localized system response. This is because nodes 3 and 4 are less likely to oscillate anti-phase with the edge addition, according to the NMA, and therefore, the perturbation decays more quickly between them. These two examples show that the NMA can illuminate differences in network structures that are influential on the spread of a perturbation through a power system, even when the dynamics of the perturbation have different effects on its spread.

Edge additions change the spread of perturbations in a power system, and the location of the new edge is also important to consider. In contrast to the star graph, Fig. 4(c) shows a chain graph structure that occurs in power systems to connect geographically distant locations together. Here, the NMA strongly depends on whether two nodes are closer to one another than they are to node 3 (e.g., nodes 1 and 2 vs nodes 1 and 4). When a new edge is added between nodes on opposite sides of the chain, as in Fig. 4(d), most NMAs move closer to $π/2$⁠, and $IPR(x^)$ increases for all nodes except node 3. However, adding an edge between two nodes on the same side of the chain, as in Fig. 4(e), creates two clusters in the heatmap of the NMA and perturbations at node 3 spread less. In other words, the location of the triangle has different effects on the NMA and $IPR(x^)$⁠.

These results illustrate that the NMA and $IPR(x^)$ are sensitive to network structures, even when networks only differ by one edge. It has been shown in other works that triangles typically decrease the spread of power fluctuations.¹⁶ However, here, we show that the local network structure around the triangle influences this trend. For example, creating a triangle in the chain graph decreased the spread of perturbations for the nodes connected by the new edge, but adding a triangle to the star graph actually increased the spread of impulse perturbations. Furthermore, these analyses provide evidence that the NMA and the spread of frequency oscillations through a power system are related. In particular, when $αi$ has many entries near $π/2$⁠, perturbations at that node tend to remain localized, but when $αi$ has a larger variance, power fluctuations tend to spread more. To further test these hypotheses, we next study how the distribution of $αi$ relates to system (de)localization for the 500-node test case shown in Fig. 1.

To relate the impact of different inertia values and network locations with the spread of power fluctuations, we now analyze the relationship between the NMA and the Delocalization Transition Integral (DTI) introduced in Sec. III. Recall that the DTI associates the inertia value of a perturbed node with the resulting system (de)localization. Specifically, a larger DTI signifies that the system response is more localized and less sensitive to $MgΔ$⁠, while a smaller DTI indicates that the system response is more delocalized and more sensitive to $MgΔ$⁠. We use the NMA to corroborate those assertions and understand how these differences in the DTI relate to the different locations of $gΔ$ and perturbation types.

As we illustrated in Sec. IV, the set of angles $αi$ is associated with the spread of perturbations that occur at node $i$⁠. We test how the distribution of $αi$ relates to the DTI by computing the dispersion index of $αi$⁠: the variance divided by the mean, $I(αi)=σ(αi)/μ(αi)$⁠. This allows us to compare the variation of the NMA distributions while also accounting for their means. In other words, $I(αi)$ accounts for the situation where $σ(αi)=σ(αj)$⁠, but $μ(αi)≠μ(αj)$⁠.

Figure 5 shows the DTI for every scenario in the representative power-system test case of Fig. 1 as a function of $I(αi)$⁠. Pink and green points indicate the impulse and stochastic perturbations, respectively. Additionally, we computed a best fit curve for each perturbation type using the method of least squares, the equations of which are shown near the curves. We chose to fit a power-law function because $⟨ΔωgΔ⟩$⁠, which influences the value of the DTI the most, has a $1/MgΔ$ relationship with inertia. The further that the fitted curves are from this relationship, the more impactful the network structure is on the value of the DTI.

These results show a clear dependence between the DTI and $I(αi)$⁠. Nodes with larger values of $I(αi)$ tend to induce a delocalized system response, while nodes with smaller $I(αi)$ tend to induce a more localized one. Furthermore, because the DTI also describes the sensitivity of system (de)localization to the inertia value of the perturbed node, larger values of $I(αi)$ suggest that the system dynamics are more sensitive to the inertia value at that location.

Both perturbation types elicit a similar functional relationship between the DTI and $I(αi)$⁠. However, the impulse perturbation has a larger margin of error in the fitted parameters. This could suggest that the individual dynamics of the perturbed generator, which are determined by its damping and inertia parameters and are not captured by $I(αi)$⁠, have more of an impact on the spread of the impulse perturbation. On the other hand, the collective dynamics of the network nodes, which is determined by the network structure and is weighted more heavily in the sensitivity vector of a node, more strongly influence the spread of the stochastic perturbation.

The relationship shown in Fig. 5 between the DTI and $I(αi)$ has important practical applications, particularly for stabilizing system oscillations and addressing inertial support. For example, Copp et al.³² showed that the effectiveness of energy storage devices in damping oscillations that occur due to a perturbation is sensitive to the location of the device with respect to the perturbed node. In this case, $I(αi)$ can provide important information on where to put a damping device, and this assessment can be made without performing any numerical simulations. Specifically, damping devices will be most effective at locations very near to the node when $I(αi)$ is small because the system response is most likely to be localized, regardless of the inertia value at that node. On the other hand, when $I(αi)$ is large, two damping devices could be necessary: one close to the node to handle the localized dynamics that occur for smaller inertia values and one farther away to address delocalized dynamics for larger inertia values.

This study shows that the NMA effectively illustrates the complex relationship between the system dynamics and inertia values, perturbation dynamics, and network locations. Furthermore, we showed that the dispersion index of the distribution of the NMA for each $gΔ$ location illuminates the relationship between all of these factors, which could help in many different areas of power-system planning and operation.

In this work, we studied the impact of a crucial dynamical parameter, the inertia value, on the spread of power fluctuations in a power system. Specifically, we analyzed a suite of numerical simulations of a power system’s response to perturbations at different nodes in the network under different inertia conditions. We studied system (de)localization, and its dependence on an inertia value, in a novel way using the inverse participation ratio, showing that smaller inertia values induce localized system dynamics and vice versa for larger inertia values. We also found that the transition between the localized and delocalized regimes reflects the sensitivity of system (de)localization to that particular inertia value. We introduced a second metric, the nodal mode angle (NMA), to relate these findings to the network structure and the location of the perturbed node. We studied a set of small networks that occur frequently in real power systems and showed that the NMA is sensitive to differences in these structures. Finally, we explored the relationship between the NMA and system (de)localization. We discussed several ways this relationship could be leveraged in power-system planning and operation, especially as these relate to inertial support.

Future work will further investigate the relationship between the NMA, network structure, and perturbation spreading. A meaningful next step is to consider the effects of multiple (possibly correlated) perturbations, as this will become more relevant with increased penetration of renewable energy sources. Additionally, applying our analyses with different inertia sources, such as grid-forming VRG inverters,³ could help further illuminate important considerations for future planning of inertial support. A thorough study of the relationship between the distribution of $αi$ and network properties (such as the number of short-cut edges, the degree distribution, or the occurrence of different network motifs) could lead to a better understanding of how different attributes of the $αi$ distribution arise (e.g., multiple peaks). To further understand the importance of network structures on the spread of power fluctuations in power systems, it would be useful to assess the descriptive differences between the NMA and distance-based measures such as the resistance distance. Finally, the NMA has many possible applications in power systems research that warrant further investigation such as the assessment of system vulnerability to cascading failures or the evaluation of a new infrastructure.

This work was authored in part by the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308 with the Alliance for Sustainable Energy, LLC, the Manager and Operator of the National Renewable Energy Laboratory. Funding is provided by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes. NREL is a national laboratory of the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, operated by the Alliance for Sustainable Energy, LLC.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors have no conflicts to disclose.