## Dynamic Vulnerability in Oscillatory Networks and Power Grids

Xiaozhu Zhang, 1, ∗ Cheng Ma, 2 and Marc Timme 1, †

1 Chair for Network Dynamics, Institute for Theoretical Physics and Center for Advancing Electronics Dresden (cfaed),

Technical University of Dresden, 01062 Dresden, Germany

2 School of Physics, Nankai University, Tianjin 300071, China

(Dated: August 5, 2019)

Recent work found distributed resonances in driven oscillator networks and AC power grids. The emerging dynamic resonance patterns are highly heterogeneous and nontrivial, depending jointly on the driving frequency, the interaction topology of the network and the node or nodes driven. Identifying which nodes are most susceptible to dynamic driving and may thus make the system as a whole vulnerable to external input signals, however, remains a challenge. Here we propose an easy-to-compute Dynamic Vulnerability Index (DVI) for identifying those nodes that exhibit largest amplitude responses to dynamic driving signals with given power spectra and thus are most vulnerable. The DVI is based on linear response theory, as such generic, and enables robust predictions. It thus shows potential for a wide range of applications across dynamically driven networks, for instance for identifying the vulnerable nodes in power grids driven by fluctuating inputs from renewable energy sources and fluctuating power output to households.

## I. INTRODUCTION

Oscillatory networks, modeling the underlying mechanisms of many real-world systems ranging from gene and neural circuits [1, 2] to AC power grids [3-8], exhibit highly nontrivial responses to external driving signals [913], due to the complexity in the underlying topology and the nonlinearity in the coupling function. Recently, growing attention has been drawn to the topic of dynamically driven networks, in part because of the important application of the second-order Kuramoto-type oscillator model in power grid operation and control [12-16]. With an increasing share of fluctuating renewable energy sources integrated in modern power grids, it is crucial for grid operators to predict the distributed frequency responses to systematic and stochastic fluctuations, to identify which units are most susceptible and may thus make the system as a whole vulnerable to dynamic inputs.

[22] and estimated via the nodal variance [23]. Yet, it is still unclear which stochastic signals may cause networkwide response patterns and how to quickly and precisely identify those nodes that potentially exhibit most severe responses and thus are most vulnerable to such perturbations.

For the example of power grid models, key aspects of network responses to dynamical perturbations have been uncovered recently: about the impact of various types of perturbation signals, including the scaling in the relaxation of power grids after pulse-like perturbations [15, 16], the differential response to static perturbations[17], the distributed dynamic patterns in response to dynamic perturbations [12], the fluctuationinduced non-Gaussian grid frequency distribution [18] and the escape of a system from an operation state if driven by white noise [14] and/or non-Gaussian noises [19, 20]. Specifically, the time-averaged nodal deviations from the network mean response was ranked using a centrality measure based on the Laplacian spectrum [21]. For lossy networks, averaged nodal sensitivity to fluctuations across network have been numerically investigated

∗ xiaozhu.zhang@tu-dresden.de

† marc.timme@tu-dresden.de

The core of the puzzle lies in an intriguing phenomena of dynamic network resonance[12], which is present in oscillator models with two (or more) variables per node (such as the second order Kuramoto model, an oscillatory power grid model) but not in the networks of phase oscillators such as in the original Kuramoto model [9-12]. While the network responses for low- and high-frequency signals are trivial thus fairly predictable-homogeneous responses for slowly-changing signals and localized responses for fast-changing signals-fluctuations in the resonance frequency regime of a network system induce complex resonance patterns in oscillatory networks [12]. The patterns are jointly determined by the perturbation frequency, the underlying network topology, the initial unperturbed network state (base operating state), and the location of the perturbation and the response of interest. Although the resonant responses can be deterministically and precisely computed for given perturbation time series deriving and evaluating a linear response theory [12], a straightforward, fast and reliable method for estimating the resonant response strengths for stochastic signals is still missing, mainly because such signals contain an extended band of frequencies.

Here, going beyond structural vulnerability in networks, we propose the Dynamic Vulnerability Index (DVI), a computationally inexpensive vulnerability measure to assess and to rank the largest possible resonant response of individual nodes in oscillatory networks. The networks are driven by stochastic perturbations containing a characteristic power spectral density (PSD) function. In power grid research and beyond, the term network vulnerability is typically used to describe the impact of purely topological changes on network performance

[24-28]. The meaning of the vulnerability of a node was extended to considering the node's transient response to a pulse-like perturbation [29, 30], and recently to the timeaveraged response to stochastic perturbations [21]. Here we propose a Dynamic Vulnerability Index (DVI) that expands the definition by considering the global maximum of a node's dynamic response to a stochastic input signal. Employing a linear response theory[12] and a frequencyspecific estimate of the resonant response strength, the DVI exhibits high prediction power and helps to identify those nodes potentially respond most strongly to a stochastic resonant perturbation and thereby posing systemic risks in power grid stability. Especially, the DVI identifies the vulnerable nodes at unexpected locations in the network not foreseeable from the topology alone.

## II. RESONANT NETWORK RESPONSE PATTERNS

Consider a network of N second-order Kuramoto-type oscillators with dynamics governed by

$$\ddot { \theta } _ { i } = P _ { i } - \alpha \dot { \theta } _ { i } + \sum _ { j = 1 } ^ { N } K _ { i j } \sin ( \theta _ { j } - \theta _ { i } ) + \delta _ { i k } D ( t ) \quad ( 1 ) \quad \begin{smallmatrix} \lambda _ { i } \end{smallmatrix} \ddot { \alpha } = \Omega \\ \quad \begin{smallmatrix} \theta _ { i } \end{smallmatrix} \\ \theta _ { \alpha \text{recoat} }$$

and driven by an external fluctuating signal δ ik D t ( ) only present at node k . Here θ i and P i denote the rotation angle and the natural acceleration of oscillator i (proportional to power input or output at i ), α &gt; 0 parametrizes the damping coefficient and K ij &gt; 0 the coupling strength of node pair ( i, j ). The model is equivalent to a coarse-grained model of AC power grids [12] enabling effective inertia for units with fluctuating power input from renewables. It describes the collective dynamics of N sub-grids with an effective damping coefficient α rotating at grid frequency Ω 0 = 2 π × 50 Hz (or 60 Hz in the US and parts of Japan) in the normal operation state. In this context θ i represents the center-of-inertia angle deviation of sub-grid i to the reference frame rotating at Ω 0 , P i the power generated ( P + &gt; 0) or consumed ( P -&lt; 0) in sub-grid i and K ij the line capacity of the power transmission between sub-grid i and j . The driving signal δ ik D t ( ) is a fluctuating time series additive to the average power generation or consumption at sub-grid k . For the purpose of the modelling setup each sub-grid is presented as one node of a graph.

How does such a network respond to dynamic input signals (Fig. 1)? Close to a normal operation state of the power grid, i.e. a stable fixed point θ ∗ := ( θ ∗ 1 , · · · , θ ∗ N ) of the oscillatory network, the collective response Θ ( ) t := θ ( ) t -θ ∗ to a perturbation vector D ( k ) ( ) defined via its t components D ( k ) i ( ) := t δ ik D t ( ) is accurately given by a linear response theory [12]

$$\ddot { \Theta } ^ { ( k ) } = - \alpha \dot { \Theta } ^ { ( k ) } - \mathcal { L } \Theta ^ { ( k ) } + \mathbf D ^ { ( k ) }, \quad \ ( 2 ) \quad ^ { \ t u t i } _ { \mathcal { + } \sigma \dots }$$

glyph[negationslash]

where L with L ij := K ij cos( θ ∗ j -θ ∗ i ) for i = j and L ii = -∑ j = i L ij is a weighted graph Laplacian. The linear glyph[negationslash]

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

FIG. 1. Complex network resonance patterns and large amplitude responses. (a 1 -a 3 ) Distinct network resonance patterns for three driving frequencies. Each node is color-coded by the relative response strength A ∗ ( k ) i := | ˙ Θ ( k ) i | / lim ω → 0 | ˙ Θ ( k ) i | and the driven node is marked with a black square. In simulation we use per-unit parameters α = Ω 0 ¯ D/ H 2 ¯ , P + = Ω 0 / H 2 ¯ , P -= -P + / 3, K ij = 2 P + with the normalized damping coefficient ¯ D = 0 02 s . 2 and the aggregated inertia constant ¯ H = 4s. The network topology is created by a random growth model of power grid networks [31].(b) The resonant response patterns varies with driving frequency in the resonance regime (shaded orange). The relative response strength of individual nodes are color-coded by its distance to the driven unit. Grey vertical lines indicate the N -1 resonance frequencies ω [ glyph[lscript] ] res .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00957v1-figure-18754666-47f1-43d4-aeae-fb8d900d14db.jpg)

/

network response is analytically solvable by projecting it to the orthogonal eigenspaces of the Laplacian matrix [12]. For a sinusoidal perturbation D t ( ) = εe i ( ωt + ) ϕ , the frequency response at unit i to a perturbation at unit k reads

$$\dot { \Theta } _ { i } ^ { ( k ) } ( \omega, t ) = \varepsilon e ^ { \imath ( \omega t + \varphi ) } \sum _ { \ell = 0 } ^ { N - 1 } \frac { i \omega v _ { k } ^ { [ \ell ] } v _ { i } ^ { [ \ell ] } } { - \omega ^ { 2 } + \imath \alpha \omega + \lambda ^ { [ \ell ] } }. \quad ( 3 )$$

Here λ [ glyph[lscript] ] and v [ glyph[lscript] ] i denote, respectively, the glyph[lscript] -th eigenvalue and the i -th component of the corresponding eigenvector. The eigenvalues are indexed as 0 = λ [0] ≤ · · · ≤ λ [ N -1] . When one of the eigenmodes is excited, that is, when the perturbation frequency ω maximizes the contribution of an eigenmode in Eq. 3 with an eigenfrequency

$$\omega = \omega _ { \text{res} } ^ { [ \ell ] } \coloneqq \sqrt { \lambda ^ { [ \ell ] } - \frac { \alpha ^ { 2 } } { 4 } }, \text{ \quad \quad } ( 4 )$$

the network response is dominated by the resonant eigenmode characterized by an overlap factor v [ glyph[lscript] ] k v [ glyph[lscript] ] i , constituting a nontrivial, highly heterogeneous dynamic pattern (Fig. 1).

We emphasize that network resonances emerge not only at a single perturbation frequency, but at all fre-

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

quencies in a wide frequency range we call the resonance regime, i.e.

$$\omega \in I _ { \text{res} } \coloneqq \left [ \sqrt { \lambda ^ { [ 1 ] } - \frac { \alpha ^ { 2 } } { 4 } }, \sqrt { \lambda ^ { [ N - 1 ] } - \frac { \alpha ^ { 2 } } { 4 } } \right ]. \quad ( 5 ) \quad \widehat { \gamma } _ { \text{i.e.} }$$

Extraordinarily high strengths of frequency response up to an order-of-magnitude (e.g. 12 times) larger than the homogeneous response strength in the low frequency limit [12] may appear across the network (Fig. 1d). Furthermore, the resonant response pattern sharply depends on the driving frequency. In an exemplary network, the response pattern appears to be distinctly different for three frequencies with no more than 1Hz apart from each other (Fig. 1a-c). Beside the heterogeneity in response amplitude, each node's response additionally exhibits heterogeneous phase delay towards the perturbation signal, due to the characteristic arguments of the complex responses (3).

## III. INDEXING RESONANT RESPONSES

Even perturbed only by a single-frequency resonant signal, the network already exhibits complex response patterns in terms of the strength and the phase delay of the sinusoidal response (3). In reality, power grids are constantly exposed to noisy fluctuations in renewable power generation and in the power consumption of households and industry, which consist of Fourier components with a wide range of frequencies in the resonance regime, stochastic magnitudes and random phases. For a given noisy perturbation time series, the network response time series is computable by summing up the linear response to each frequency (3). However, the influence of future remains unknown. Making general predictions for a network's complex resonant response to noisy fluctuations, and particularly, identifying the most susceptible nodes still remains a challenge.

We propose an index of the vulnerability of individual nodes in a network under resonant perturbations, the Dynamic Vulnerability Index (DVI), which helps to rank the maximum resonant response magnitude for perturbations with a characteristic PSD S ω ( ). Measurements of wind and solar power systems indicate that both of the strongly fluctuating renewable power sources are characterized by a power-law PSD with the Kolmogorov exponent -5 / 3, see Ref. 32. Such a characteristic PSD allows for estimating Fourier components' amplitudes through the relation ε ω ( ) ∝ S ω ( ) 1 2 . We thus define the DVI for node i given a noisy perturbation with PSD S ω ( ) driving node k as

$$D V I _ { i } ^ { ( k ) } = \sum _ { \omega \in I _ { r e s } } S ( \omega ) ^ { \frac { 1 } { 2 } } \left | \sum _ { \ell = 0 } ^ { N - 1 } \frac { \omega v _ { k } ^ { [ \ell ] } v _ { i } ^ { [ \ell ] } } { - \omega ^ { 2 } + \imath \alpha \omega + \lambda ^ { [ \ell ] } } \right |. \quad ( 6 ) \quad \text{matic} }$$

Here we assume that the response to a noisy fluctuation, as the real part of the sum of the complex responses

˙ Θ ( k ) i ( ω ) (3) to each Fourier component of frequency ω , approaches the sum of the magnitudes of Θ ˙ ( k ) i ( ω ) for sufficiently long time series with length T .

)

(

FIG. 2. The ranking of DVI well predicts the ranking of maximum resonant responses. A time series of a resonant driving signal (a) is extracted from a colored noise with a PSD exponent -5 / 3 (b). According to the PSD of the original noise (black dots), the resonant signal is obtained by filtering the frequencies in the original signal and keeping only the ones (orange dots) in the resonance regime (shaded in orange). The ranking of the maximum frequency response max t ∈ [0 ,T ] ∣ ∣ ∣ ˙ Θ ( k ) i ( ) t ∣ ∣ ∣ in the simulation T = 100 s is shown by a heat map in (c) and the ranking of the DVI (6) with S ω ( ) ∝ ω -5 / 3 also by a heat map in (d). The perturbed node is marked with a black square. Three exemplary unexpected vulnerable nodes (marked with circles) are successfully predicted by the DVI. The network settings is the same as in Fig. 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00957v1-figure-088a1471-3944-4565-97fb-dbfab155662d.jpg)

Note that the PSD S ω ( ) gives only the scaling information of ε ω ( ), so the relative value of DVI among all nodes within a network is more relevant rather than its absolute value. The ranking of DVI thus provides information about which nodes are most susceptible rather than predicting the actual response magnitudes. For instance, in a 100-second simulation (Fig. 2), the ranking of the numerically determined maximum frequency response appears to be highly similar to the ranking given

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

(

n

u

)

&gt;

/

ranking of NRI

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

Z

/

I

J

D

r

X

Y

0

P

C

z

V

w

A

y

8

&gt;

E

n

c

j

N

d

9

M

F

W

2

q

H

L

S

U

m

o

u

k

p

5

R

7

B

+

G

K

O

3

f

v

Q

T

g

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

z

0

E

T

R

F

u

m

Z

8

P

+

V

2

g

w

&gt;

A

n

c

j

N

d

9

M

W

K

I

H

X

L

S

O

Q

f

D

C

J

v

U

3

y

Y

k

/

5

B

G

o

r

7

q

p

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

S

z

C

+

7

D

g

U

p

5

R

n

O

/

V

M

K

&gt;

A

E

c

j

N

d

9

w

F

P

W

2

q

8

I

H

X

L

m

Y

0

o

u

k

B

r

3

f

J

Z

y

Q

G

v

T

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

+

m

W

B

M

v

H

N

y

J

7

5

R

j

Z

O

2

z

p

c

Y

&gt;

A

E

g

n

V

d

9

w

F

P

K

8

I

X

L

S

0

Q

f

D

C

U

3

k

/

u

o

T

G

r

q

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

r

K

T

3

o

2

V

M

q

0

C

7

L

O

B

&gt;

A

E

c

j

N

d

9

w

F

P

W

v

5

S

k

X

I

u

+

8

p

Q

/

G

J

R

f

n

z

Z

H

g

D

U

m

Y

y

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

5

F

w

d

8

r

2

m

z

3

U

c

S

W

9

v

A

M

&gt;

E

n

j

V

N

X

u

y

Z

7

I

D

+

K

0

J

k

O

B

R

o

f

T

H

q

p

P

g

Y

G

C

L

/

Q

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

u

q

/

g

+

m

I

7

S

5

F

C

0

c

U

p

v

2

o

N

D

k

&gt;

A

E

n

j

V

d

9

M

w

P

W

K

8

H

X

y

Q

R

f

r

Z

O

T

L

Y

B

3

z

J

G

&lt;

l

a

t

e

x

i

s

h

1

\_

b

6

4

=

"

r

K

T

3

o

2

V

M

q

0

C

7

L

O

B

&gt;

A

E

c

j

N

d

9

w

F

P

W

v

5

S

k

X

I

u

+

8

p

Q

/

G

J

R

f

n

z

Z

H

g

D

U

m

Y

y

FIG. 3. Fast-converging, robustly high prediction performance of DVI ranking. (a) The prediction error drops sharply at first, decreases over time and rests at a level about 85% lower than the error of random guesses. The Gaussian error distribution of random guesses [33] is indicated by various shades of grey: black line for the expectation value of E for random guesses, areas colored with different shades of grey for intervals 1 ± ρ, 1 ± 2 ρ, 1 ± 3 ρ respectively with ρ being the ratio between the standard deviation and the expectation value of E for random guesses. (b) A snapshot of the prediction error's dependence on signal's PSD exponent at T = 100 . The s prediction error increases slightly for larger PSD exponents of the driving signal. In both panels the prediction error is color-coded by the PSD exponent of the driving signal. (ce) At T = 100 , the DVI exhibits a high correlation s with the maximal frequency response for various PSD exponents b ∈ { 0 1 2 , , } .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00957v1-figure-fa5c3fca-5c13-40a1-8e09-5b18f2a1a757.jpg)

by DVI. Particularly, it gives warnings about some particular nodes (marked in circles) which are far away from the perturbation, unexpected based solely on topology (e.g. not dead-ends [34]), but especially vulnerable due to resonances.

We further quantitatively investigate the prediction performance of DVI in terms of its robustness over time and over the stochastic features of fluctuations. We measure DVI's prediction error with a normalized Spearman's footrule distance [33] between the predicted ranking of the maximal ˆ σ DVI ( ) and the actual ranking i σ i ( )

$$E \coloneqq \frac { 1 } { E _ { \text{rand} } } \sum _ { i = 1 } ^ { N } | \hat { \sigma } _ { \text{DIV} } ( i ) - \sigma ( i ) | \,. \quad \quad ( 7 ) \quad \text{cab} | \, \underset { F } { \text{of} }$$

The prediction error is normalized by the expectation value of the Spearman's footrule distance E rand = N / 2 3 between two random rankings chosen independently and uniformly in the set S N of permutation of N elements, see Ref. 33. Numerical results show that the ranking σ of the maximum response from direct simulation converges fast to the a priori DVI ranking ˆ σ DVI (Fig. 3a). For a 100 second perturbation time series, we measure the true ranking σ every 0 1 . second and compute the footrule distance E . For a sample power grid with the eigenfrequencies ∼ 1 Hz ( I res = [0 32 Hz 3 74 Hz]), . , . the prediction error drops about 80% in the first 10 seconds and continues to decrease slowly. At T = 100s, the prediction error drops to about 15% of random guess error level and the DVI is highly correlated to the maximal frequency responses with Pearson's correlation coefficient r larger than 0 985 (Fig. 3c-e). . Furthermore, we find the prediction performance of DVI is quite robust over time and for different types of colored noise with the powerlaw exponent b ∈ [0 , 2], from white noise to brown noise. The prediction error remains at almost the same level and shows only a mild increase with growing b (Fig. 3b).

## IV. CONCLUSION

We presented a measure of dynamic node vulnerability to predict the most resonant nodes in stochastically driven oscillator networks and specifically AC power grid models. Based on a linear response theory of the network's resonance patterns for a single frequency, we propose to estimate the susceptibility of a node to stochastic driving signals by i) estimating the driving signal's Fourier spectrum by its PSD characteristics and ii) accumulating the response amplitudes to each Fourier components. Numerical results indicate strong prediction power of the proposed DVI in identifying the most resonant nodes. The true ranking of maximum response from direct simulation converges fast to the ranking prediction given by DVI, thus revealing the most vulnerable nodes. The prediction performance is robust not only over time, but also for various types of colored noise sources. For all tested settings, ranging from white noise to brown noise b ∈ [0 , 2], the prediction error stays at a low level and the DVI ranking highly correlates with the true response ranking, with a Pearson's correlation coefficient larger than 0 985. .

Given the position and the characteristic PSD of the driving signal, for instance the location of a wind farm in a power grid network, the ranking of the DVI obtained from (6) may help to identify which stations in the power grid would particularly be influenced by the resonant signals carried by the fluctuating wind power input, allowing precautionary measures to be taken. As the method is robust and computationally fast, it might also be applicable ad hoc if the network changes after failures or other unforeseen events such as load shedding.

Furthermore, the DVI may support optimizing future power grid planning. For instance, new stations or new lines should be built in a way that important units in the network would not suffer from severe resonant disturbances in the altered network topology.

Taken together, the proposed dynamic vulnerability index provides a powerful tool to rank the nodal resonance level in networks of dynamical units. Whereas the presentation above is focused on the second order Kuramoto model to overcome previous analytic limitations, the index is readily generalized, both to phase oscillator networks (with single variable nodes) and more complex systems such as networks of the third order model of power grids and more generally, networks of oscillatory and non-oscillatory dynamic units.

- [1] F. Hussain, C. Gupta, A. J. Hirning, W. Ott, K. S. Matthews, K. Josic, and M. R. Bennett, Proc. Nat. Acad. Sci. 111 , 972 (2014).
- [2] S. Jahnke, R.-M. Memmesheimer, and M. Timme, PLoS Comput. Biol. 10 , e1003940 (2014).
- [3] P. Kundur, N. J. Balu, and M. G. Lauby, Power System Stability and Control, ser. The EPRI Power System Engineering Series , Vol. 7 (McGraw-Hill New York, 1994).
- [4] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, Eur. Phys. J. B 61 , 485 (2008).
- [5] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Phys. Rev. Lett. 109 , 064101 (2012).
- [6] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Nature Physics 9 , 191 (2013).
- [7] F. D¨ orfler, M. Chertkov, and F. Bullo, Proceedings of the National Academy of Sciences 110 , 2005 (2013).
- [8] T. Coletta and P. Jacquod, Physical Review E 93 , 032222 (2016).
- [9] D. H. Zanette, Europhys. Lett. 68 , 356 (2004).
- [10] D. H. Zanette, Eur. Phys. J. B 43 , 97 (2005).
- [11] M. Tyloo, T. Coletta, and P. Jacquod, Phys. Rev. Lett. 120 , 084101 (2018).
- [12] X. Zhang, S. Hallerberg, M. Matthiae, D. Witthaut, and M. Timme, Science Advances 5 , eaav1027 (2019).
- [13] H. Haehne, K. Schmietendorf, S. Tamrakar, J. Peinke, and S. Kettemann, Physical Review E 99 , 050301 (2019).
- [14] B. Sch¨ afer, M. Matthiae, X. Zhang, M. Rohden, M. Timme, and D. Witthaut, Phys. Rev. E 95 , 060203 (2017).
- [15] S. Kettemann, Phys. Rev. E 94 , 062311 (2016).
- [16] S. Tamrakar, M. Conrath, and S. Kettemann, Sci. Rep. 8 (2018).
- [17] D. Manik, M. Rohden, H. Ronellenfitsch, X. Zhang, S. Hallerberg, D. Witthaut, and M. Timme, Physical Review E 95 , 012319 (2017).
- [18] B. Sch¨ afer, C. Beck, K. Aihara, D. Witthaut, and M. Timme, Nature Energy 3 , 119 (2018).

## ACKNOWLEDGMENTS

We gratefully acknowledge support from the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under Germanys Excellence Strategy EXC-2068 390729961 Cluster of Excellence Physics of Life and the Cluster of Excellence Center for Advancing Electronics at TU Dresden, the German Federal Ministry for Research and Education (BMBF grants no. 03SF0472F and 03EK3055F), and the Helmholtz Association (via the joint initiative Energy System 2050A Contribution of the Research Field Energy and grant no. VH-NG-1025).

- [19] M. Tyloo, R. Delabays, and P. Jacquod, Phys. Rev. E 99 , 062213 (2019).
- [20] J. Hindes, P. Jacquod, and I. B. Schwartz, arXiv:1904.12174 (2019).
- [21] M. Tyloo, L. Pagnier, and P. Jacquod, arXiv:1810.09694 (2018).
- [22] S. Auer, F. Hellmann, M. Krause, and J. Kurths, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 , 127003 (2017).
- [23] A. Plietzsch, S. Auer, J. Kurths, and F. Hellmann, arXiv:1903.09585 (2019).
- [24] S. Arianos, E. Bompard, A. Carbone, and F. Xue, Chaos: An Interdisciplinary Journal of Nonlinear Science 19 , 013119 (2009).
- [25] I. Simonsen, L. Buzna, K. Peters, S. Bornholdt, and D. Helbing, Physical Review Letters 100 , 218701 (2008).
- [26] R. Albert, I. Albert, and G. L. Nakarado, Physical review E 69 , 025103 (2004).
- [27] H. Song and M. Kezunovic, in 2005 IEEE/PES Transmission &amp; Distribution Conference &amp; Exposition: Asia and Pacific (IEEE, 2005) pp. 1-7.
- [28] G. Qu, J. Rudraraju, R. Modukuri, S. Hariri, and C. S. Raghavendra, in IASTED International Conference on Communications, Internet, and Information Technology, November 18-20, 2002, St. Thomas, US Virgin Islands, USA (2002) pp. 289-294.
- [29] R. Guti´ errez, F. Del-Pozo, and S. Boccaletti, PloS one 6 , e20236 (2011).
- [30] M. Tyloo and P. Jacquod, arXiv:1905.03582 (2019).
- [31] P. Schultz, J. Heitzig, and J. Kurths, Eur. Phys. J. Spec. Top. 223 , 2593 (2014).
- [32] M. Anvari, G. Lohmann, M. W¨ achter, P. Milan, E. Lorenz, D. Heinemann, M. R. R. Tabar, and J. Peinke, New J. Phys. 18 , 063027 (2016).
- [33] P. Diaconis and R. L. Graham, J. Roy. Stat. Soc. Ser. B , 262 (1977).
- [34] P. J. Menck, J. Heitzig, J. Kurths, and H. J. Schellnhuber, Nature Comm. 5 , 3969 (2014).
Xiaozhu Zhang Cheng Ma Marc Timme
Dynamic Vulnerability in Oscillatory Networks and Power Grids
Recent work found distributed resonances in driven oscillator networks and AC
power grids. The emerging dynamic resonance patterns are highly heterogeneous
and nontrivial, depending jointly on the driving frequency, the interaction
topology of the network and the node or nodes driven. Identifying which nodes
are most susceptible to dynamic driving and may thus make the system as a whole
vulnerable to external input signals, however, remains a challenge. Here we
propose an easy-to-compute Dynamic Vulnerability Index (DVI) for identifying
those nodes that exhibit largest amplitude responses to dynamic driving signals
with given power spectra and thus are most vulnerable. The DVI is based on
linear response theory, as such generic, and enables robust predictions. It
thus shows potential for a wide range of applications across dynamically driven
networks, for instance for identifying the vulnerable nodes in power grids
driven by fluctuating inputs from renewable energy sources and fluctuating
power output to households.