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$ .

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 [Reference Hale12, Reference Larter, Speelman and Worth18, Reference Postnov, Ryazanova, Mosekilde and Sosnovtseva27, Reference Casagrande5, Reference Stankovski, Pereira, McClintock and Stefanovska32], in particular also appearing in the Kuramoto model [Reference Kuramoto17] and its variations [Reference Filatrella, Nielsen and Pedersen10, Reference Acebrón, Bonilla, Pérez Vicente, Ritort and Spigler1]. 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}$ .

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 [Reference Bapat3]. 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.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 [Reference Milano23]. By linearising the AC power flow equations at an operation point, the DC power flow modelFootnote ¹ 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 [Reference Machowski, Bialek and Bumby19].

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

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}$ .

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 [Reference Filatrella, Nielsen and Pedersen10, Reference Rohden, Sorge, Timme and Witthaut28] and [Reference Motter, Myers, Anghel and Nishikawa24]. This model allows for analytical understanding of the dynamics of power grids and has yielded fruitful research results over the past decade [Reference Rohden, Sorge, Timme and Witthaut28, Reference Motter, Myers, Anghel and Nishikawa24, Reference Dörfler, Chertkov and Bullo7, Reference Witthaut, Rohden, Zhang, Hallerberg and Timme42, Reference Tyloo, Pagnier and Jacquod37]. 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 [Reference Kundur, Balu and Lauby16, Reference Machowski, Bialek and Bumby19]:

(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}$ .

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. [Reference Manik, Witthaut, Schäfer, Matthiae, Sorge, Rohden, Katifori and Timme22]]. 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.

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}

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).

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.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 [Reference Schultz, Heitzig and Kurths31] (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 cFootnote ² .

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 [Reference Zhang, Hallerberg, Matthiae, Witthaut and Timme45]]. 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 convenienceFootnote ³ . 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 [Reference Zhang, Witthaut and Timme47].

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})$ .

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 [Reference Zhang, Hallerberg, Matthiae, Witthaut and Timme45, Reference Zhang, Witthaut and Timme47], 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$ [[Reference Anvari, Lohmann, Wächter, Milan, Lorenz, Heinemann, Tabar and Peinke2]]. 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 otherFootnote ⁴ 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 [Reference Dumas8]. 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}

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 [Reference Zhang, Witthaut and Timme47], 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 [Reference Zhang, Witthaut and Timme47]. Related recent works [Reference Wolter, Lünsmann, Zhang, Schröder and Timme43, Reference Schröder, Zhang, Wolter and Timme30] 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 [Reference Zhang, Hallerberg, Matthiae, Witthaut and Timme45] and [Reference Zhang, Ma and Timme46]]. 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 [Reference Zhang, Hallerberg, Matthiae, Witthaut and Timme45]. 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 [Reference Manik, Witthaut, Schäfer, Matthiae, Sorge, Rohden, Katifori and Timme22], 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 [Reference Zhang, Hallerberg, Matthiae, Witthaut and Timme45] for quantitative results of the LRT errors].