Fluctuation response patterns of network dynamics – An introduction

Abstract

Networked dynamical systems, i.e., systems of dynamical units coupled via nontrivial interaction topologies, constitute models of broad classes of complex systems, ranging from gene regulatory and metabolic circuits in our cells to pandemics spreading across continents. Most of such systems are driven by irregular and distributed fluctuating input signals from the environment. Yet how networked dynamical systems collectively respond to such fluctuations depends on the location and type of driving signal, the interaction topology and several other factors and remains largely unknown to date. As a key example, modern electric power grids are undergoing a rapid and systematic transformation towards more sustainable systems, signified by high penetrations of renewable energy sources. These in turn introduce significant fluctuations in power input and thereby pose immediate challenges to the stable operation of power grid systems. How power grid systems dynamically respond to fluctuating power feed-in as well as other temporal changes is critical for ensuring a reliable operation of power grids yet not well understood. In this work, we systematically introduce a linear response theory (LRT) for fluctuation-driven networked dynamical systems. The derivations presented not only provide approximate analytical descriptions of the dynamical responses of networks, but more importantly, also allow to extract key qualitative features about spatio-temporally distributed response patterns. Specifically, we provide a general formulation of a LRT for perturbed networked dynamical systems, explicate how dynamic network response patterns arise from the solution of the linearised response dynamics, and emphasise the role of LRT in predicting and comprehending power grid responses on different temporal and spatial scales and to various types of disturbances. Understanding such patterns from a general, mathematical perspective enables to estimate network responses quickly and intuitively, and to develop guiding principles for, e.g., power grid operation, control and design.

1 Introduction

Networked dynamical systems abound around and in us. From the circuits supporting metabolism and gene regulation in our cells to the neural networks in our brains, from the power grids that supply electric energy to most of our technical infrastructure to the internet that connects our computers, all of these systems are driven externally, often by irregular, time-dependent and spatially heterogeneous signals. How networked dynamical systems respond to such perturbations, driving signals or other types of time-dependent inputs is hardly understood to date. In this article, we offer a general introduction to the basic theory of analysing response features of networked dynamical systems by linear response theory (LRT) and focus on applications in the realm of power grids.

A reliable supply of electric power fundamentally underlies most aspects of modern society. As the shares of renewable electric energy supply grow and consumer dynamics is increasingly influenced by digital technologies, fluctuations impinge on power grids, making them intrinsically driven, non-equilibrium systems, with distributed and often non-stationary response dynamics [41]. If fluctuations become large, they may destabilise grid dynamics, affect normal operations of parts of the grid and cause cascades of failures or even total blackouts [38, 40, 9]. To predict, control or mitigate the fluctuating and distributed responses of such networks resulting from input (and output) power fluctuations, we need to understand their nature in network dynamical systems.

How can we systematically characterise fluctuating responses that are distributed across meshed networks? How can we predict their influences and at which nodes in a network is their impact most profound? LRT relates sufficiently small time-dependent driving signals to time-dependent responses. The theory approximates the resulting system dynamics near some operating point to first order in the strength of the driving signals. It has traditionally found applications across physics, chemistry and engineering [15, 4, 13, 29, 20, 25], and recently also in power grid models [7, 42, 14, 21, 33, 11, 37]. However, a framework of LRT uncovering spatiotemporal response patterns in systems with multi-dimensional dynamics of units that simultaneously interact via intricate network topologies is not yet well established.

In this article, we introduce a general formalism of LRT for networked dynamical systems and demonstrate its applications in stationary and non-stationary models of power grids. Specifically, in Section 2, we present the main ideas of the LRT for networks by first analysing the linear responses of networked dynamical systems with the most general settings (Section 2.1) and then boiling down step-by-step to a specific LRT which provides a direct link between the temporal and the spatial features of the dynamic network responses (Section 2.2). In Section 3, we demonstrate the LRT for networks by applying it to two models of power grids: a stationary model for the DC approximation of the AC power flow distributions in power grids (Section 3.1), and a non-stationary model, the (second-order) oscillator model, commonly used to describe the dynamics of the high-voltage AC power transmission networks (Section 3.2). The next two sections, Sections 4 and 5, focus on the extraction and the interpretation of spatiotemporal response patterns of power grids from the LRT, as well as the analytical techniques needed therein. In Section 4, we elaborate how distinctive steady-state response patterns emerge in three frequency regimes (Section 4.1), how the transient spreading pattern of a perturbation entangles with the underlying network topology (Section 4.2), and how LRT helps to estimate the long-term risk of individual nodes from external fluctuations (Section 4.3). In Section 5, we summarise the role of LRT in uncovering such patterns by comparing the analyses for patterns in the steady-state response versus the transient responses (Section 5.1), the patterns in the deterministic responses to a given perturbation versus the cumulative responses to a random signal (Section 5.2), and the small responses close to the base operation state versus the large responses further away (Section 5.3). In the last section (Section 6), we point out several directions for extending the present theory for a better understanding of the response dynamics of networked dynamical systems and for a safer and more reliable operation of future power grids.

2 Main ideas of the LRT for networks

In this section, we introduce the main ideas of the LRT framework. Section 2.1 formulates the LRT on a general model of networks of one-dimensional dynamical units. Section 2.2 highlights the linear operators, represented as matrices, arising in the LRT for specific network settings and how basic information of the underlying network system such as its topology and features defining the dynamically determined operating point enter those operators. In Section 2.2, we also derive the explicit linear responses of the general network model of one-dimensional units introduced in Section 2.1 and extend the results for networks of higher-order nodal dynamics.

2.1 General formulation of LRT for networked dynamical systems

We now illustrate the main idea underlying LRT by starting with a general setting of networked dynamical systems where each unit’s dynamics is one-dimensional. We consider a dynamical process involving N interacting variables, whose state is represented by a vector $\boldsymbol{x}=(x_1,\cdots,x_N)\in \mathbb{R}^N$ . The autonomous dynamics of the N-dimensional dynamical system is governed by

(2.1) \begin{equation} \dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x}),\end{equation}

where $\boldsymbol{f}:\mathbb{R}^N\rightarrow\mathbb{R}^N$ is a function that in general depends on the states $\boldsymbol{x}$ of all units and does not explicitly depend on time. Let us consider a system that exhibits a fixed point at $\boldsymbol{x}=\boldsymbol{x}^*$ where $\boldsymbol{f}(\boldsymbol{x}^*)=\textbf{{0}}$ , and is influenced by an external dynamic driving vector $\boldsymbol{D}(t)\in\mathbb{R}^N$ at the fixed point $\boldsymbol{x}^*$ . We investigate how the system dynamically responds to $\boldsymbol{D}(t)$ , i.e. how the system’s deviation $\boldsymbol{X}(t):=\boldsymbol{x}(t)-\boldsymbol{x}^*$ from the fixed point evolves in time. In general, the dynamics of the system’s response $\boldsymbol{X}$ follows

(2.2) \begin{equation} \dot{\boldsymbol{X}}=\boldsymbol{f}(\boldsymbol{x}^*+\boldsymbol{X})+\boldsymbol{D}(t),\end{equation}

where both functions $\boldsymbol{f}$ and $\boldsymbol{D}$ can be highly nonlinear so that an exact analytical solution of the system’s response $\boldsymbol{X}$ typically does not exist or is unknown.

Important information about the response dynamics (2.2) is given by the stability operator at the fixed point obtained from the linearisation of the function $\boldsymbol{f}$ at $\boldsymbol{x}=\boldsymbol{x}^*$ , i.e. the Jacobian matrix $\mathcal{J}$ with $\mathcal{J}_{ij}:=\frac{\partial f_i}{\partial x_j}|_{\boldsymbol{x}=\boldsymbol{x}^*}$ . The system’s response dynamics thus approximately follows the linearised differential equation

(2.3) \begin{equation} \dot{\boldsymbol{X}}=\mathcal{J}\boldsymbol{X}+\boldsymbol{D}(t).\end{equation}

The signs of the real parts of the Jacobian eigenvalues $\boldsymbol{w}^{[\ell]}$ (with $\ell$ being the index) indicate whether the fixed point $\boldsymbol{x}^*$ is linearly stable $\left(\hbox{Re}\left[\boldsymbol{w}^{[\ell]}\right]<0\right)$ , unstable $\left(\hbox{Re}\left[\boldsymbol{w}^{[\ell]}\right]>0\right)$ or neutrally stable $\left(\hbox{Re}\left[\boldsymbol{w}^{[\ell]}\right]=0\right)$ in the eigendirection/eigenspace corresponding to the eigenvalue $\boldsymbol{w}^{[\ell]}$ . Below, we typically focus on the dynamics near stable fixed points, where all $\hbox{Re}\left[\boldsymbol{w}^{[\ell]}\right]<0$ (or marginally stable ones where one $\hbox{Re}\left[\boldsymbol{w}^{[\ell]}\right]=0$ if $\mathcal{J}$ is a Laplacian operator) to justify the assumption that the solution $\boldsymbol{X}(t)$ of the linearised equation (2.3) reasonably well approximates the full solution of (2.2).

Let us now focus on systems with pairwise interactions between the N units such that the function $f_i(\boldsymbol{x})$ controlling the time evolution of unit i can be written as

(2.4) \begin{equation} f_i(\boldsymbol{x})=h_i(x_i)+\sum_{j=1; \ j\neq i}^N K_{ij}g_{ij}(x_i,x_j),\end{equation}

where $h_i$ denotes the intrinsic dynamics depending on the variable $x_i$ itself and the coupling term $K_{ij}g_{ij}(x_i,x_j)$ represents the pairwise interaction between variable $x_i$ and $x_j$ with $i\neq j$ . Here, $K_{ij}\in\mathbb{R}^+_0$ is the non-negative coupling strength and $g_{ij}$ is the coupling function depending on the state of the pair of variables $(x_i, x_j)$ and $g_{ij}(x_i,x_j)\not \equiv 0$ .

Remark 2.1 Pairwise interactions (2.4) induce an interaction structure of the dynamical system (2.1) that can be represented by a graph G(V,E) where the set of vertices V consists of N variables $x_1,\cdots, x_N$ and the set of edges E consists of all node pairs with the pairwise coupling strength being nonzero, i.e. $E=\{(i,\, j)\in V^2 | \ (i\neq j) \land (K_{ij}\neq0)\}$ .

The linearisation of the pairwise interaction (2.4) yields the Jacobian matrix of the response dynamics (2.2) of the driven networked system

(2.5a) \begin{eqnarray} \mathcal{J}_{ij} &=& \frac{\partial}{\partial x_j}\left.\left({h}_i({x}_i)+\sum_{k\neq i} K_{ik}{g}_{ik}\left({x}_i,{x}_k\right)\right)\right|_{\boldsymbol{x}=\boldsymbol{x}^*}\quad\end{eqnarray}

(2.5b) \begin{eqnarray} \quad\qquad &=& \left\{\begin{array}{l@{\quad}l} \left.\dfrac{\mathrm{d} {h}_i}{\mathrm{d} {x}_i}\right|_{\boldsymbol{x}=\boldsymbol{x}^*}+\displaystyle\sum_{k\neq i} K_{ik}\left.\dfrac{\partial{g}_{ik}}{\partial {x}_i}\right|_{\boldsymbol{x}=\boldsymbol{x}^*} &\mathrm{for}\quad j=i \\[18pt] K_{ij}\left.\dfrac{\partial{g}_{ij}}{\partial {x}_j}\right|_{\boldsymbol{x}=\boldsymbol{x}^*} &\mathrm{for}\quad j\neq i. \end{array}\right.\end{eqnarray}

It is clear that the topology of the interaction network G explicitly enters the Jacobian matrix through the matrix of the coupling strength $K\in\mathbb{R}^{N\times N}$ with its ij-th element being defined as $K_{ij}$ , which is effectively a weighted adjacency matrix of graph G.

The combination of the linearised dynamics of a general N-dimensional networked dynamical system (2.1) near a fixed point (2.3) and the arising Jacobian matrix that explicitly depends on the system’s topology (2.5) due to pairwise interactions (2.4) provides a general form of LRT for networked dynamical systems. For a specific system with given forms of or constraints for the intrinsic nodal dynamics $h_i(x_i)$ , coupling strengths $K_{ij}$ and coupling function $g_{ij}(x_i,x_j)$ , the solution of the matrix equation (2.3), X(t), can be characterised by evaluating the specific spectral properties of the Jacobian matrix $\mathcal{J}$ . Thereby the LRT provides a powerful tool to describe the dynamic response of a networked system temporally and spatially at once: the solution X(t) explicitly depends on time, and at the same time its relation to the topology-dependent Jacobian matrix can be exploited to reveal the temporally and spatially distributed patterns of the network response.

2.2 Arising linear operators and network topology

As discussed in the last section, the dependence of the Jacobian matrix of network response dynamics on a weighted adjacency matrix of the underlying interaction network gives the first hint on how temporal and spatial features of dynamic network responses are intertwined. In this section, we further show that, under a few commonly satisfied conditions such as diffusive coupling between units, interesting results of network response dynamics emerge. Especially, another important graph-theoretical matrix, the Laplacian matrix, arises in the network response dynamics (2.3), providing us with powerful tools for characterising the spatiotemporal patterns in dynamic network responses.

Diffusive coupling is a very common type of coupling present in many physiological and chemical systems [12, 18, 27, 5, 32], in particular also appearing in the Kuramoto model [17] and its variations [10, 1]. A diffusive coupling function $\tilde{g}_{ij}$ mediating the interaction between a pair of nodes (i, j) is characterised by its dependence on the state difference $x_j-x_i$ of the node pair, i.e. $g_{ij}(x_i,x_j)=\tilde{g}_{ij}(x_j-x_i)$ in (2.5). For notational simplicity, we again denote the functions $\tilde{g}_{ij}$ just by $g_{ij}$ .

Proposition 2.1 (Stability of diffusively coupled networks: a special case) A networked dynamical system with evolution function (2.4) and interaction network G(V,E) is at least neutrally stable at its fixed point $\boldsymbol{x}=\boldsymbol{x}^*$ if a) the intrinsic nodal dynamics $h_i$ satisfies $\frac{\mathrm{d} h_i}{\mathrm{d} x_i}|_{\boldsymbol{x}=\boldsymbol{x}^*}\leqslant 0$ for all nodes $i\in V$ , b) the coupling function $g_{ij}$ is diffusive and satisfies $\frac{\mathrm{d} g_{ij}}{\mathrm{d}(x_j-x_i)}|_{\boldsymbol{x}=\boldsymbol{x}^*}\geqslant 0$ for all node pairs $(i, j)\in E$ .

Proof. The diffusive form of the coupling function $g_{ij}(x_j-x_i)$ yields a useful relation $\frac{\partial g_{ij}}{\partial x_j}=-\frac{\partial g_{ij}}{\partial x_i}=\frac{\mathrm{d} g_{ij}}{\mathrm{d}(x_j-x_i)}$ . With this particular symmetry, the Jacobian matrix (2.5) of the system at the fixed point $\boldsymbol{x}=\boldsymbol{x}^*$ takes the following form

(2.6) \begin{equation} \mathcal{J}_{ij}=\left\{\begin{array}{l@{\quad}l} -\beta_i-\displaystyle\sum_{k\neq i} K_{ik}\gamma_{ik} &\mathrm{for}\quad i=j \\ K_{ij}\gamma_{ij} &\mathrm{for}\quad i\neq j, \end{array}\right.\end{equation}

where $\beta_i:=-\left.\frac{\mathrm{d} h_i}{\mathrm{d} x_i}\right|_{\boldsymbol{x}=\boldsymbol{x}^*}$ and $\gamma_{ij}:=\left.\frac{\mathrm{d} g_{ij}}{\mathrm{d}(x_j-x_i)}\right|_{\boldsymbol{x}=\boldsymbol{x}^*}$ . Given that $\beta_i\geqslant 0$ , $\gamma_{ij}\geqslant 0$ and $K_{ij}\geqslant 0$ by definition, the Jacobian $\mathcal{J}$ is diagonally dominant:

(2.7) \begin{equation} \left|\mathcal{J}_{ii}\right|=\left|-\beta_i-\displaystyle\sum_{k\neq i} K_{ik}\gamma_{ik}\right| =\beta_i+\displaystyle\sum_{k\neq i} K_{ik}\gamma_{ik} \geqslant \displaystyle\sum_{k\neq i} K_{ik}\gamma_{ik} =\displaystyle\sum_{k\neq i} \left|K_{ik}\gamma_{ik}\right| =\displaystyle\sum_{j\neq i} \left|\mathcal{J}_{ij}\right|.\end{equation}

According to the Gershgorin circle theorem, every eigenvalue of the Jacobian matrix lies within at least one of the N Gershgorin discs $D_i(\mathcal{J}_{ii},r_i)=\left\lbrace z\in\mathbb{C} |\, \,\left|z-\mathcal{J}_{ii}\right|\leqslant r_i \right\rbrace$ with $r_i=\sum_{k\neq i}|\mathcal{J}_{ki}|$ in the complex plane. Relation (2.7) indicates that all N Gershgorin discs lie in the left half of the complex plane, i.e. in $\left\lbrace z\in \mathbb{C}\left|\,\hbox{Re}\left(z\right)\leqslant 0\right.\right\rbrace$ , because the centers of the discs $(\mathcal{J}_{ii},0)$ lie on the negative real axis and the radius of the discs $r_i=\sum_{k\neq i}|\mathcal{J}_{ki}|\leqslant |\mathcal{J}_{ii}|$ . The discs touch the imaginary axis from the half plane $\left\lbrace z\in \mathbb{C}\left|\,\hbox{Re}\left(z\right)\leqslant 0\right.\right\rbrace$ only if $\beta_i=0$ . Therefore, all Jacobian eigenvalues can only have non-positive real parts, and consequently, the networked dynamical system is at least neutrally stable at the fixed point $\boldsymbol{x}=\boldsymbol{x}^*$ .

Remark 2.2 (Homogeneous nodal dynamics and graph Laplacian) We remark that in case of identical intrinsic nodal dynamics at all nodes, a weighted Laplacian matrix $\mathcal{L}$ of the interaction graph explicitly enters the linearised response dynamics of the network (2.3). Assuming $\beta_i=\beta\in\mathbb{R}$ for all nodes i, we can express the Jacobian matrix as

(2.8) \begin{equation} \mathcal{J}=- \beta \boldsymbol{1} - \mathcal{L},\end{equation}

where the weighted graph Laplacian $\mathcal{L}$ is defined as

(2.9) \begin{equation} \mathcal{L}_{ij}:=\left\{\begin{array}{l@{\quad}l} \displaystyle\sum_{k\neq i} K_{ik}\gamma_{ik} &\mathrm{for}\quad i=j \\[6pt] -K_{ij}\gamma_{ij} &\mathrm{for}\quad i\neq j. \end{array}\right.\end{equation}

Here $K_{ij}\gamma_{ij}$ is considered a weight of edge (i, j), containing the coupling strength $K_{ij}$ and the linearised coupling function $\gamma_{ij}$ at the fixed point.

Remark 2.3 (Symmetry and linear responses in Laplacian eigenbasis) In general, the interaction network can be directed, meaning that for an edge (i, j) the coupling strength $K_{ij}$ and the derivative of the coupling function $\gamma_{ij}$ at the fixed point, i.e. the sensitivity of the diffusive coupling function $g_{ij}(x_j-x_i)$ to a change in the state difference $x_j-x_i$ , can differ from their counterparts $K_{ji}$ and $\gamma_{ji}$ for the edge (j,i) with the opposite direction. This asymmetry leads to an asymmetric weighted graph Laplacian (2.9). However, for undirected networks with symmetric strengths ( $K_{ij}=K_{ji}$ ) and symmetric sensitivities of coupling functions ( $\gamma_{ij}=\gamma_{ji}$ ), or more generally, a symmetric combination $K_{ij}\gamma_{ij}=K_{ji}\gamma_{ji}$ , the weighted graph Laplacian $\mathcal{L}$ is symmetric. Its eigenvectors thus form an orthogonal basis which allows us to solve for the linear network responses in (2.3) and (2.8) by expressing the response vector $\boldsymbol{X}(t)$ in terms of the eigenvalues and eigenvectors of the Laplacian.

We first analyse a system that is perturbed by only one sinusoidal signal with magnitude $\varepsilon>0$ and frequency $\omega>0$ at node k, i.e. $D^{(k)}_i(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$ , where $\delta_{ik}$ is the Kronecker delta with $\delta_{ik}=1$ if $i=k$ and $\delta_{ik}=0$ otherwise. We solve for the linear network response vector $\boldsymbol{X}^{(k)}(\omega,t)$ governed by

(2.10) \begin{equation} \dot{\boldsymbol{X}}^{(k)}=-\beta\boldsymbol{X}^{(k)}-\mathcal{L}\boldsymbol{X}^{(k)}+\boldsymbol{D}^{(k)}(t).\end{equation}

Expressing the response vector in the constant eigenbasis of Laplacian

(2.11) \begin{equation} \boldsymbol{X}^{(k)}(t)=\sum_{\ell=0}^{N-1}c^{[\ell]}(t)\boldsymbol{v}^{[\ell]}\end{equation}

and the time-dependent coefficients $c^{[\ell]}(t)$ and exploiting orthogonality of the Laplacian eigenvectors, we obtain the ordinary differential equations

(2.12) \begin{equation} \dot{c}^{[\ell]}=\left(-\beta-\lambda^{[\ell]}\right)c^{[\ell]}+\varepsilon v^{[\ell]}_k e^{\imath(\omega t+\varphi)}\end{equation}

for the coefficients. Here, the N Laplacian eigenvalues $\lambda^{[\ell]}$ and eigenvectors $\boldsymbol{v}^{[\ell]}$ are indexed according to $0=\lambda^{[0]}\leqslant \lambda^{[1]}\leqslant \cdots\leqslant \lambda^{[N-1]}$ . If the graph is connected, the zero eigenvalue is unique, such that all other eigenvalues are positive real numbers [3]. The solution of differential equation (2.12) for the coefficients is

(2.13) \begin{equation} c^{[\ell]}(\omega,t)=\left(\boldsymbol{X}_0^{(k)}\cdot\boldsymbol{v}^{[\ell]}\right)e^{\left(-\beta-\lambda^{[\ell]}\right)t}+\dfrac{\varepsilon v^{[\ell]}_k e^{\imath\varphi}}{\beta+\lambda^{[\ell]}+\imath\omega}\left(-e^{\left(-\beta-\lambda^{[\ell]}\right)t}+e^{\imath\omega t}\right),\end{equation}

where $\boldsymbol{X}_0^{(k)}$ denotes the initial response vector at $t=0$ given node k is perturbed. The linear network response is thus given by (2.11) and (2.13).

Remark 2.4 (LRT for distributed arbitrary perturbations) In case perturbations with arbitrary temporal structures are distributed across the network, that is, multiple elements of the perturbation vector $\boldsymbol{D}(t)$ are arbitrary time series, the linear network response can be obtained by summing up the single-signal single-frequency response (2.11 and 2.13) over all frequency components $\omega$ and all perturbation signals k by resorting to the linearity of the dynamics (2.3)

(2.14) \begin{equation} \boldsymbol{X}(t)=\sum_{k}\left(\boldsymbol{X}^{(k)}_{\omega=0}(t)+\sum_{\omega>0}\boldsymbol{X}^{(k)}(\omega,t)\,\mathrm{d}\omega\right).\end{equation}

Here, $\boldsymbol{X}^{(k)}_{\omega=0}(t)$ denotes the linear response to a constant ( $\omega=0$ ) perturbation $\varepsilon$ , which shares the same form as $\boldsymbol{X}^{(k)}(\omega,t)$ with $\omega=0$ for $\beta>0$ but not for $\beta=0$ . Moreover, the sum becomes an integral for continuously distributed frequencies, with $\boldsymbol{X}^{(k)}(\omega,t)$ replaced by the associated response density per unit frequency.

For systems where nodal damping vanishes ( $\beta=0$ ), the coefficient for the 0-th eigenmode $c^{[0]}$ diverges for a constant perturbation as $\omega \rightarrow 0$ , so it needs to be solved separately and takes the form of $\boldsymbol{X}_0^{(k)}\cdot\boldsymbol{v}^{[0]}+\varepsilon v_k^{[0]} t$ , inducing a linear drift and thus unbounded growth with time. Since an unbounded growth typically induces the approximation of the linear response to the real system to break down, this solution would become useless in practice due to larger errors between the approximation and the exact solution.

Remark 2.5 (LRT for higher-order nodal dynamics) In the above paragraphs, we discussed the main ideas of the LRT for networked dynamical systems with first-order nodal dynamics. For more general systems with second- or higher-order nodal dynamics, the straightforward relation (2.8) between the Jacobian matrix and a weighted graph Laplacian does not hold any more. Nevertheless, a symmetric weighted graph Laplacian still arises in the linearised response dynamics for diffusively coupled undirected networks with symmetric coupling strengths and symmetric sensitivities of coupling functions as discussed in Remark 2.3. If the higher-order time derivatives of the state variables has homogeneous coefficients for individual nodes, i.e. the response dynamics to a perturbation $\boldsymbol{D}(t)$ has the form of

(2.15) \begin{equation} \sum_d\kappa_dD_t^{d}\boldsymbol{X}=-\mathcal{L}\boldsymbol{X}+\boldsymbol{D}(t), \end{equation}

an explicit solution of the linear responses in the eigenbasis of $\mathcal{L}$ can still be obtained following the routine in Remark 2.3, if the corresponding ODEs for the time-dependent coefficients

(2.16) \begin{equation} \sum_d\kappa_dD_t^{d}{c}^{[\ell]}=-\lambda^{[\ell]}c^{[\ell]}+\boldsymbol{D}(t)\cdot\boldsymbol{v}^{[\ell]} \end{equation}

are solvable. Here, $\kappa_d\in\mathbb{R}$ are constant coefficients and we use Euler’s notation for derivatives, where $D_t^{d}x$ denotes the d-th time derivative of variable x.

In summary, if a networked dynamical system consisting of N units

1. (1) has a fixed point at $\boldsymbol{x}=\boldsymbol{x}^*$ ,

2. (2) in the neighbourhood of $\boldsymbol{x}=\boldsymbol{x}^*$ the intrinsic nodal dynamics $h_i(x_i)$ gives homogeneous non-positive feedback to the respective nodes, i.e. $\beta_i:=-\left.\frac{\mathrm{d} h_i}{\mathrm{d} x_i}\right|_{\boldsymbol{x}=\boldsymbol{x}^*}=\beta\geqslant 0$ for all i, and

3. (3) the coupling function $g_{ij}(x_j-x_i)$ is diffusive and the coupling term’s sensitivity to small changes at $\boldsymbol{x}=\boldsymbol{x}^*$ is symmetric and non-negative, i.e. $K_{ij}=K_{ji}$ and $\gamma_{ij}:=\left.\frac{\mathrm{d} g_{ij}}{\mathrm{d}(x_j-x_i)}\right|_{\boldsymbol{x}=\boldsymbol{x}^*}=\gamma_{ji}\geqslant 0$ for all edges (i, j),

then i) the dynamical system is at least neutrally stable (Proposition 2.1) and ii) a symmetric weighted graph Laplacian (2.9) arises in the network response dynamics (2.10) and enables the expression of linear network responses in the Laplacian eigenbasis (2.11) and (2.13) (Remark 2.2 and Remark 2.3). As we will show in Section 4, the explicit dependence of the linear network response on Laplacian eigenvalues and eigenvectors provides a powerful tool to reveal how the dynamic responses spatially distributed across the network.

3 LRT for power grid models

We now discuss how the LRT for general networked dynamical systems introduced in Section 2 applies to stationary and non-stationary models of power grids and helps reveal static and dynamic responses of power grid systems to external perturbations.

3.1 LRT for the DC power flow model

We first demonstrate how LRT works in a minimal model, the DC power flow model and how it helps to compute the systemic stationary response of a power transmission network to perturbations in power injections and withdrawals.

For common AC power grids, the full power flow analysis poses several challenges such as possible difficulties in finding a solution in ill-conditioned cases and the existence of multiple solutions due to the inherent nonlinearities [23]. By linearising the AC power flow equations at an operation point, the DC power flow model ¹ provides a relatively simple and computational inexpensive way to compute the power flows.

In an AC power transmission grid, the total complex power flow $S_j$ from unit j to unit i reads

(3.1) \begin{equation} S_{ij}=U_{j} I^*_{ij} =U_j\left(\dfrac{U_j-U_i}{Z_{ij}}\right)^* =\dfrac{|U_j|e^{\imath\theta_j}\left(|U_j|e^{-\imath\theta_j}-|U_i|e^{-\imath\theta_i}\right)}{R_{ij}-\imath X_{ij}},\end{equation}

where $U_j=|U_j|e^{\imath\theta_j}$ and $I_{ij}=(U_j-U_i)/Z_{ij}$ denote the voltage at node j and the current between nodes j and i, respectively. Both are expressed as complex numbers to reflect the oscillating nature of AC power generation. The asterisks in, e.g., $I_{ij}^*$ indicate the complex conjugates (e.g., of $I_{ij}$ ). Moreover, $Z_{ij}=R_{ij}+\imath X_{ij}$ denotes the impedance of the transmission line (i, j) between unit i and j, with $R_{ij}$ denoting the resistance and $X_{ij}$ the reactance of (i, j). The complex power flow $S_{ij}=P_{ij}+\imath Q_{ij}$ consists of two parts: the active power $P_{ij}$ that results in net energy transfer and the reactive power $Q_{ij}$ that returns to the source in each cycle, not doing any work, but supporting the voltage stability of the power system [19].

The DC power flow model is based on the following assumptions on the parameters and the operating state of the power grid systems:

Assumption 3.1 (Perfect voltage support) The voltage amplitude is constant and identical for each node in the power grid network, $|U_i|\equiv U$ for all i, and the management of the reactive power is not considered.

Assumption 3.2 (Lossless lines) Transmission losses on the lines are neglected, implying that line resistances are negligible compared to line reactances: $R_{ij}/ X_{ij}\rightarrow 0$ for all lines (i, j).

Assumption 3.3 (Low line loads) Loads on all transmission lines are low, that is, the voltage angle differences between all neighbouring nodes are much smaller in magnitude than $\pi/2$ such that $\sin(\theta_j-\theta_i)\approx\theta_j-\theta_i$ and $\cos(\theta_j-\theta_i)\approx1$ .

With the above-mentioned assumptions in mind, the complex power flow (3.1) between neighbouring nodes simplifies to

(3.2) \begin{equation} S_{ij}=\dfrac{U^2(1-e^{\imath(\theta_j-\theta_i)})}{-\imath X_{ij}} =\dfrac{U^2}{X_{ij}}\left(\theta_j-\theta_i\right) = P_{ij}.\end{equation}

Here, the complex power flow $S_{ij}$ naturally reduces to the active power flow $P_{ij}$ since the imaginary part vanishes. The equation (3.2) resembles the expression of the direct current carried by a ‘resistor’ $X_{ij}/U^2$ influenced by a ‘voltage drop’ $\theta_j-\theta_i$ according to Ohm’s law, hence the name ‘DC power flow model’.

For an AC power grid system consisting of N units, the active power flow $P_i$ injected at unit i is the sum

(3.3) \begin{equation} P_i=\sum_{j=1}^{N} P_{ij}=\sum_{j=1}^N K_{ij}\left(\theta_i-\theta_j\right),\end{equation}

over all connected units. Equation (3.3) is the core of the DC power flow model as it yields the pattern of power flows $K_{ij}\theta_j$ across the grid network. Here, we follow the notation introduced in Section 2.1 and define the coupling strength as $K_{ij}=U^2/X_{ij}$ if there exists a transmission line between unit i and j and $K_{ij}=0$ if there is not. Denoting the nodal active power injections and nodal voltage angles as vectors, $\boldsymbol{P}:=(P_1,\cdots,P_N)$ and $\boldsymbol{\theta}:=(\theta_1,\cdots,\theta_N)$ , we express the linear relation between them by a weighted graph Laplacian $\mathcal{L}$ introduced in Section 2.2, i.e.

(3.4) \begin{equation} \boldsymbol{P}=\mathcal{L}\boldsymbol{\theta},\end{equation}

Here, $\mathcal{L}$ is defined similarly as in (2.9), only with $\gamma_{ij}\equiv 1$ for all edges (i, j).

Assuming that the power transmission network runs at a normal operation state where the voltage angles $\boldsymbol{\theta}^*$ are stationary at all nodes, the fixed voltage angle differences determine a specific power flow pattern $\boldsymbol{P}^*$ across the network through the linear operator $\mathcal{L}$ such that $\boldsymbol{P}^*=\mathcal{L}\boldsymbol{\theta}^*$ . If the nodal power injections are perturbed as $\boldsymbol{P}(t)=\boldsymbol{P}^*+\boldsymbol{D}(t)$ by a shift vector $\boldsymbol{D}$ (t) that in general is time dependent and reflects an increase or decrease of power generation for consumption at some of the nodes, the nodal voltage angles $\boldsymbol{\theta}=\boldsymbol{\theta}^*+\boldsymbol{\Theta}$ change accordingly through $\mathcal{L}$ due to the linear operation in (3.4). The response vector of the voltage angles $\boldsymbol{\Theta}$ is given by

(3.5) \begin{equation} \boldsymbol{\Theta}(t)=\mathcal{L}^{+}\boldsymbol{D}(t),\end{equation}

where $\mathcal{L}^{+}$ denotes the Moore–Penrose inverse of the weighted graph Laplacian $\mathcal{L}$ .

Remark 3.1 The weighted graph Laplacian $\mathcal{L}$ defined in (2.9) is singular since it always has an eigenvalue $\lambda_{0}=0$ with the corresponding eigenvector $\boldsymbol{v}_0=(1,\cdots,1)$ satisfying $\mathcal{L}\boldsymbol{v}_0=\boldsymbol{0}$ by construction. Therefore, to compute the voltage angle shifts (3.5), we need the generalised inverse matrix $\mathcal{L}^{+}$ , which can be computed by, e.g., the singular value decomposition. Alternatively, we can remove one dimension of the system by treating the phase $\theta_k$ of one of the units k as the reference for voltage angle, i.e. by setting $\theta_i\rightarrow\theta_i-\theta_k$ for all i. If the network is connected, then the $(N-1)$ -dimensional Laplacian matrix is invertible and both of the matrix equations, (3.4) and (3.5), have a unique solution, respectively.

3.2 LRT for the oscillator model of AC power grid dynamics

In this section, we discuss how LRT applies for the oscillator model of AC power grids, a widely used model for analysing the dynamics of AC power grids, and thereby provides a way to accurately determine the high-dimensional dynamic responses of an arbitrary power grid network to fluctuating power injections and withdrawals.

The dynamics of the high-voltage AC power transmission networks is essentially captured by an oscillator model (or second-order model) of AC power grids, of which synchronisation in terms of networked dynamical systems have been initially studied in references [10, 28] and [24]. This model allows for analytical understanding of the dynamics of power grids and has yielded fruitful research results over the past decade [28, 24, 7, 42, 37]. As the name suggests, in the oscillator model, each unit of AC power grids, a synchronous machine, is represented by an oscillator and the power transmission lines are represented by the pairwise couplings between the oscillators. The normal operation state of a power grid corresponds to the synchronisation of all oscillators, where all units rotate at the same frequency $\Omega_0^{m}$ corresponding to the nominal grid frequency $\Omega_0=2\pi\times50$ or $2\pi\times60$ Hz.

For each unit in the oscillator model, a synchronous machine, any change of the angular velocity of rotation results from the imbalance of the torques acting on the rotor operated at the nominal grid frequency. Its dynamics is governed by the so-called swing equation [16, 19]:

(3.6) \begin{equation} I\ddot{\theta}^{\text{m}}+D^{\text{m}}\dot{\theta}^{\text{m}}=T^{\text{m}}-T^{\text{e}},\end{equation}

where $\theta^{\text{m}}$ denotes the mechanical rotor angle deviation from the rotating reference frame $\Omega_0 t$ , I denotes the moment of inertia of the rotor and the connected turbine, $D^{\text{m}}$ denotes the coefficient of the damping torque resulting from the velocity-dependent friction at the air gap between the rotor and the stator in the synchronous machine. $T^{\text{m}}$ and $T^{\text{e}}$ denote, respectively, the net mechanical torque and the counteracting electromagnetic torque acting on the rotor.

For deviations of the angular velocity $\dot{\theta}^{\text{m}}$ , the local frequency deviation is small compared to the nominal grid frequency $\Omega_0^{\text{m}}$ , i.e. $\left(\dot{\theta}^{\text{m}}+\Omega_0^{\text{m}}\right)^{-1}\approx\left(\Omega_0^{\text{m}}\right)^{-1}$ , so that a torque T acting on the rotor can be expressed in terms of the power P of the synchronous machine as $T=\left(\dot{\theta}^{\text{m}}+\Omega_0^{\text{m}}\right)^{-1}P\approx\left(\Omega_0^{\text{m}}\right)^{-1}P$ . Also introducing the effective quantities $\theta(t)=\theta^{\text{m}}(t)/(p/2)$ and $\Omega_0=\Omega^{\text{m}}_0/(p/2)$ between the mechanical quantities and their electrical counterparts for a synchronous machine with p poles per phase, we obtain the more common version of the swing equation describing the dynamical relation between the rate of change of the electrical load angle and the power transmission between units:

(3.7) \begin{equation} M\ddot{\theta}+\tilde{D}\dot{\theta}=P^{\text{m}}-P^{\text{e}}.\end{equation}

Here, $M:=I\Omega_0/(p/2)^2$ and $\tilde{D}:=D^{\text{m}}\Omega_0/(p/2)^2$ are, respectively, the angular momentum of the rotor operated at the nominal grid frequency and the damping coefficient of the machine. On the right hand side of (3.7), $P^{\text{m}}$ denotes the net injected mechanical power (positive when the machine is operated as a generator and negative when operated as a motor), and $P^{\text{e}}$ denotes the electrical power injected to the grid by the synchronous machine. In the oscillator model, we again assume perfect voltage support (Assumption 3.1) and lossless transmission lines (Assumption 3.2), which yield $P^{\text{e}}=\sum_{j=1}^NK_{ij}(\theta_j-\theta_i)$ [cf. $P_i$ in the DC power flow model (3.3)]. Putting everything together, we obtain the governing equations of the oscillator model of AC power grids

(3.8) \begin{equation} \ddot{\theta}_i=P_i-\alpha_i\dot{\theta}_i+\sum_{j=1}^NK_{ij}\sin(\theta_j-\theta_i),\quad \text{for } i\in\{1,\cdots,N\}\end{equation}

with the parameters $P_i:=P_i^{\text{m}}/M_i$ , $\alpha_i:=\tilde{D}_i/M_i$ , and $K_{ij}:=U^2/X_{ij}$ .

Proposition 3.1 (Linear stability of the oscillator model) AC power grid systems described by the oscillator model (3.8) with underlying interaction topology G(V,E) is at least neutrally stable at a fixed point $\boldsymbol{\theta}^*$ , if all edges are not overloaded, i.e. $|\theta_j^*-\theta_i^*|\leqslant \frac{\pi}{2}$ for all $(i,\, j)\in E$ .

Proof. At a fixed point of the system $\boldsymbol{\theta}=\boldsymbol{\theta}^*$ , a small deviation of the oscillators’ angles $\boldsymbol{\Theta}:=\boldsymbol{\theta}-\boldsymbol{\theta}^*$ follows the linear dynamics

(3.9) \begin{align} \dfrac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix} \boldsymbol{\Theta} \\ \dot{\boldsymbol{\Theta}} \end{pmatrix} =\mathcal{J}\begin{pmatrix} \boldsymbol{\Theta} \\ \dot{\boldsymbol{\Theta}} \end{pmatrix}, \end{align}

where the Jacobian matrix $\mathcal{J}\in\mathbb{R}^{2N\times2N}$ of the 2N-dimensional dynamical system is given by

(3.10) \begin{align} \mathcal{J}= \begin{pmatrix} \boldsymbol{0}_N & \quad\boldsymbol{I}_N\\[3pt] -\mathcal{L} & \quad-\mathcal{A} \end{pmatrix}. \end{align}

Here, $\mathcal{L}$ is a weighted graph Laplacian as defined in (2.9) with $\gamma_{ij}=\cos\left(\theta^*_j-\theta^*_i\right)$ , $\mathcal{A}$ is an $N\times N$ diagonal matrix with $\mathcal{A}_{ii}:=\alpha_i$ , and $\boldsymbol{0}_N$ and $\boldsymbol{I}_N$ are, respectively, the $N\times N$ zero matrix and identity matrix.

Let $\boldsymbol{w}=(\boldsymbol{w}_1,\boldsymbol{w}_2)\in\mathbb{C}^{2N}$ with $\boldsymbol{w}_1,\boldsymbol{w}_2\in\mathbb{C}^{N}$ be an eigenvector of $\mathcal{J}$ corresponding to eigenvalue $\mu\in\mathbb{C}$ . By definition, we have $\mathcal{J}\boldsymbol{w}=\mu\boldsymbol{w}$ , which by writing (3.9) as a second-order differential equation implies

(3.11) \begin{equation} \mu^2\boldsymbol{w}_1+\mu\mathcal{A}\boldsymbol{w}_1+\mathcal{L}\boldsymbol{w}_1=\boldsymbol{0}. \end{equation}

Multiplying both sides of the equation above with the conjugate transpose, $\boldsymbol{w}_1^{\dagger}$ of $\boldsymbol{w}_1$ from the left, we obtain an expression of the eigenvalue $\mu=\left(-\chi_2\pm\sqrt{\chi_2^2-4\chi_1\chi_3}\right)/2\chi_1$ , with $\chi_1=\boldsymbol{w}^{\dagger}_1\boldsymbol{w}_{\boldsymbol{1}}\geqslant 0$ , $\chi_2=\boldsymbol{w}^{\dagger}_1\mathcal{A}\boldsymbol{w}_1\geqslant 0$ and $\chi_3=\boldsymbol{w}^{\dagger}_1\mathcal{L}\boldsymbol{w}_1\geqslant 0$ . $\chi_2$ and $\chi_3$ are non-negative since $\mathcal{A}_{ii}=\alpha_i>0$ and $\mathcal{L}$ is positive semi-definite because $\gamma_{ij}\geqslant 0$ is ensured through $|\theta_j^*-\theta_i^*|\leqslant \frac{\pi}{2}$ for all $(i, j)\in E$ (see Proposition 2.1 and Remark 2.2) [cf. [22]]. Therefore, the eigenvalue always has a non-positive real part, implying that the networked dynamical system is at least neutrally stable at the fixed point.

Remark 3.2 (Neutral stability and global phase shift) Connected AC power grids described by the oscillator model (3.8) is neutrally stable at a fixed point only when the deviation $\boldsymbol{\Theta}$ is a global phase shift. Because (3.11) indicates that the Jacobian eigenvalue $\mu=0$ only when $\chi_3=0$ , which implies $\mathcal{L}\boldsymbol{w}_1=\boldsymbol{0}$ and thus $\boldsymbol{w}_1$ lies in the Laplacian eigenspace corresponding to the eigenvalue $\lambda^{[0]}=0$ . For connected graphs, there is only one eigenvector $\boldsymbol{v}^{[0]}=(1,\cdots,1)$ corresponding to $\lambda^{[0]}=0$ , therefore the system is only neutrally stable when all nodes undergo a phase shift with the same magnitude. A global phase shift has no effects on the power flow pattern in the network since pairwise phase differences across edges remain the same.

Proposition 3.2 (LRT of the oscillator model and homogeneous nodal damping) Consider an AC power grid oscillator model with arbitrary topology (3.8) with homogeneous nodal damping $\alpha_i=\alpha\geqslant 0$ for all nodes i. Then the network-wide linear responses to arbitrary external perturbations near a normal operation state $\boldsymbol{\theta}^*$ can be expressed explicitly in the eigenbasis of a weighted graph Laplacian: i) The network response to time-independent distributed perturbations $\boldsymbol{D}^*$ is

(3.12) \begin{equation} \boldsymbol{\Theta}(t)=\boldsymbol{D}^*\cdot\boldsymbol{v}^{[0]}\left(\dfrac{e^{-\alpha t}}{\alpha^2}-\dfrac{1}{\alpha^2}+\dfrac{t}{\alpha}\right)\boldsymbol{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\boldsymbol{D}^*\cdot\boldsymbol{v}^{[\ell]}}{\lambda^{[\ell]}}\left(\dfrac{\Delta_-^{[\ell]}e^{\Delta_+^{[\ell]}t}-\Delta_+^{[\ell]}e^{\Delta_-^{[\ell]}t}}{2\eta^{[\ell]}}+1\right)\boldsymbol{v}^{[\ell]}, \end{equation}

and ii) the network response to a single sinusoidal perturbation given by $\boldsymbol{D}(t)$ with $D_{i}(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$ is

(3.13) \begin{equation} \boldsymbol{\Theta}^{(k)}(t)=\sum_{\ell=0}^{N-1}\dfrac{\varepsilon v_k^{[\ell]}e^{\imath\varphi}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}\left[ \dfrac{\left(\Delta^{[\ell]}_--\imath\omega\right)e^{\Delta^{[\ell]}_+t}-\left(\Delta^{[\ell]}_+-\imath\omega\right)e^{\Delta^{[\ell]}_-t}}{2\eta^{[\ell]}}+e^{\imath\omega t}\right]\boldsymbol{v}^{[\ell]} \end{equation}

with $\Delta^{[\ell]}_{\pm}:=-\alpha/2\pm\eta^{[\ell]}$ and $\eta^{[\ell]}:=\sqrt{\alpha^2/4-\lambda^{[\ell]}}$ .

Proof. Vectorising the linear response $\boldsymbol{\Theta}(t)$ of the system (3.8) to a perturbation vector $\boldsymbol{D}(t)$ , we obtain the matrix equation describing the response dynamics of the oscillator model of AC power grids

(3.14) \begin{equation} \ddot{\boldsymbol{\Theta}}+\boldsymbol{\alpha}\circ\dot{\boldsymbol{\Theta}}=-\mathcal{L}\boldsymbol{\Theta}+\boldsymbol{D}(t) , \end{equation}

where $\boldsymbol{\alpha}:=(\alpha_1,\cdots,\alpha_N)$ denotes the vector of damping parameters and ‘ $\circ$ ’ denotes the Schur (element-wise) product of two vectors. Let $\alpha_i=\alpha$ for all i, the term $\boldsymbol{\alpha}\circ\dot{\boldsymbol{\Theta}}$ reduces to a scalar multiplication $\alpha\dot{\boldsymbol{\Theta}}$ , thereby all terms involving the variable $\boldsymbol{\Theta}$ in equation (3.14) can be expressed as linear combinations of Laplacian eigenvectors. Because the $\mathcal{L}$ here is real and symmetric so that we can write $\boldsymbol{\Theta}(t)=\sum_{\ell=0}^{N-1}c^{[\ell]}(t)\boldsymbol{v}^{[\ell]}$ . Using the same trick as in Remark 2.3, we obtain equations for the time-dependent coefficients $c^{[\ell]}(t)$ given by

(3.15) \begin{equation} \ddot{c}^{[\ell]}+\alpha\dot{c}^{[\ell]}+\lambda^{[\ell]} c^{[\ell]}=\boldsymbol{D}(t)\cdot\boldsymbol{v}^{[\ell]}\quad \text{for }\ell\in\{0,\cdots,N-1\}. \end{equation}

Assuming at $t=0$ the AC power grid system operates normally at the fixed point $\boldsymbol{\theta}^*$ , we have initial conditions for the coefficients $c^{[\ell]}(0)=0$ and $\dot{c}^{[\ell]}(0)=0$ for all $\ell$ , thereby obtain explicitly solutions for the coefficients and for the network linear responses.

For perturbations independent of time, $\boldsymbol{D}(t)=\boldsymbol{D}^*$ , such as constant shifts in power injection and consumption, the linear response of a power grid system is given by directly solving for the coefficients in (3.15):

(3.16) \begin{equation} \boldsymbol{\Theta}(t)=\boldsymbol{D}^*\cdot\boldsymbol{v}^{[0]}\left(\dfrac{e^{-\alpha t}}{\alpha^2}-\dfrac{1}{\alpha^2}+\dfrac{t}{\alpha}\right)\boldsymbol{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\boldsymbol{D}^*\cdot\boldsymbol{v}^{[\ell]}}{\lambda^{[\ell]}}\left(\dfrac{\Delta_-^{[\ell]}e^{\Delta_+^{[\ell]}t}-\Delta_+^{[\ell]}e^{\Delta_-^{[\ell]}t}}{2\eta^{[\ell]}}+1\right)\boldsymbol{v}^{[\ell]}, \end{equation}

with $\Delta^{[\ell]}_{\pm}:=-\alpha/2\pm\eta^{[\ell]}$ and $\eta^{[\ell]}:=\sqrt{\alpha^2/4-\lambda^{[\ell]}}$ . For time-dependent perturbations $\boldsymbol{D}(t)$ , such as fluctuating power injections from renewables, we obtain the network linear response based on the responses to each single frequency components at each perturbed nodes, as discussed in Remark 2.4. Similarly, we let $D_{i}(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$ and obtain the oscillator model’s linear response to a sinusoidal signal at node k as

(3.17) \begin{equation} \boldsymbol{\Theta}^{(k)}(t)=\sum_{\ell=0}^{N-1}\dfrac{\varepsilon v_k^{[\ell]}e^{\imath\varphi}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}\left[ \dfrac{\left(\Delta^{[\ell]}_--\imath\omega\right)e^{\Delta^{[\ell]}_+t}-\left(\Delta^{[\ell]}_+-\imath\omega\right)e^{\Delta^{[\ell]}_-t}}{2\eta^{[\ell]}}+e^{\imath\omega t}\right]\boldsymbol{v}^{[\ell]}. \end{equation}

Remark 3.3 (Low dissipation regime and grid eigenfrequencies) In case the dissipation in the system is low enough such that $\alpha<2\sqrt{\lambda^{[\ell]}}$ for the $\ell$ -th eigenvalue, in the solution of the linear response (3.13) the corresponding $\eta^{[\ell]}$ for the same eigenmode becomes imaginary, suggesting this mode is oscillating under-damped in the power grid system with an exponentially decaying amplitude proportional to $e^{-\frac{\alpha}{2}t}$ , i.e., with a time constant $\tau=2/\alpha$ . Such intrinsic oscillation modes can also be identified by looking at the eigenvalues of the Jacobian matrix at the fixed point. Let $\boldsymbol{w}_1$ in (3.11) be the $\ell$ -th eigenvector of the Laplacian $\mathcal{L}$ , we can see the corresponding Jacobian eigenvalue $\mu^{[\ell]}=-\alpha/2\pm\sqrt{\alpha^2/4-\lambda^{[\ell]}}=\Delta^{[\ell]}_{\pm}$ , which indicates the corresponding eigenfrequency $\omega^{[\ell]}_{\text{eigen}}:=\hbox{Im}[\mu^{[\ell]}]=\sqrt{\lambda^{[\ell]}-\alpha^2/4}$ of the power grid system. Since for a connected networked system with N nodes has $N-1$ positive Laplacian eigenvalues, it also has a band of $N-1$ eigenfrequencies if the dissipation is sufficiently low satisfying $\alpha<2\sqrt{\lambda^{[\ell]}}$ for all $N-1$ Laplacian eigenvalues.

4 Emerging network response patterns from LRT

In Section 3, we applied LRT on power grid models and obtained explicit solutions for linear network responses to perturbations at a normal operation state. The solutions are expressed in terms of the eigensystem of a weighted graph Laplacian. These Laplacians and thus their eigensystems contain information about the underlying network topology as well as the base operating state of the system, enabling us to understand and manipulate how complex networked systems such as power grids collectively respond to external perturbation signals.

In this section, we focus on the dynamic responses of AC power grids to time-varying perturbations based on the oscillator model (see Section 3.2) and explicate how steady-state response patterns constituted by the set of nodal response magnitudes as well as transient response patterns describing the spatio-temporal spreading of a perturbation in power grids are mathematically extracted from the solution given by the LRT.

4.1 Frequency regimes of steady-state response patterns

After a transient phase characterised by a dissipation-related time constant $\tau=2/\alpha$ (cf. Remark 3.3), the perturbed power grid systems reside in a second regime of network responses, where the entire network respond periodically to perturbation signals for $t\gg\tau$ (see 3.12 and 3.13). We thus call the network responses for such large times steady-state responses. We remark that the steady-state responses here are not necessarily stationary, meaning that the nodal responses themselves can vary with time, but their characteristics, such as the amplitude and the phase of sinusoidal responses, constitute network-wide response patterns that are time-independent.

Proposition 4.1 (Steady-state response pattern for a constant perturbation) For AC power grids with arbitrary topologies (3.8), the steady-state responses to time-independent perturbations $\boldsymbol{D}(t)=\boldsymbol{D}^*$ , i.e. with a perturbation frequency $\omega=0$ , near a normal operation state $\boldsymbol{\theta}^*$ are constituted by a homogeneous shift of grid frequency

(4.1) \begin{equation} \delta\dot{\theta}_i=\dfrac{1}{N\alpha}\displaystyle\sum_{j=1}^N D^*_j\quad\text{for }i\in\{1,\cdots,N\}, \end{equation}

and a topology-dependent phase shift

(4.2) \begin{equation} \delta\boldsymbol{\theta}=-\dfrac{1}{N\alpha^2}\displaystyle\sum_{i=1}^N D^*_i\boldsymbol{1}+\sum_{\ell=1}^{N-1}\dfrac{\boldsymbol{D}^*\cdot\boldsymbol{v}^{[\ell]}}{\lambda^{[\ell]}}\boldsymbol{v}^{[\ell]}. \end{equation}

Proof. By definition, the steady-state response patterns become clear by investigating the asymptotic behaviour of the responses as $t\rightarrow\infty$ . For responses to a constant perturbation vector $\boldsymbol{D}(t)=\boldsymbol{D}^*$ (3.12), the steady-state response reads

(4.3) \begin{equation} \boldsymbol{\Theta}(t)\overset{t\rightarrow\infty}{\sim} \boldsymbol{D}^*\cdot\boldsymbol{v}^{[0]}\left(-\dfrac{1}{\alpha^2}+\dfrac{t}{\alpha}\right)\boldsymbol{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\boldsymbol{D}^*\cdot\boldsymbol{v}^{[\ell]}}{\lambda^{[\ell]}}\boldsymbol{v}^{[\ell]},\end{equation}

which consists of two characteristic patterns: i) the phases of all units drift away from the normal operation state with a constant angular velocity $\left(\boldsymbol{D}^*\cdot\boldsymbol{v}^{[0]}\right)v^{[0]}_i/\alpha$ , and ii) a time-independent and unit-specific phase shift. Since $\boldsymbol{v}^{[0]}=\frac{1}{\sqrt{N}}\boldsymbol{1}$ , the former pattern represents the constant homogeneous shift of grid frequency (4.1) and the latter the topology-dependent phase shift (4.2).

Remark 4.1 The global grid frequency shift (4.1) is a consequence of the imbalance between power injected into and drawn from the power grid system, which is imposed by the constant perturbation $\boldsymbol{D}^*$ . The topology-dependent phase shifts (4.2) at all units suggest a network-wide redistribution of power flows on the transmission lines. The first-order approximation of the change of the load $L_{ij}:=\sin(\theta_j-\theta_i)$ on line (i, j), $\delta L_{ij}=\cos\big(\theta_j^*-\theta_i^*\big)\left(\delta\theta_j-\delta\theta_i\right)$ , provides an indicator for the emerging risks such as overheating for heavily-loaded lines with $L_{ij}$ approaching the upper limit of its safety range.

Proposition 4.2 (Steady-state response pattern for a sinusoidal perturbation) For AC power grids with arbitrary topologies (3.8), the steady-state responses to a sinusoidal perturbation at node k, i.e. $\boldsymbol{D}(t)$ with $D_i(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$ , near a normal operation state $\boldsymbol{\theta}^*$ are constituted by a homogeneous phase shift

(4.4) \begin{equation} \delta\theta_i=\dfrac{\imath\varepsilon e^{\imath\varphi}}{\alpha\omega N}\quad\text{for }i\in\{1,\cdots,N\}, \end{equation}

and sinusoidal responses at each node with the same frequency $\omega$ and a characteristic complex amplitude

(4.5) \begin{equation} R_i^{(k)}(\omega):=\sum_{\ell=0}^{N-1}\dfrac{ v_k^{[\ell]}v_i^{[\ell]}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}. \end{equation}

for node i.

Proof. The steady-state response to a sinusoidal perturbation at a single node k is obtained by studying the asymptotic behaviour of the responses (3.13) has the form of

(4.6) \begin{equation} \boldsymbol{\Theta}^{(k)}(t)\overset{t\rightarrow\infty}{\sim}\dfrac{\imath\varepsilon e^{\imath\varphi}}{\alpha\omega N}\boldsymbol{1}+e^{\imath(\omega t+\varphi)}\sum_{\ell=0}^{N-1}\dfrac{\varepsilon v_k^{[\ell]}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}\boldsymbol{v}^{[\ell]},\end{equation}

which is composed of a homogeneous phase shift $\frac{\imath\varepsilon e^{\imath\varphi}}{\alpha\omega N}$ and a driven oscillation at each node. Each node’s angular variable $\theta_i$ changes at the same frequency as the perturbation frequency $\omega$ , but with a complex amplitude

(4.7) \begin{equation} R_i^{(k)}(\omega)=\sum_{\ell=0}^{N-1}\dfrac{ v_k^{[\ell]}v_i^{[\ell]}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}.\end{equation}

Remark 4.2 (Characterisation of the steady-state response patterns to a sinusoidal perturbation) The complex nature of the amplitude suggests shifts in the amplitude and in the phase between the perturbation signal and the response signal, which are both topology-dependent and node-specific. Hence, we also refer to $R_i^{(k)}$ as the nodal response factor of node i to a sinusoidal perturbation at node k.

The homogeneous phase shift contributes to neither the change of grid frequency nor the overall power flow pattern in the network, while the absolute value of $R_i^{(k)}$ determines the maximal deviation of the local grid frequency at node i caused by a perturbation at node k through $|\delta\dot{\theta}_i|=|\dot{\Theta}^{(k)}_i|=\varepsilon\omega|R_i^{(k)}|$ . If it exceeds the safety range of normal operation, the local frequency deviation may cause damage the synchronous machine and other related grid components such as the turbine. In the rest of the subsection, we focus on the steady-state response pattern constituted by the set of nodal response strengths

(4.8) \begin{equation} A_i^{(k)}:=\omega\left|R_i^{(k)}\right|, \end{equation}

for each node i and discuss in detail its distinctive spatial distributions in different frequency regimes (see Fig. 1 for an example).

Figure 1. Steady-state response patterns exhibit three frequency regimes. (a) Relative response strength $A_i^{*(k)}$ (see Remark 4.4 for definition) of all 80 nodes in an example power grid network across all three frequency regimes (homogeneous bulk, resonant and localised responses). Vertical grey lines represent the $N-1$ eigenfrequencies. (b1,c1,d1) Qualitatively different dependencies of $A_i^{*(k)}$ on the graph-theoretic distance $d:=d(k,i)$ with three representative driving frequencies $\omega/2\pi\in\{0.1,2,10\}$ Hz of three frequency regimes. The exponential dependence of $A_i^{*(k)}$ on d is illuatrated in the inset of d1. (b2,c2,d2) Distinctive response patterns for the three driving frequencies, corresponding to (b1,c1,d1). The curves in (a) are colour coded with the distance d, and the discs in (b1-b2,c1-c2,d1-d2) are colour coded with the relative response strength $A_i^{*(k)}$ . The black square marks the perturbed node. Network settings are the same as Figure 2 in [45].

4.1.1 Regime of grid resonances

As suggested in Remark 3.3, the dynamics of the perturbed oscillator model of AC power grids can be understood in comparison with the dynamics of a driven damped harmonic oscillator. For each intrinsic under-damped oscillation mode corresponding to a non-zero Laplacian eigenvalue $\lambda^{[\ell]}>\alpha^2/4$ , the oscillatory power grid system resonates if the perturbation frequency matches the corresponding eigenfrequency $\omega_{\text{eigen}}^{[\ell]}=\sqrt{\lambda^{[\ell]}-\alpha^2/4}$ . Driven at frequencies close to an eigenfrequency $\omega_{\text{eigen}}^{[\ell]}$ , the network responses exhibit large amplitudes since the rationalised denominator $\left(\omega^2-\lambda^{[\ell]}\right)^2-\alpha^2\omega^2$ for the corresponding $\ell$ -th oscillation mode is minimised. However, the response amplitudes vary greatly for different nodes in the network due to the factor $v_k^{[\ell]}v_i^{[\ell]}$ . If the system dissipation is sufficiently low such that there exist $N-1$ under-damped oscillation modes, the corresponding $N-1$ eigenfrequencies form a resonance regime where the power grid system can potentially exhibit strong distributed responses across the network (cf. Remark 3.3).

We emphasise that the spatiotemporal resonance pattern is characteristic of each perturbation frequency within the resonance regime, of each specific network topology including the prior perturbation base state and of each location of perturbation. Therefore, a power grid system’s responses to real-world fluctuating perturbations containing a collection of frequency components within the resonance regime are quite irregular both temporally and spatially, which makes it a non-trivial task to evaluate the risks in perturbed power grids induced by network resonances (see Section 4.3 for further discussions).

4.1.2 Homogeneous responses: The low-frequency regime

The network response pattern for lower perturbation frequencies, i.e. the ones lower than the smallest eigenfrequency $\omega_{\text{eigen}}^{[1]}=\sqrt{\lambda^{[1]}-\alpha^2/4}$ , can be understood by investigating the asymptotic behaviour of the nodal response strength $A_i^{(k)}=\omega|R_i^{(k)}|$ as the $\omega\rightarrow 0$ .

Proposition 4.3 (Homogeneous responses at the low frequency limit) As the perturbation frequency $\omega\rightarrow0$ , the steady-state response strength at each node of an AC power grid system with an arbitrary topology (3.8) approaches a constant value, i.e.

(4.9) \begin{equation} A_i^{(k)}\overset{\omega\rightarrow 0}{\sim}\frac{1}{N\alpha}. \end{equation}

Proof. Considering (4.5) and $\lambda^{[0]}=0$ , we derive the asymptotic behaviour of the real part and the imaginary part of the nodal response factor $R_i^{(k)}$ as $\omega\rightarrow 0$ :

(4.10) \begin{align} \hbox{Re}\left[R_i^{(k)}\right]&=\sum_{\ell=0}^{N-1}\frac{v^{[\ell]}_kv^{[\ell]}_i\left(-\omega^2+\lambda^{[\ell]}\right)}{\left(-\omega^2+\lambda^{[\ell]}\right)^2+\alpha^2\omega^2}\overset{\omega\rightarrow 0}{\sim}-\dfrac{1}{N\alpha^2}+\displaystyle\sum_{\ell=1}^{N-1}\dfrac{v^{[\ell]}_kv^{[\ell]}_i}{\lambda^{[\ell]}},\end{align}

(4.11) \begin{align} \,\,\,\hbox{Im}\left[R_i^{(k)}\right]&=\sum_{\ell=0}^{N-1}\frac{v^{[\ell]}_kv^{[\ell]}_i\left(-\alpha\omega\right)}{\left(-\omega^2+\lambda^{[\ell]}\right)^2+\alpha^2\omega^2}\overset{\omega\rightarrow 0}{\sim}-\dfrac{1}{N\alpha\omega}.\qquad\qquad\qquad \end{align}

As the asymptotic behaviour of the response strength $A_i^{(k)}=\omega|R_i^{(k)}|$ is dominated by the imaginary part, we have

(4.12) \begin{equation} A_i^{(k)}=\omega\left|R_i^{(k)}\right|\overset{\omega\rightarrow 0}{\sim}\omega\cdot\dfrac{1}{N\alpha\omega}=\dfrac{1}{N\alpha}. \end{equation}

Remark 4.3 (Consistency with the homogeneous grid frequency shift at $\omega=0$ ) The homogeneous response strength at each node as $\omega\rightarrow 0$ suggests a global shift of grid frequency inversely proportional to the network size N and the system dissipation parameter $\alpha$ . This result is quantitatively consistent with the homogeneous grid frequency shift induced by constant perturbations, as discussed in Proposition 4.1 with $D_i^*=\delta_{ik}$ if node k is the perturbed one.

Remark 4.4 (Relative nodal response strength) The homogeneous nodal response at the low frequency limit (4.9) serves as a reference value for the nodal response strengths of a sinusoidally driven network. Therefore, by normalising $A_i^{(k)}$ with its low frequency limit $\frac{1}{N\alpha}$ , we define the relative nodal response strength

(4.13) \begin{align} A_i^{*(k)}:=\dfrac{A_i^{(k)}}{\displaystyle\lim_{\omega\rightarrow0}A_i^{(k)}}=N\alpha A_i^{(k)}, \end{align}

so that the nodal response strengths can be compared across networks with different sizes and dissipation values.

4.1.3 Localised responses: The high-frequency regime

Similarly, we investigate the network response pattern in the high-frequency regime where $\omega>\omega_{\text{eigen}^{[N-1]}}$ by observing the asymptotic behaviour of the nodal response strengths $A_i^{(k)}$ as the perturbation frequency becomes sufficiently large and approaches infinity.

Proposition 4.4 (Localised response patterns in the high-frequency regime) As the perturbation frequency $\omega\rightarrow\infty$ , the steady-state nodal response strength $A_i^{(k)}$ in an AC power grid system with an arbitrary topology (3.8) decays as a power-law of $\omega$ with an exponent depending on the graph-theoretic distance d between node k and i, i.e.

(4.14) \begin{equation} A_i^{(k)}\overset{\omega\rightarrow\infty}{\sim}\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1}, \end{equation}

where $\left|\Phi^{[d]}_{ki}\right|$ is a distance- and node-specific prefactor but independent on the perturbation frequency.

Proof. To determine the asymptotic behaviour of the response strength $A_i^{(k)}=\omega|R_i^{(k)}|$ , we first reduce $\hbox{Re}[R_i^{(k)}]$ and $\hbox{Im}[R_i^{(k)}]$ from (4.10) and (4.11) to a common denominator $M(\omega)$ and obtain the respective numerators $N_{\text{Re}}(\omega)$ and $N_{\text{Im}}(\omega)$ as polynomials of $\omega$ ,

(4.15) \begin{align} \ \ M(\omega)&:=\displaystyle\prod_{\ell=0}^{N-1}\left[\left(-\omega^2+\lambda^{[\ell]}\right)^2+\alpha^2\omega^2\right],\quad\qquad\end{align}

(4.16) \begin{align} N^{\text{Re}}_{ki}(\omega)&:=\displaystyle\sum_{\ell=0}^{N-1}v_k^{[\ell]}v_i^{[\ell]}\left(-\omega^2+\lambda^{[\ell]}\right)Q^{[\ell]}(\omega), \text{ and}\end{align}

(4.17) \begin{align} \ \ N^{\text{Im}}_{ki}(\omega)&:=\displaystyle\sum_{\ell=0}^{N-1}v_k^{[\ell]}v_i^{[\ell]}\left(-\alpha\omega\right)Q^{[\ell]}(\omega) \quad\text{with}\qquad\end{align}

(4.18) \begin{align} Q^{[\ell]}(\omega)&:=\prod_{\substack{\ell'=0,\ell'\neq \ell}}^{N-1}\left[\left(-\omega^2+\lambda^{\big[\ell'\big]}\right)^2+\alpha^2\omega^2\right]. \end{align}

The asymptotic behaviour of $M(\omega)$ , $N^{\text{Re}}_{ki}(\omega)$ and $N^{\text{Im}}_{ki}(\omega)$ as $\omega\rightarrow\infty$ is dominated by the respective leading terms with the highest power of $\omega$ . The denominator scales asymptotically as

(4.19) \begin{equation} M(\omega)\overset{\omega\rightarrow\infty}{\sim}\omega^{4N}.\end{equation}

For the numerators, the leading term depends on a common product $Q^{[\ell]}(\omega)$ . As shown in Appendix A, $Q^{[\ell]}(\omega)$ can be written in terms of $\lambda^{[\ell]}$ and other variables that are dependent on the underlying matrix $\mathcal{L}$ but independent of the choice of $\ell$ . Thus, we define

(4.20) \begin{equation} Q(\lambda^{[\ell]},\omega):=\sum_{j=0}^{2N-2}C^{[j]}\left(\lambda^{[\ell]}\right)\omega^{4N-4-2j},\end{equation}

where the coefficients $C^{[j]}\left(\lambda^{[\ell]}\right)$ are polynomials in $\lambda^{[\ell]}$ of degree j. Inserting the expression of $Q(\lambda^{[\ell]},\omega)$ (4.20) to the numerators (4.16) and (4.17), we obtain

(4.21) \begin{align} N^{\text{Re}}_{ki}&=\sum_{\ell=0}^{N-1}v_k^{[\ell]}v_i^{[\ell]}\sum_{j=0}^{2N-1}F^{[j]}\left(\lambda^{[\ell]}\right)\omega^{4N-2-2j}\text{ and } N^{\text{Im}}_{ki}&=\sum_{\ell=0}^{N-1}v_k^{[\ell]}v_i^{[\ell]}\sum_{j=0}^{2N-2}G^{[j]}\left(\lambda^{[\ell]}\right)\omega^{4N-3-2j},\end{align}

where $F^{[j]}\left(\lambda^{[\ell]}\right)$ and $G^{[j]}\left(\lambda^{[\ell]}\right)$ are also polynomials in $\lambda^{[\ell]}$ of degree j and can be written in terms of $C^{[j]}\left(\lambda^{[\ell]}\right)$ as

(4.22) \begin{align} F^{[j]}\left(\lambda^{[\ell]}\right)&= \begin{cases} -C^{[j]}\left(\lambda^{[\ell]}\right) & j=0\\[3pt] -C^{[j]}\left(\lambda^{[\ell]}\right)+\lambda^{[\ell]}C^{j-1}\left(\lambda^{[\ell]}\right) & 1\leqslant j \leqslant 2N-2,\\[3pt] \lambda^{[\ell]}C^{j-1}\left(\lambda^{[\ell]}\right) & j=2N-1 \end{cases} \end{align}

(4.23) \begin{align} \,\,G^{[j]}\left(\lambda^{[\ell]}\right)&=-\alpha C^{[j]}\left(\lambda^{[\ell]}\right),\quad j\in\{0,1,\cdots,2N-2\}. \qquad\qquad \end{align}

Considering the numerators $N^{\text{Re}}_{ki}(\omega)$ and $N^{\text{Im}}_{ki}(\omega)$ as the ki-th elements of numerator matrices $\mathcal{N}^{\text{Re}}$ and $\mathcal{N}^{\text{Im}}$ , we can conveniently write (4.21) in a matrix form

(4.24) \begin{equation} \mathcal{N}^{\text{Re}}=\sum_{j=0}^{2N-1}{\Phi}^{[j]}\omega^{4N-2-2j},\quad \mathcal{N}^{\text{Im}}=\sum_{j=0}^{2N-2}{\Gamma}^{[j]}\omega^{4N-3-2j}, \end{equation}

with the coefficient matrices ${\Phi}^{[j]}:=\mathcal{V}\mathcal{F}^{[j]}\mathcal{V}^{\text{T}}$ and ${\Gamma}^{[j]}:=\mathcal{V}\mathcal{G}^{[j]}\mathcal{V}^{\text{T}}$ . Here $\mathcal{V}:=\left(v^{[0]},\cdots,v^{[N-1]}\right)$ and $\mathcal{F}^{[j]}$ , $\mathcal{G}^{[j]}$ are diagonal matrices with $\mathcal{F}^{[j]}_{ii}:=F^{[j]}\left(\lambda^{[i]}\right)$ and $\mathcal{G}^{[j]}_{ii}:=G^{[j]}\left(\lambda^{[i]}\right)$ , respectively. By spelling out the polynomials

(4.25) \begin{equation} F^{[j]}\left(\lambda^{[\ell]}\right)=\sum_{m=0}^jf^{[j]}_m\cdot\left(\lambda^{[\ell]}\right)^m,\quad G^{[j]}\left(\lambda^{[i]}\right)=\sum_{m=0}^jg^{[j]}_m\cdot\left(\lambda^{[\ell]}\right)^m, \end{equation}

with coefficients $f_m^{[j]},g_m^{[j]}\in\mathbb{R}$ , we can see the diagonal matrices $\mathcal{F}^{[j]}$ and $\mathcal{G}^{[j]}$ are polynomials of a diagonal matrix $\Lambda$ with $\Lambda_{ii}:=\lambda^{[i]}$

(4.26) \begin{equation} \mathcal{F}^{[j]}=\sum_{m=0}^jf^{[j]}_m\Lambda^m,\quad \mathcal{G}^{[j]}=\sum_{m=0}^jg^{[j]}_m\Lambda^m, \end{equation}

so that the coefficient matrices $\Phi^{[j]}$ and $\Gamma^{[j]}$ can be written as

(4.27) \begin{equation} {\Phi}^{[j]} = \mathcal{V}\left(\sum_{m=0}^{j}f^{[j]}_m\Lambda^m\right)\mathcal{V}^{\text{T}}=\sum_{m=0}^{j}f^{[j]}_m\left(\mathcal{V}\Lambda^m\mathcal{V}^{\text{T}}\right)=\sum_{m=0}^{j}f^{[j]}_m\mathcal{L}^m,\end{equation}

(4.28) \begin{equation} {\Gamma}^{[j]} = \mathcal{V}\left(\sum_{m=0}^{j}g^{[j]}_m\Lambda^m\right)\mathcal{V}^{\text{T}}=\sum_{m=0}^{j}g^{[j]}_m\left(\mathcal{V}\Lambda^m\mathcal{V}^{\text{T}}\right)=\sum_{m=0}^{j}g^{[j]}_m\mathcal{L}^m, \end{equation}

indicating that they are polynomials of the weighted graph Laplacian matrix $\mathcal{L}$ of degree j, as $\mathcal{V}\Lambda\mathcal{V}^{\text{T}}\equiv\mathcal{L}$ . In short, the numerators $\mathcal{N}^{\text{Re}}$ and $\mathcal{N}^{\text{Im}}_{ki}$ are in fact polynomials in $\omega$ , with its coefficient being also a polynomial in $\mathcal{L}$ ,

(4.29) \begin{equation} \mathcal{N}^{\text{Re}}=\sum_{j=0}^{2N-1}\left(\sum_{m=0}^{j}f^{[j]}_m\mathcal{L}^m\right)\omega^{4N-2-2j},\quad \mathcal{N}^{\text{Im}}=\sum_{j=0}^{2N-2}\left(\sum_{m=0}^{j}g^{[j]}_m\mathcal{L}^m\right)\omega^{4N-3-2j}. \end{equation}

We note that as the index j increases, the powers of $\omega$ decrease and the degrees of the polynomials $\Phi^{[j]}(\mathcal{L})$ and $\Gamma^{[j]}(\mathcal{L})$ , i.e. the coefficients of $\omega$ , increase. For sufficiently large perturbation frequency $\omega$ , the leading terms in the numerators which dominate the asymptotic behaviours would be the ones with the highest degrees of $\omega$ with nonzero coefficients. We know from graph theory that $(\mathcal{L}^m)_{ij}\neq0$ only for node pair (i, j) with the graph theoretic distance d(i, j) between them, i.e. the length of the shortest path from j to i on the unweighted graph defined by $\mathcal{L}$ , satisfying $d(i, j)\leqslant m$ . Therefore, the first terms in $N^{\text{Re}}_{ki}(\omega)$ and $N^{\text{Im}}_{ki}(\omega)$ with the highest degrees of $\omega$ have exactly zero coefficients, i.e. $\Phi^{[j]}=0$ and $\Gamma^{[j]}=0$ , for all $j<d(k,i)$ . Consequently, the leading term in the numerators are

(4.30) \begin{equation} N^{\text{Re}}_{ki}(\omega)\overset{\omega\rightarrow\infty}{\sim}\Phi^{[d]}_{ki}\omega^{4N-2-2d},\quad N^{\text{Im}}_{ki}(\omega)\overset{\omega\rightarrow\infty}{\sim}\Gamma^{[d]}_{ki}\omega^{4N-3-2d}, \end{equation}

with $d:=d(k,i)$ . Together with the asymptotic behaviour of the denominator (4.19), we obtain

(4.31) \begin{equation} A_i^{(k)}=\omega\left|R_i^{(k)}\right|\overset{\omega\rightarrow \infty}{\sim}\omega\cdot\left|\dfrac{\Phi^{[d]}_{ki}\omega^{4N-2-2d}}{\omega^{4N}}\right|=\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1}. \end{equation}

Remark 4.5 (Localised response patterns in the high-frequency limit) Proposition 4.4 implies that the response amplitude of the grid frequency $A_i^{(k)}$ decays exponentially with the graph-theoretic distance d between the perturbed node and the responding node. The response amplitude also decays as a power law in $\omega$ for fixed d. Visualizations of such localised response patterns in three networks are given in Fig. 2. In the high-frequency limit, a network’s response to the perturbation is restricted to the perturbed node, i.e.

(4.32) \begin{align} \lim_{\omega\rightarrow\infty}\dfrac{A_i^{(k)}}{A_i^{(k)}} \overset{\omega\rightarrow \infty}{\sim}\lim_{\omega\rightarrow\infty}\dfrac{\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1}}{\left|\Phi^{[0]}_{ki}\right|\omega^{-1}} =\lim_{\omega\rightarrow\infty}\left|\Phi^{[d]}_{ki}\right|\omega^{-2d}=\delta_{ki}, \end{align}

with ${\delta}_{ki}$ being the Kronecker delta function.

Figure 2. Topological localisation of network responses. For frequencies larger than all eigenfrequencies and across network types (row 1), the response amplitudes (4.13) decays exponentially with shortest-path distance d (row 2) and algebraically with driving frequency $\omega$ (row 3) (cf. Proposition 4.4). Dashed vertical lines in row 3 indicate the displayed frequency responses in row 2. Columns display graphs and responses for a random tree (column a), the topology of the British high voltage transmission grid (column b) and a random power grid network topology generated according to [31] (column c). Network settings: $\left(N,N_g,P_g,P_c,K_g,K_c,\alpha\right)=\left(264,24,10\text{ s}^{-2},-1\text{ s}^{-2},200\text{ s}^{-2},20\text{ s}^{-2},1\text{ s}^{-1}\right)$ for column a, $\left(120,30,39\text{ s}^{-2},-13\text{ s}^{-2},390\text{ s}^{-2},390\text{ s}^{-2},1\text{ s}^{-1}\right)$ for column b, and $(80,20,39\text{ s}^{-2},-13\text{ s}^{-2},390\text{ s}^{-2},$ $390\text{ s}^{-2},1\text{ s}^{-1})$ for column c ² .

Remark 4.6 (Localised response patterns in networks of multi-dimensional dynamical systems) We consider networks of N diffusively coupled identical units governed by n-dimensional dynamics. Each unit i has n state variables $\left(x_i^{[0]},x_i^{[1]}\cdots,x_i^{[n-1]}\right)$ and is governed by dynamics

\begin{align} \begin{cases} \dot{x}_i^{[0]}&=x_i^{[1]}\nonumber\\[3pt] \dot{x}_i^{[1]}&=x_i^{[2]}\nonumber\\ \cdots\nonumber\\ \dot{x}_i^{[n-1]}&=f\left(x_i^{[0]},x_i^{[1]}\cdots,x_i^{[n-1]}\right)+\displaystyle\sum_{j=1}^N g\left(x_j^{[0]}-x_i^{[0]}\right),\nonumber \end{cases} \end{align}

where $f:\mathbb{R}^N\rightarrow\mathbb{R}$ and $g:\mathbb{R}\rightarrow\mathbb{R}$ are functions, respectively, representing the intrinsic and the coupling dynamics of units and allowing for a stable fixed point of the system. If unit k is sinusoidally driven with frequency $\omega\rightarrow\infty$ , we conjecture that the amplitude $\tilde{A}_{i,m}^{(k)}$ of the sinusoidal response in state variable $D_t^m x_i$ of unit i is given by

(4.33) \begin{equation} \tilde{\!A}_{i,m}^{(k)}\overset{\omega\rightarrow \infty}{\sim}\left|\Psi^{[d]}_{ki}\right|\omega^{-n(d+1)+m},\end{equation}

where $|\Psi^{[d]}_{ki}|$ is a distance- and node-specific prefactor but independent on the perturbation frequency.

The response of a sinusoidally driven damped harmonic oscillator with $\ \tilde{\!A}_{i,m}^{(k)}\overset{\omega\rightarrow \infty}{\sim}\omega^{-2}$ can be seen as a special case of (4.33) with $n=2,d=0$ and $m=0$ . For networks of Kuramoto phase oscillators with $n=1$ , (4.33) is proven to be valid [see the Supplementary of [45]]. For the oscillator model of power grid with $n=2$ and $m=1$ , (4.33) reduces to Proposition 4.4.

Remark 4.7 (Generalisability of the steady-state response patterns in three frequency regimes) In the above discussions of the steady responses patterns in three frequency regimes in Section 4.1, we do not make any assumptions on the network topology. Therefore, our results on the characteristics of the homogeneous, the resonance and the localised response patterns in three regimes hold for arbitrary network topologies. Nevertheless, the evaluation of the parameters and the prefactors, such as $\Phi_{ki}^{[d]}$ in (4.14), is network and topology dependent by definition.

4.2 Topological factor of transient spreading dynamics

In the following, we focus on the transient response of AC power grids to external perturbations, which is referred to as the network responses close to the time of perturbation and thus describes the spatiotemporal pattern in the perturbation spreading process across the network. Particularly, we demonstrate how to extract the role of the network topology in the spreading pattern based on the LRT of the oscillator model.

To investigate the transient response close to the onset of perturbation at $t=0$ , we Taylor-expand the linear response (3.13) at node i to a sinusoidal perturbation at node k in powers of t as

(4.34) \begin{equation} \Theta_i^{(k)}(t)=\sum_{n=0}^{\infty}\dfrac{D_t^n{\Theta}_i^{(k)}(0)}{n!}t^n,\end{equation}

around $t=0$ , which is characterised by the time derivatives of the linear response at $t=0$ . Here, $D^n_t:=\frac{\mathrm{d}^n}{\mathrm{d}t^n}$ is Euler’s notation for differential operator. The n-th order time derivative of the linear response at $t=0$ is

(4.35) \begin{equation} D^n_t\Theta_i^{(k)}(0)=\sum_{\ell=0}^{N-1} \dfrac{v_k^{[\ell]}v_i^{[\ell]}\varepsilon e^{\imath\varphi}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}} \left[\dfrac{\left(\Delta^{[\ell]}_+\right)^n\left(\Delta^{[\ell]}_--\imath\omega\right)-\left(\Delta^{[\ell]}_-\right)^n\left(\Delta^{[\ell]}_+-\imath\omega\right)}{2\eta^{[\ell]}}+(\imath\omega)^n\right].\end{equation}

The transient response of a power grid network can thus be estimated by the first non-zero term in the power series of t (4.34). Interestingly, the resulting series does not start with low powers of t such as $t^0$ or $t^1$ as typical for common Taylor expansions. Instead, it typically starts with large powers of t as the following proposition illustrates.

Proposition 4.5 (Leading-term approximation of transient response) The transient response at node i in an AC power grid network (3.8) to a sinusoidal perturbation of frequency $\omega$ at node k with an onset at $t=0$ is approximated by the $(2d+2)$ -nd term in the Taylor expansion of the linear response (3.13) around $t=0$ ,

(4.36) \begin{equation} \Theta_i^{(k)}(t)=\dfrac{\varepsilon e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}{\left(2d+2\right)!}t^{2d+2}+O\left(t^{2d+3}\right).\end{equation}

Here $d:=d(k,i)$ denotes the graph-theoretic distance between node k and node i.

Proof. To find out the leading term in the Taylor series of the linear response (4.34), we study the summand of the $\ell$ -th eigenmode in the derivative of the linear response $D^n_t\Theta_i^{(k)}(0)$ (4.35), which we denote as $v_k^{[\ell]}v_i^{[\ell]}\varepsilon e^{\imath\varphi}F_n\left(\lambda^{[\ell]}\right)$ for convenience ³ . In the summand, the function $F_n\left(\lambda^{[\ell]}\right)$ appears to be a division of two polynomials of $\lambda^{[\ell]}$ . In fact, it can be shown that the leading term of $F_n\left(\lambda^{[\ell]}\right)$ (denoted as $\mathrm{LT}\left[{F_n\left(\lambda^{[\ell]}\right)}\right]\big)$ , i.e. the term with the highest order of $\lambda^{[\ell]}$ is

(4.37) \begin{align} \mathrm{LT}\left[{F_n\left(\lambda^{[\ell]}\right)}\right]=\left\lbrace \begin{array}{l@{\quad}l} (-1)^{\frac{n-1}{2}}\left(-\imath\omega+\dfrac{n-1}{2}\alpha\right)\left(\lambda^{[\ell]}\right)^{\frac{n-3}{2}}& \quad \text{if } n \text{ is odd,}\\[11pt] (-1)^{\frac{n-2}{2}}\left(\lambda^{[\ell]}\right)^{\frac{n-2}{2}}& \quad \text{if } n \text{ is even.} \end{array}\right.\end{align}

A derivation of the result (4.37) is given in Appendix B. We note that, for $n=0$ and $n=1$ , $F_n\left(\lambda^{[\ell]}\right)=0$ , which is a consequence of the choice of initial condition: the linear response and its first derivative are supposed to be zero at $t=0$ , as the responses $\Theta_i^{(k)}$ and the frequency response $\dot{\Theta}_i^{(k)}$ are zero at the onset of perturbation. For $n\geqslant 2$ , (4.37) indicates a monotonic relation between the degree of $F_n\left(\lambda^{[\ell]}\right)$ as a polynomial of $\lambda^{[\ell]}$ and the order of derivative n: $n=2\deg\left[F_n\left(\lambda^{[\ell]}\right)\right]\,{+}\,2$ .

As we have shown in Section 4.1.3, sums of the form of $\sum_{\ell=0}^{N-1}v_k^{[\ell]}v_i^{[\ell]}P^j\left(\lambda^{[\ell]}\right)$ , where $P^j\left(\lambda^{[\ell]}\right)$ represents a general polynomial of $\lambda^{[\ell]}$ of degree j, can be seen as $[P^j(\mathcal{L})]_{ki}$ , the ki-th element of the matrix $P^j(\mathcal{L})$ , a polynomial of $\mathcal{L}$ with degree j. Applying this result to the derivatives of the linear response (4.35), we find that $D^n_t\Theta_i^{(k)}(0)$ can be considered as the ki-th element of matrix $F_n(\mathcal{L})$ , a polynomial of $\mathcal{L}$ . The leading term of $D^n_t\Theta_i^{(k)}(0)$ is thus given by

(4.38) \begin{equation} \mathrm{LT}\left[{D^n_t\Theta_i^{(k)}(0)}\right]=\varepsilon e^{\imath\varphi} \mathrm{LT}\left[{F_n(\mathcal{L})}\right]_{ki}, \end{equation}

which contains $(\mathcal{L}^m)_{ki}$ with $m=\frac{n-3}{2}$ if n is odd and $m=\frac{n-2}{2}$ if n is even. We notice that for a given node pair (k, i) at distance d, we have $(\mathcal{L}^m)_{ki}=0$ for all $m<d$ because no path of length $m<d$ can connect nodes k and i. Therefore, all terms in the Tayler series ( $4.34$ ) with the leading term’s degree lower than d vanish. The first non-zero term in the series thus equals the leading term of the $(2d+2)$ -th derivative of the linear response,

(4.39) \begin{equation} D^{2d+2}_t\Theta_i^{(k)}(0)=\varepsilon e^{\imath\varphi} \mathrm{LT}\left[{F_{2d+2}(\mathcal{L})}\right]_{ki}=\varepsilon e^{\imath\varphi}(-1)^d\left(\mathcal{L}^d\right)_{ki}, \end{equation}

because all other terms contain $(\mathcal{L}^m)_{ki}$ with $m<d$ and thus vanish. Taken together, the transient linear response near $t=0$ can be approximated as

(4.40) \begin{equation} \Theta_i^{(k)}(t)=\sum_{n=2d+2}^{\infty}\dfrac{D_t^n\Theta_i^{(k)}(0)}{n!}t^n=\dfrac{\varepsilon e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}{\left(2d+2\right)!}t^{2d+2}+O\left(t^{2d+3}\right). \end{equation}

Remark 4.8 (Topological factor in perturbation spreading) The leading-term approximation of the linear response (Proposition 4.5) provides a way to disentangle the impact of various factors on the dynamical spreading process in power grid networks. Specifically, the impact of the specific network topology, including the interaction structure between units and the system’s base state at $t=0$ , is reflected in the factor $\left(\mathcal{L}^d\right)_{ki}$ in the leading-term approximation. It satisfies

(4.41) \begin{equation} \left(\mathcal{L}^{d}\right)_{ki}=\sum_{\mathcal{P}_{k\rightarrow i}^{d}}\ \prod_{(u,v)\in\mathcal{P}_{k\rightarrow i}^{d}}\mathcal{L}_{uv}, \end{equation}

suggesting that it can be interpreted as the product of the edge weights along a shortest path $\mathcal{P}_{k\rightarrow i}^{d}$ between node k and i, summed over all shortest paths. This insight provides guidelines for manipulating the perturbation spreading dynamics in power grid networks through changing the underlying topology. Numerical evidences show that the topological factor that revealed by the leading-term approximation also enables a master function approach to accurately predict threshold-crossing arrival times in power grid networks [47].

Remark 4.9 (Scaling behaviours in transient responses) The leading-term approximation of the transient response (4.36) reveals two scaling behaviours as $t\rightarrow0$ (cf. Figure 3). First, the transient response grows algebraically in time with a distance-dependent exponent: $\Theta_i^{(k)}(t)\sim C_d t^{2d+2}$ . Here, $C_d:={\varepsilon e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}/{\left(2d+2\right)!}$ is a time independent prefactor but depends on signal magnitude, topology, base operating state and inter-node distance d. Second, the transient response decays nearly exponentially with distance d, since the factor $t^{2d+2}$ dominates the asymptotic behaviour of the response at large but finite distances as $t\rightarrow0$ .

Figure 3. Transient Network Response Dynamics exhibits algebraic growth with time and exponential decay with shortest-path distance. (a) Basic network of $N=6$ units illustrates (b,c) transient algebraic responses [colour-coded as units in (a)] to a sinusoidal perturbation at node 1 that increase like $\Theta_i^{(k)}(t)\sim C_d t^{2d+2}$ as $t \rightarrow 0 $ , with time independent constant $C_d$ that depends on signal magnitude, topology, base operating state and inter-node distance d, see (4.36). Thus, responses (b) algebraically increase with time t at any given unit and (c) at any given time, they decay nearly exponentially with shortest-path distance $d=d(k,i)$ between the perturbed unit k and the observed unit i. The grey dotted lines in (b) indicate the leading-term approximations. Network settings: $\left(N,N_g,P_g,P_c,K_g,K_c,\alpha\right)=(6,3,1\text{ s}^{-2},-1\text{ s}^{-2},10\text{ s}^{-2},10\text{ s}^{-2},1\text{ s}^{-1})$ . For the perturbation signal $(\varepsilon,\omega/2\pi,\varphi)=(1,1\text{ Hz},0\text{ rad})$ .

Remark 4.10 (Generalisability of the transient spreading patterns) The above discussions on the spatio-temporal pattern of transient spreading do not involve any assumptions on the underlying network topology. Thus, the form of the leading term approximation of the transient network response (Proposition 4.5) does not depend on the specific choice of network topology. We underline that the evaluation of the topological factor in perturbation spreading (Remark4.8) is indeed topology- and node-specific, as it captures the local interaction structure and consists of all shortest paths between the perturbed node and the responding node.

4.3 Nodal vulnerability to unpredictable fluctuations

In Sections 4.1 and 4.2, we show how the distributed response patterns of power grid systems, in the steady state ( $t\rightarrow\infty$ ) and in the transient stage ( $t\rightarrow0$ ), are analytically extracted from the LRT of the oscillator model. The response patterns are numerically proven to be highly accurate [45, 47], but the results are valid only for given perturbation signals, hence deterministic. Meanwhile, real-world power grid systems are perturbed by power fluctuations whose exact time series are hardly predictable. In this section, we discuss how LRT helps to estimate network responses to random perturbations.

As discussed in Section 4.1, power grid systems respond resonantly to perturbations with frequencies falling in the band of network’s eigenfrequencies $I_{\text{res}}:=\left[\omega_{\text{eigen}}^{[1]},\omega_{\text{eigen}}^{[N-1]}\right]$ , exhibiting the most irregular network-wide spatiotemporal patterns, compared to the almost homogeneous pattern for lower frequencies and the localised pattern for higher frequencies. Moreover, the resonance response pattern varies drastically for different perturbation frequencies and for different locations of perturbation. Therefore, estimating the nodal responses for random perturbation signals involving frequency components within $I_{\text{res}}$ is a task not only of practical significance regarding the operational safety of power grid systems, but also with a high theoretical complexity.

Definition 4.1 (Dynamic vulnerability index (DVI) for random network resonances) In an AC power grid system (3.8) perturbed by a random fluctuation at node k, characterised by a power spectral density $S(\omega)$ with frequency components $\omega\in I_{\mathrm{res}}=\left[\omega_{\mathrm{eigen}}^{[1]},\omega_{\mathrm{eigen}}^{[N-1]}\right]$ , the nodal Dynamic Vulnerability Index (DVI) is defined as

(4.42) \begin{equation} \mathrm{DVI}_i^{(k)}:=\int_{I_{\mathrm{res}}}{S(\omega)}^{\frac{1}{2}}\left|\sum_{\ell=0}^{N-1}\dfrac{\imath\omega v^{[\ell]}_kv^{[\ell]}_i}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}\right|\mathrm{d}\omega.\end{equation}

Proposition 4.6 (DVI estimates ranking of the nodal all-time-high steady frequency responses) Let an AC power grid system (3.8) be perturbed by a time-dependent fluctuation at node k. A random signal time series is characterised by a power spectral density $S(\omega)$ , i.e., the strength $\varepsilon$ of its frequency component $\varepsilon e^{\imath(\omega t+\varphi)}$ is frequency dependent and follows $\varepsilon(\omega)\propto S(\omega)^{\frac{1}{2}}$ , and the corresponding phase $\varphi$ is randomly drawn from the uniform distribution on $[0,2\pi)$ , independently for each realisation of the fluctuation time series. Suppose the fluctuation signal is composed of frequency components with $\omega\in I_{\mathrm{res}}=\left[\omega_{\mathrm{eigen}}^{[1]},\omega_{\mathrm{eigen}}^{[N-1]}\right]$ , the ranking $\sigma_{\mathrm{ATH}}$ of the nodal all-time-high steady frequency response magnitude $\max_{t\in[0,T]}|\dot{\Theta}^{(k)}_i(t)|$ in an observation window T approaches the ranking $\sigma_{\mathrm{DVI}}$ of the DVI defined in Definition 4.1 as the observation time window goes to infinity, i.e.

(4.43) \begin{equation} \lim_{T\rightarrow\infty}\sigma_{\mathrm{ATH}}(i)=\sigma_{\mathrm{DVI}}(i) \quad \text{for all }i\in\{1,\cdots,N\}. \end{equation}

Proof. To analyse the steady network response, i.e. the response in a time window $T\rightarrow\infty$ , to a perturbation signal composed of a range of frequencies, we make use of the steady-state response of the oscillator model to a sinusoidal signal (4.6). The steady-state response $\dot{\Theta}_i(t)$ of the frequency at node i to a perturbation signal $\varepsilon e^{\imath(\omega t+\varphi)}$ at node k is given by

(4.44) \begin{equation} \dot{\Theta}_i^{(k)}(t)\overset{t\rightarrow\infty}{\sim}\sum_{\ell=0}^{N-1}\dfrac{\imath\varepsilon\omega v_k^{[\ell]}{v}_i^{[\ell]}}{-\omega^2+\imath\alpha\omega+\lambda^{[\ell]}}e^{\imath(\omega t+\varphi)}=\imath\varepsilon\omega R_i^{(k)}e^{\imath(\omega t+\varphi)},\end{equation}

which can be seen as a driven oscillation with a complex amplitude $\imath\varepsilon\omega R_i^{(k)}$ . Here, $R_i^{(k)}$ is the response factor defined in (4.5) characterising the response strength at individual nodes in a network. The complex amplitude gives rise to an amplitude shift $\big|\imath\varepsilon\omega R_i^{(k)}\big|$ and a phase shift $\arg\left(\imath\varepsilon\omega R_i^{(k)}\right)$ , which both are specific to the perturbation strength $\varepsilon$ and the perturbation frequency $\omega$ .

As a consequence of the linear nature of LRT, the nodal frequency response to a temporally fluctuating perturbation signal containing a spectrum of frequency components is obtained by summing the response $\dot{\Theta}_i^{(k)}(\varepsilon,\omega,t)$ (4.44) over all frequency components (see also Remark 2.4). Particularly, for perturbation signals which are characterised by a specific power spectral density (PSD) S(w), the strength of its frequency component $\varepsilon e^{\imath(\omega t+\varphi)}$ can be expressed in terms of the frequency as $\varepsilon(\omega)\propto S(\omega)^{\frac{1}{2}}$ while no assumptions is made on the choice of its phase $\varphi$ . For instance, in modern power grids integrated with renewable energies, the power fluctuations from wind and solar energy are characterised by a power law PSD with the Kolmogorov exponent $-5/3$ [[2]]. For such resonant perturbations, the nodal frequency response reads

(4.45) \begin{equation} \dot{\Theta}_i^{(k)}(S(\omega),t)\overset{t\rightarrow\infty}{\sim}\displaystyle\int_{I_{\text{res}}}\imath cS(\omega)^{\frac{1}{2}}\omega R_i^{(k)}e^{\imath(\omega t+\varphi)}\mathrm{d}\omega. \end{equation}

Here, c is a scaling factor between the power distribution over frequencies of a specific signal and the (normalised) PSD. Independent of the specific realisation of the fluctuation signal, the largest possible magnitude of the frequency response (4.45) is reached only when the involved frequencies are finite in number and non-resonant to each other ⁴ so that all oscillating responses for each frequency component $\omega$ (4.44) would eventually align. The largest possible magnitude of frequency response time series is given by the sum over $\omega$ of the magnitude of each frequency response $|\imath cS(\omega)^{\frac{1}{2}}\omega R_i^{(k)}|$ . For finite time series with M data points and time step $\Delta T$ , its frequency components $\omega_n=\frac{2\pi n}{M\Delta T}$ with $n \in \left\{1,\cdots,\frac{M}{2}\right\}$ given by discrete Fourier transform are apparently resonant to each other. However, for time series with a fixed sampling rate, the longer the observation time window T, the smaller the frequency interval $\Delta \omega=2\pi/T$ and thus the more frequency components exist within the given interval of interest $I_{\text{res}}$ . As T approaches infinity, the order of the resonant frequencies $\{\omega_n\}$ becomes sufficiently high so that the alignment of the phases can be attained at a finite rate [8]. Therefore, as $T\rightarrow\infty$ , the all-time-high frequency response approaches the sum of the frequency-specific response magnitudes, which is proportional to the DVI given in (4.42):

(4.46) \begin{equation} \lim_{T\rightarrow\infty}\max_{t\in[0,T]}\left|\dot{\Theta}^{(k)}_i(S(\omega),t)\right|=\int_{I_{\text{res}}}cS(\omega)^{\frac{1}{2}}\left|\imath \omega R_i^{(k)}\right|\mathrm{d}\omega=c\,\text{DVI}_i^{(k)}.\end{equation}

If one only considers the relative ranking of the all-time-high frequency response of nodes in a given network, but not their absolute values, the overall scaling factor c in (4.46) does not play a role in the ranking and we finally arrive at

(4.47) \begin{equation} \lim_{T\rightarrow\infty}\sigma_{\mathrm{ATH}}(i)=\sigma_{\mathrm{DVI}}(i) \quad \text{for all }i\in\{1,\cdots,N\}.\end{equation}

Remark 4.11 (Generalisation of DVI and convergence time) The DVI defined in (4.42) provides a measure to estimate the relative nodal risk from the power grid network resonances induced by a unpredictable perturbation signal and thus helps to identify the most vulnerable nodes in networks with arbitrary topologies exhibiting particularly strong resonant responses. The integration interval of the frequencies in DVI is chosen to be the band of eigenfrequencies $I_{\text{res}}$ for a specific network so that its irregular resonance patterns are covered, however it is not a must. In principle, any frequency range can be chosen for DVI to estimate the relative strength of the all-time-high responses in the specific frequency range. However, one should note that the timescale for the ranking of all-time-high responses to converge to the ranking given by DVI (in 4.46 and 4.47) depends on the chosen frequency range. For instance, the convergence time would be longer if lower frequencies are included in the integration interval of DVI.

5 Role of LRT in uncovering response patterns

In the previous section, we elaborated how network-wide dynamic response patterns of power grid systems can be extracted from the explicit solution of nodal responses given by the LRT of the oscillator model (see Section 3). In this section, we summarise and compare the role of LRT in revealing different categories of the dynamic response patterns, such as the patterns emerging on different timescales, and in responses to perturbations with different levels of randomness and magnitude.

5.1 Transient versus steady-state responses

As discussed in Sections 4.1 and 4.2, power grid transmission networks exhibit distinctive spatiotemporal response patterns on different timescales. The transient spreading pattern (see Section 4.2) of an external perturbation signal in a normally operating power grid network appears close to the time of impact $t=0$ . It is characterised by a set of points in time, at which the impact arrives at individual units in the network. To a large extent, the topological dependence of the arrival times can be captured by a topological factor which arises from the leading-term approximation of the linear response. The steady-state response patterns, in contrast, emerge as $t\rightarrow\infty$ (see Section 4.1), where the nodal responses to a sinusoidal perturbation converge to sinusoidal oscillations as well, but with various amplitudes. Consequently the set of the nodal response amplitudes, a time-invariant but frequency-dependent feature of the oscillating responses, constitute the steady-state response patterns characterising three frequency regimes. The time scale separating transient from steady-state regimes is not a universal constant but intricately depends on several factors including network size N and network topology, damping constant $\alpha$ , and the specific node location we are interested in within the network.

Both transient and steady-state response patterns have been revealed and characterised through asymptotic analyses of the explicit solution of the linear nodal responses (3.13). The solution depends explicitly on time while the dependence on the network topology is implicit, embedded in the eigenvalues and eigenvectors of the weighted graph Laplacian matrix $\mathcal{L}$ . Through asymptotic analyses, either with respect to time t or the perturbation frequency $\omega$ , the lengthy solution of the linear nodal response is reduced to one term per eigenmode that dominates the asymptotic behaviour as the variable t or $\omega$ approaches its corresponding limit (see Proposition 4.4 and Proposition 4.5). As the contribution of each eigenmode contains the ‘overlap factor’ $v_i^{[\ell]}v_k^{[\ell]}$ of the perturbed node k and the responding node i, the powers of the Laplacian eigenvalue $\left(\lambda^{[\ell]}\right)^m$ that is involved in the dominating term in each eigenmode $\ell$ translate to elements of the power of the Laplacian matrix $(\mathcal{L}^m)_{ki}$ through the summation over all N eigenmodes $\ell\in\{0,\cdots,N-1\}$ . Furthermore, with the help of the result from graph theory that $(\mathcal{L}^m)_{ki}\equiv0$ if $m\in\mathbb{N}_0$ is smaller than d(k, i), the shortest path distance between node k and i, the dependence of the nodal linear responses on graph-theoretic distance emerges. In this way one obtains the asymptotic spatiotemporal response patterns that depends explicitly on distance.

A major difference between the patterns we uncover in steady-state responses and in transient responses is, the asymptotic behaviours of the steady-state response patterns are exact, in the sense that the higher order terms are negligible at the observed limits of the perturbation frequency $\omega$ (see Section 4.1). Meanwhile, the contribution of the higher order terms are not negligible in the patterns in transient responses: the leading-term approximation of the response exhibits a diverging error as the perturbation spreads further and the arrival time grows larger. Numerical simulations show that the higher order terms accounts for about 10% of the actual arrival time of perturbations [47], which is significant in predicting the perturbation spreading behaviour in real-world power grid systems. However, by means of numerical techniques, we can still use the topological factor proposed in Remark 4.8 to estimate the contribution of higher order terms $O(2d+3)$ in the Taylor series (4.36) in a specific network ensemble and give accurate predictions for the actual arrival times of the impact of a perturbation [47]. Related recent works [43, 30] studied transient propagation of perturbations in networked systems consisting of one-dimensional dynamical units. One main finding is a similar scaling of the unit’s state variables $x_i(t)$ (or their deviations from a base state) with time t as $x_i(t)\sim t^d$ where d is the shortest-path distance between perturbed node and the node i the response in measured at.

5.2 Responses to deterministic perturbations versus responses to unpredictable perturbations

The power grid response patterns we discuss in this article can be classified into two categories: the ones emerging in the deterministic responses to a given perturbation signal, such as the transient responses (Section 4.2) and the steady-state responses (Section 4.1) to a given signal, and the ones that estimate the cumulative impact of an unpredictable signal on a power grid network, such as the DVI measuring the nodal risk of network resonances (Section 4.3). Both categories of power grid response patterns are discovered based on the explicit solution of the linear nodal responses (3.13) that derived from the LRT.

Looking more closely, one finds that the estimated patterns in the cumulative responses can be seen as a result built upon the deterministic steady-state linear responses with one more dimension, i.e. the specific properties of the perturbation signals. At its core, the steady-state network frequency responses $\dot{\boldsymbol{\Theta}}$ (4.44) to a single-frequency sinusoidal oscillation $D_k(t)$ at a given node k is a mapping $f_{G,\boldsymbol{\theta}^*}:\mathbb{R}\rightarrow\mathbb{R}^N$ with $f_{G,\boldsymbol{\theta}^*}$ depending on the underlying network topology G and network’s base state $\boldsymbol{\theta}^*$ prior to perturbations. The fluctuating nature of the perturbation signal adds another dimension to the responses: the amplitude of $D_k$ is no longer a given constant $\varepsilon$ , but becomes further dependent of the frequency $\omega$ through the PSD $S(\omega)$ of the signal, i.e. $\varepsilon(\omega)=S(\omega)^{\frac{1}{2}}$ . In this way, the unknown or unpredictable temporally detailed features of the perturbations are integrated into the LRT framework for networked dynamical systems, which extends standard LRT statements and enables us to estimate features of network responses beyond the deterministic realm of systems driven with known signals.

We emphasise that, due to the intrinsic irregularity of fluctuating perturbation signals, the errors of the estimates for the associated responses appear to be significantly higher than the ones for the deterministic responses [compare [45] and [46]]. In a finite signal time series characterised by a PSD function, the contained frequencies are finite in number in the considered frequency interval (such as the resonance regime $I_{\text{res}}$ ), and apparently resonant to each other, which leads to a deviation in the ranking of the all-time-high nodal responses to the ranking given by indices computed a priori (such as the DVI discussed in Section 4.3). Additionally, realistic perturbation signals do not follow exactly the characteristic PSD, such that randomness also exists in the amplitudes of the frequency components. Nevertheless, compared to the deterministic patterns which gives only a posteriori information of network responses, estimates such as the DVI given by the extended LRT may provide a useful guiding tool for risk assessments in real-world power grid systems.

5.3 Small responses versus large responses

As the name suggests, LRT provides the linear approximation of a system’s response to a perturbation close to a fixed point of a networked dynamical system. Therefore, the solution given by LRT intrinsically deviates from the actual system responses due to the neglected higher order terms in the system’s collective nonlinear dynamics. As the system being driven further and further away from the fixed point, the responses typically increase and so do the estimation error of the LRT. However, the range of validity of the LRT, as well as how its error grows with the perturbation, is usually nontrivial and system-dependent.

For the oscillator model of power grid networks, the error of the solution given by LRT (3.12 and 3.13) follows the same trend and grows with an increasingly stronger perturbation signal. However, numerical evidence shows that the error increases mildly with the magnitude of the perturbation until it blows up close to a bifurcation point [45]. In the linearised dynamics of the oscillator model at a fixed point $\boldsymbol{\theta}^*$ (3.14), the deviation of the nonlinear coupling terms $\sin(\theta_j-\theta_i)$ for all edges $(i, j)\in E$ to their values $\sin\left(\theta^*_j-\theta^*_i\right)$ at the fixed point are represented by the first-order approximations $\cos\left(\theta^*_j-\theta^*_i\right)\left(\Theta_j-\Theta_i\right)$ , vectorised as the term $-\mathcal{L}\boldsymbol{\Theta}$ in (3.14). For power grid systems working at a stable operation state without any transmission line overloaded ( $|\theta_j^*-\theta_i^*|\leqslant \frac{\pi}{2}$ for all edges $(i, j)\in E$ , see Proposition 2.1), the linear approximation breaks down only when the system is driven far enough from the fixed point $\boldsymbol{\theta}^*$ and goes close to the point where one of the lines (i, j) is fully loaded, i.e. $\sin(\theta_j-\theta_i)=1$ . In this regime, the linear approximation diverges from the actual state of the system and the error grows explosively. As elaborated in an article by [22], when one of the transmission line is fully loaded, the power grid system reaches a bifurcation point where the initial stable fixed point is lost and the oscillating units in the system becomes desynchronised. Therefore, the LRT of the oscillator model of the power grid systems, together with all of the derived response patterns, are generally valid as long as none of the lines become overloaded and the entire system becomes unstable [see [45] for quantitative results of the LRT errors].

6 Conclusions and outlook

6.1 Conclusions and discussions

In this work, we systematically discuss how LRT may shed light on the spatio-temporal response patterns emerging in networked dynamical systems under time-dependent perturbations. We exemplify a full analysis for model dynamics of power grid systems which are inevitably exposed to fluctuating power injections from renewable energy sources. Beyond previous works, we integrate and present all details required for a full mathematical analysis, specifically demonstrate how to evaluate the generally intricate, multiple-sum expressions determining spatio-temporal response patterns in a useful way and highlight how different results interconnect, for instance between transient and long-term dynamics or between different types of perturbations and across topologies. We introduce the main ideas of LRT and its general requirements for applicability on system settings (cf. Section 2). We explicate various aspects of application to models of power grid systems, such as i) the solution of linear responses of the stationary DC power flow model and of the dynamic oscillator model of AC power grids (cf. Section 3), and ii) how it helps to identify dynamic patterns in network-wide responses which provide theoretical guidelines for power grid design, control and risk assessments.

Although LRT has been widely used as a powerful tool in analysing various response dynamics of many complex networked dynamical systems, the works presented in this article provide a fresh methodological angle to approach the problem. For power grid systems, LRT has been used to estimate, e.g. quadratic performance measures for the network’s overall excursion away from synchrony [36, 37, 26, 6] and the variance of the frequency response increment distribution [11]. For connectome dynamics in brain, LRT has also been used to analytically estimate the covariance of Gaussian linear model of the stochastically perturbed system [35], which links the structural and the functional connectivities between brain regions [44, 39].

The work presented in this article approaches the dynamic network responses in a way different from the above-mentioned works: instead of quantifying the stochastic features of the overall or distance-specific responses directly based on the linear responses, we start from explicating the deterministic solution of network-wide responses to a single-frequency signal and using methods from graph theory and asymptotic analysis to extract spatiotemporal response patterns. These may be interpreted in a physically intuitive way. Specifically, the three frequency regimes of steady-state response patterns (Section 4.1) and the master curve of transient perturbation spreading (Section 4.2) are entirely deterministic. Especially, the former work on fluctuation-induced network resonances can be seen as a direct generalisation of the classical resonance phenomenon of a single driven damped harmonic oscillator to oscillators interacting on networks. The emerging pattern constituted by the estimations of the all-time-high nodal response magnitudes to a irregularly varying perturbation (Section 4.3), is also a straightforward result derived from the spatial patterns in the network-wide responses to a deterministic (periodic) signal.

6.2 Challenges and future work

Future work regarding the response theory for networked dynamical systems may follow several directions.

First, exact analytical solutions, and even many asymptotic results of the linear responses of general networked dynamical systems, as well as the response patterns emerging from these solutions, remain unknown to date. Applying the LRT presented in this work on a networked dynamical system, we employed several conditions on the system’s dynamics (see Section 2) to explicate a full analysis without too many notational and other complications. The conditions include homogeneous nodal dynamics, a diffusive coupling term $g_{ij}(x_j-x_i)$ and the evenness of the coupling function’s sensitivity $\frac{\mathrm{d} g_{ij}}{\mathrm{d}(x_j-x_i)}$ to small changes in the difference of nodal states $x_j-x_i$ , such that a symmetric weighted Laplacian matrix arises in the linearised response dynamics of the system at the fixed (operating) point (2.10). The presence of a symmetric Laplacian matrix in the linearised dynamics ensures the option to express the linear responses in the Laplacian eigenbasis, which plays a critical role in linking the response at a specific node to the graph-theoretic distance to the perturbation. Thereby it is also critical in uncovering the topological structure of the dynamical response patterns across the network. However, for many networked dynamical systems, such as the third-order model of power grid dynamics including the voltage dynamics [19], such preconditions are not fulfilled. One way to overcome this theoretical barrier and to extend LRT to such networks of dynamical systems is to transform the system’s state variables to another coordinate system where the Jacobian matrix $\mathcal{J}$ in (2.3) is diagonal or almost diagonal (such as in the Jordan normal forms of $\mathcal{J}$ ). In this way, explicit solutions for linear responses can be obtained in the new coordinate system where dynamic response patterns can be identified in similar ways presented in this work.

Second, one could use LRT to develop strategies to control the impact of fluctuations on networked dynamical systems such as power grids. So far, we gained insights into the spatiotemporal structure of the responses across networks and developed indices to estimate the nodal risks against external fluctuations. The next step towards more reliable and more robust power grid systems would be to utilise the obtained understanding to develop countermeasures against the risks, e.g. to suppress the potentially dangerous responses such as network resonances and to slow down the spreading of the impact of a sudden drop of injected power. A potential way to achieve such tasks could be to manipulate discovered response patterns by changing the interaction structure of the power grid.

Third, research on how different classes of network topologies potentially impact response patterns through their specific characteristics of their eigensystems. Progress in this direction seems hard, because one would need to be able to characterise, e.g., eigenvectors of graph ensembles such that they directly help to extract useful information from complex expressions like (4.5) or even (3.13) specifically for that ensemble.

Fourth, the LRT per se could be extended by considering also the higher-order approximations of the system’s responses, cf. [34]. In the current work, we demonstrated that many features of collective response patterns of networked dynamical systems with nonlinear couplings, such as the oscillator model of power grids, are dominantly captured by the first-order (i.e. linear) approximation of the system’s responses, yet nonlinear effects may also play a role for other systems with certain forms of the intrinsic nodal dynamics or the interaction dynamics between nodes, specifically if we ask for the loss of solutions near operating states [34]. Therefore, it would be desirable if the contributions of the higher-order approximations of the system’s responses can be estimated. An open question also here is how such nonlinear effects depend on the interaction topology of the network.

We conjecture that general network dynamical systems, also beyond power grids, similarly respond in characteristic ways to external input signals, making the systems non-equilibrium and often non-stationary, and to be described by non-autonomous deterministic or stochastic evolution equations. Several of the analysis steps presented above hint that the key methodological tools are either readily transferable to more general systems’ settings or may be adapted to such settings. Candidate classes of systems include networks of multi-dimensional units, with discrete or hybrid dynamics, with delayed interactions or with spatially or temporally correlated stochastic inputs. Application areas may range from gene regulatory networks and metabolic circuits in cell biology to the controlled self-organised dynamics of engineered systems with feedback, from complex mechatronic systems to swarms of autonomous aerial vehicles.