Linear Response Theory for Renewable Fluctuations in Power Grids with Transmission Losses

We study the spreading of renewable power fluctuations through grids with Ohmic losses on the lines. By formulating a network adapted linear response theory, we find that vulnerability patterns are linked to the left Laplacian eigenvectors of the overdamped eigenmodes. We show that for tree-like networks fluctuations are amplified in the opposite direction of the power flow. This novel mechanism explains vulnerability patterns that were observed in previous numerical simulations of renewable micro-grids. While exact mathematical derivations are only possible for tree like networks with homogeneous response, we show that the mechanisms discovered also explain vulnerability patterns in realistic heterogeneous meshed grids by studying the IEEE RTS-1996 test system.


I. INTRODUCTION
A fundamental challenge for the operation and control of power grids is to maintain the balance between power production and power demand.In AC power systems that are dominated by conventional generators, the power fluctuations on the demand side are balanced by the control schemes of the production side to maintain a stable frequency at 50 or 60 Hz, respectively.However, with the ongoing integration of highly intermittent renewable energy sources such as wind and solar, there are not only fluctuations on the demand side but also on the generation side. 1 Demand fluctuations are typically uncorrelated and can, therefore, average out for a large number of consumers.In contrast, the power fluctuations in large wind and solar farms stem from the same meteorological conditions and can, therefore, be highly correlated.As a result, these fluctuations add up and can lead to large fluctuations of power production at single nodes in the network system. 2 The impact of noise on the stability of the synchronous state in complex dynamical systems has been intensively studied with the method of linear response theory.6][7] The spreading of intermittency from fluctuations to the frequency response throughout a lossless network was calculated by Haehne et al. 8 Zhang et al. 29 identified three frequency regimes of the network response networks: a bulk, a resonant, and a local regime.The bulk regime covers low-frequency perturbations, and the network responds as a whole.In contrast, as already pointed out by Kettemann et al., 9 high-frequency perturbations stay localized at the fluctuating node and decay exponentially.They constitute the

ARTICLE
scitation.org/journal/chalocal regime.The resonant regime is where the fluctuation spectrum and the oscillatory dynamics of the network overlap and produces complex resonant response patterns.
To our knowledge, all prior analytic works on linear response in power systems consider rather simple power system models.In this work, we want to transfer this theoretical knowledge and develop a linear response theory that is well suited to also describe more realistic power systems, including higher order dynamical models of inverters or generators.Of particular importance is that our approach is capable of dealing with transport losses on the lines.Mathematically such systems can be represented by an asymmetric effective network Laplacian.We derive an approximation of the response that is highly predictive of the actual behavior of the system in many cases of interest.The theory is used to explain key features of the complex phenomenology that was numerically observed for renewable fluctuations in an AC microgrid model. 10The major finding is that auto-correlated power fluctuations are enhanced in the opposite direction of the power flow due to Ohmic losses on the lines.With simulations in the IEEE RTS-1996 11 test case, we are able to show that this mechanism is relevant also for more realistic systems with a meshed topology and heterogeneous parameters.The fact that renewable fluctuations in load-heavy regions of the grid have a higher impact on the frequency stability is of high relevance, for example, for the connection of new wind parks to a grid.

II. POWER FLOW NETWORKS WITH LOSSES
We first give the general form of the power grid equations we want to consider and derive some general properties of the linearization that help characterize the system.
The network structure of the grid can be represented by a graph G = (N , E), with a set of N nodes corresponding to generators and loads and a set of E edges corresponding to the transmission lines that carry the power flow.The dynamical state of each node i is represented by x i (t) : R → R D i .In the following, we will use the notation that node indices are denoted by superscripts and variable indices are denoted by subscripts, such that x i l is the lth variable of the ith node.The state of the entire dynamical system x(t) : R → R S contains the states of all components, with a total system size S = N i=1 D i .We assume that every node i is coupled with its adjacent nodes by a power flow P i (x) : R S → R that depends only on the state difference ẋi = f(x i , P i (x)), Here, P ij (•) is the signed power flowing on the line ij as a function of the node states, and as seen from node i.If no power is lost on the line, we have P ij = −P ji .In Eq. (1), we made two additional assumptions, 1.The node dynamics is homogeneous, i.e., f i = f j = : f. 2. The power flow depends only on one internal state variable x θ , e.g., the voltage phase angle for AC power grids and the absolute voltage for DC power grids.
Furthermore, we require some natural properties to hold for the power flow.
• Positive power flow increases with state difference: • Losses increase with increasing power flow: From this, it follows that if P ij > 0, that is, we have power flowing from i to j, we have We assume the system has a stable stationary state x = ξ .The change of the power flow at node i for a deviation δx from the stationary state is given by We define a matrix L ij := ∂P i ∂x j θ (ξ θ ).Inserting the power flow equation, we see that this has the form of a weighted Laplacian matrix with weights Usually, Laplacians are defined to be symmetric matrices, with the underlying assumption being the conservation of flow on the links.However, if we consider transport losses, the Laplacian matrix describing the diffusion dynamics on the linear level is asymmetric, i.e., w ij = w ji .
Similar to the symmetric case, asymmetric Laplacians always have an eigenvalue λ 1 = 0.The corresponding right eigenvector is homogeneous v (1) r,i = v (1) r,j for all j.The corresponding left eigenvector, however, is generally heterogeneous.It is determined by the equation In tree-like networks, this gives a strict relation for the eigenvector entries of two neighboring nodes.We can see this by starting at the nodes with degree one.At these nodes, there is only one summand F ij , which, therefore, has to be zero.All the summands are by definition antisymmetric F ij = −F ji and, therefore, we know that the corresponding summand F ji in the condition for the neighboring node is also zero.By going up the tree structure to nodes of a higher degree we can eliminate the summands corresponding to all previously visited nodes.Doing this, we see that in the above equations not only the sum is equal to zero but every single summand has to Chaos ARTICLE scitation.org/journal/chabe zero itself and, therefore, v (1) l,i v (1) l,j = w ji w ij .
Assume the power is flowing from node i to node j.From Eq. ( 2), it follows that w ij > w ji and, hence, in a tree network the entries of the left eigenvector corresponding to λ 1 = 0 of the asymmetric Laplacian are increasing along the power flow in the network Since we assume that the nodes are homogeneous, we can factorize the linearization of the dynamical system (1) into a network part and a local part.In Appendix A it is shown that the Jacobian of this system can be written in the form with matrices A, B ∈ R D×D , the weighted Laplacian L ∈ R N×N as defined in Eq. (3) and ⊗ denoting the Kronecker product.The eigenvectors and eigenvalues of this Jacobian take the form where λ, v are the eigenvalues and eigenvectors of the Laplacian and µ(λ), u(λ) are the eigenvalues and eigenvectors of the matrix C(λ) = A + λB.In the following, we will use n and m to denote the multi-index (a, b).

III. LINEAR RESPONSE THEORY
In the following, we want to calculate the response of (1) to an additive fluctuation η(t) : R → R S .We assume the system has a stable fixed point x(t) = ξ .The linear response of the deviation δx(t) : = x(t) − ξ is then given by where χ (t) = θ (t)e Jt is the response function of the system.In the Fourier space, the convolution reduces to a simple product We quantify the response signal by applying the L 2 norm.For a single system variable, it is defined as In the following, we will restrict ourselves to the case of single node fluctuations, where there is only one non-zero entry η j (t) driving the system.In principle, it is straightforward to generalize our approach to fluctuations at multiple nodes.In that case, not only the autocorrelation but also the cross-correlation of fluctuations has to be taken into account.However, the focus on single node fluctuations is sufficient to understand the effects of auto-correlated fluctuations and transport losses.For single node fluctuations, the L 2 norm of the response is given by where S η j η j (ν) = | ηj (ν)| 2 is the power spectrum of the fluctuation.
The response matrix can be decomposed into the response of the single eigenmodes χ (ν) = n χ (n) (ν).In Appendix B, it is shown that these mode response functions are given by where σ n are the eigenvalues of the Jacobian and v (n) r v (n) l is the outer product of the corresponding right and left eigenvectors.Inserting the mode expansion into Eq.( 7) yields a sum of single mode terms the single mode terms are given by where L (n) are Lorentzian functions with width γ n and maximum at ν = ν n , For auto-correlated perturbation signals, the integral ( 7) is generally hard to solve, particularly, when the power spectral density is not known analytically.In that case, we could only determine the power spectral density from measurements and compute the L 2norm semi-analytically.An analytical approximation for the single mode terms can be calculated for small γ , the low damping regime.In the limit γ n → 0, the Lorentzian function converges toward a Dirac delta distribution and, hence, for small γ n we can approximate the integral as This approximation is valid if the spectral density does not vary much over width of the Lorentzian functions S η j η j (ν) ≈ S η j η j (ν + γ n ).In Appendix C, it is further shown that for small damping parameters γ n the cross-mode terms are suppressed.Neglecting these terms is valid if the mode dampings are much smaller than their spectral distance γ n , γ m |ν n − ν m |.The L 2 norm of the response can then be approximated as which we call the peak approximation.We see that the response is given by a superposition of the different mode contributions.How strongly a certain mode is excited depends on the power spectral density at the eigenfrequency ν n and the entry of the left corresponding eigenvector at the perturbed node, whereas the response Chaos ARTICLE scitation.org/journal/chastrength at different nodes is determined by the entries of the right corresponding eigenvector.

IV. EXAMPLE: AC MICROGRID MODEL
We now will turn toward an example system to demonstrate the power of our theoretical approach.We analyze the impact of turbulent fluctuations on power grids by characterizing the vulnerability and excitability of the different nodes in the network.Here, we define the vulnerability of a node as the strength of the network response to power fluctuations at this particular node.Conversely, the excitability of a node is the strength of the response at this node given power fluctuations at another node.Using our approximation of the L 2 norm, we are able to explain three previously observed properties of renewable power grids. 10There is a pronounced fine structure in both vulnerability and excitability of nodes.2. Losses on the lines lead to a pronounced network structure in the vulnerability of the nodes, but not in which nodes are excited.3. The vulnerability appears high in parts of the network that are consumer-heavy, and within these areas tends to rise the further away from the center of the network the node is.
In the following, we will provide an analytical explanation for these observations for an islanded microgrid with fluctuating renewable in-feed.Following Schiffer et al., 12 droop-controlled inverters with virtual inertia in their simplest form can be modeled by the swing equation We include resistive losses of the power flow on the lines via the conductance matrix G ik Furthermore, we assume that the power in-feed at each node is composed of a constant and a small fluctuating part We simulate this system on a 100-node network generated by a random growth model for power grids. 13The power fluctuation signal is generated by a combination of stochastic wind and solar power fluctuation models. 14,15The power spectrum of the resulting signal is power-lawed with the Kolmogorov exponent of turbulence.Further details on the parametrization and the fluctuation modeling can be found in Appendix D. The response at each node is quantified in terms of the L 2 norm of the frequency deviation from the nominal grid frequency ω s The response of the entire dynamical system S can then be quantified in terms of the L 2 norm of the average deviation from the nominal grid frequency It should be noted that this measure is different from the synchronization norm that has been used in most of the studies on the response of swing equations in the linear regime 5,7,16,17 While being useful to study the synchronicity in the network, this measure by definition omits any fluctuation of the bulk of synchronous frequencies.However, as will be shown in the following, this bulk behavior turns out to be the most dominant mode in the frequency response to renewable power fluctuations.Furthermore, due to the presence of losses, this mode is no longer homogeneous throughout the network.As we will see, it can completely dominate the effect of network structure on the systems node-wise vulnerability.

A. Fine structure of network responses
Following Auer et al., 10 we simulate single node fluctuations in the full nonlinear system (9) for every pair of perturbed (input) and observed (output) nodes and depict the response strength as a color-coded matrix plot (Fig. 1).Here, the horizontal and vertical lines correspond to nodes with large vulnerability or excitability, respectively.Comparing the simulation of the linearized system with the full nonlinear system shows that the main response pattern also remains in the linearized dynamics.However, some nonlinear artifacts are not captured by the linearized model.These nonlinear effects cannot be analyzed with our linear theory.For a given time series of a power fluctuation δP(t), we can numerically determine its power spectral density and thereby semi-analytically compute the peak approximation (8) for the L 2 norm of the frequency responses.In Fig. 1, it can be seen that this approximation can reproduce the response pattern of the linearized system.This means that by only knowing the eigenvectors of the system and the spectral density of the power fluctuations at the eigenfrequencies of the system, we can analytically predict the response strength at every node in the network.The full network response is a superposition of the responses of every single mode.The contribution of each eigenmode is determined by the spectral excitation factor and the response pattern by the left and right eigenvectors of that mode.When a single mode is dominating the response of the network, the vulnerability and excitability of the nodes can be linked to the left and right eigenvectors of this mode.

B. Line losses and the bulk mode
In the following, we will assume homogeneous damping and inertia parameters D i = D, M i = M.The Jacobian of the system (9) can then be written in the form ( 5), with matrices

Chaos
and the Laplacian weights From Eq. ( 6), it follows that For the Laplacian eigenvalue λ 1 = 0, we have two Jacobian eigenvalues, σ (+,1) = 0 and σ (−,1) = − D M .The eigenvalue σ (+,1) corresponds to the symmetry of homogeneous phase shifts that do not contribute to the dynamics, whereas σ (+,1) corresponds to homogeneous frequency shifts leading to an exponentially decaying response of the nodes' frequencies with rate D M .When the algebraic connectivity of the network, i.e., the second smallest eigenvalue of the Laplacian, fulfills the condition λ 2 > D 2 4M , the square root term in (13) will be imaginary and, therefore, σ (−,1) is the only overdamped mode in the system.In this case, the mode fully determines the behavior of the system in the bulk regime and we, therefore, refer to it as the bulk mode.
When the algebraic connectivity is significantly larger than the threshold λ 2 D 2 4M , the eigenfrequencies of all the other system modes are rather high.For correlated fluctuations, the power spectrum at high frequencies is suppressed 2,18 and, therefore, we find that the network response in this regime is entirely dominated by the bulk mode.The right Laplacian eigenvector of this mode is homogeneous, whereas the left eigenvector has heterogeneous entries.This means that all nodes are equally excited but certain nodes have much higher vulnerability to power fluctuations.The resulting dynamical asymmetry corresponds to the continuous horizontal lines in the bulk mode plot (Fig. 1) and has its origin in the Ohmic losses of the lines.

C. Vulnerability in tree-networks
When the bulk mode is dominating the response of the system, the difference of the vulnerability of nodes to power fluctuations is entirely determined by the left eigenvector of this mode.From Eq. ( 4) it follows that entries of this eigenvector are increasing along the power flow.This means that fluctuations are amplified in the opposite direction of the power flow and the nodes located at the sinks of the flow are the most vulnerable.In Fig. 2, it can clearly be seen that the network branch where the power is flowing from the center toward the outlying nodes is much more vulnerable than the network branches where the power is flowing toward the center.This explains the observation by Auer et al. 10 that the vulnerability of nodes and the closeness centrality of the network are closely related.

V. EXAMPLE: IEEE RELIABILITY TEST SYSTEM-1996
As the last example, we simulate the IEEE Reliability Test System-1996 11 to show that our findings are still valid in a much more realistic test case, including a meshed grid topology, (algebraically modeled) load buses, and heterogeneous generator parameters.Similar to the example in Sec.IV, we simulate single node power fluctuations using a combination of stochastic wind and solar power fluctuation models. 14,15Imagine a large solar or wind farm connected to one of the buses in the system.Since the fluctuations at single solar panels and wind turbines can be highly correlated within a farm, these fluctuations will not average out but rather add up to a large power fluctuation at the respective node.We simulate single node power fluctuations at the load buses and measure the frequency response at the generator buses in the system using the deviation norm (12).The results are depicted in Fig. 3.It can be seen that the vulnerability toward power fluctuations at a bus increases in the direction of the power flow.Accordingly, the most vulnerable buses are located in the parts of the grid with a very high share of loads.This indicates that our findings concerning the joined effect of auto-correlated power fluctuations and line losses are also valid in a more realistic system setup.

VI. CONCLUSION
In this paper, we presented a linear response theory for correlated fluctuations in power grids that was used to derive an approximation for the frequency response that turned out to be highly accurate in the regime studied.We have shown that important features of the response can be understood as a consequence of the Laplacian that describes the dynamical coupling of nodes in the network being asymmetric in the presence of Ohmic losses on the power lines.In particular, we were able to fully explain the structures in the node vulnerability that have been observed in numerical simulations of an islanded microgrid with high renewable penetration. 10or tree-like grids, the location of vulnerable nodes is related to the power flow throughout the network.In particular, fluctuations at net power flow sinks result in a strong frequency response at all nodes in the network.For tree-like grids with very unbalanced power production, we find that the consumer heavy-branches are much more vulnerable to turbulent in-feed of renewable power.This effect is a direct consequence of both the losses on the power lines and the correlated nature of renewable power fluctuations.Considering the generality of our theoretical approach, it should be mentioned that the application to power grids is not limited to the fluctuation of renewable energy sources but might also be the basis for studying the impact of demand fluctuations on grids.Also, this work was focused on single node fluctuations.However, the formulation of the response in terms of power spectra also provides an elegant starting point for understanding correlated multi-node fluctuations.In that case, not only the auto-correlation but also the cross-correlation of fluctuations will play a crucial role for understanding dynamical interactions.Furthermore, the formulation used here is suited to study active power flow variations.To understand the response of the full active and reactive power flow, as well as voltage variations, e.g., for the study of fluctuations in effective models, 19 the approach needs to be generalized to more general models than Eq.(1).Finally, in this work, we only focused on explaining the linear response patterns in the network.However, in the simulations of the full nonlinear system, we observed nonlinear resonance effects at a few nodes, which are not fully understood an need to be investigated in future research.
Taking the total derivative of the kth component of the righthand-side function at the ith node with respect to the lth variable at the jth node yields ∂p δ lθ and using the fact that the definition of the weighted Laplacian (3) yields the Jacobian (5).For a right eigenvector of this Jacobian, we have (a,b) .
The proof for the left Jacobian eigenvectors can be done similarly.Both the PDF of the fluctuations and their increment time series are fat-tailed (the tails are not exponentially bounded 25 ).The power generation from wind and solar power plants has a power spectrum that is power-lawed with the Kolmogorov exponent of turbulence. 2,18hus, these time series show long-term temporal correlations.

APPENDIX E: CIGRE MV GRID SIMULATIONS
As a second realistic example, we simulate renewable power fluctuations in the CIGRE MV grid. 26Following Ref. 27, we slightly modify the original model such that it can be operated in the islanded mode (disconnected from the higher grid level).For this, we scale the power ratings of all distributed generation units by a factor of 2 and reduce the load at node 1 such that the grid is self-supplied.Furthermore, we assume that all generation units are controllable and equipped with frequency and voltage droop control and can, therefore, be modeled as θi = ω i , The simulation results are depicted in Fig. 4. It can also be seen that for a more detailed dynamical model including voltage dynamics, the vulnerability to renewable fluctuations increases along the power flow and that power fluctuations at the node with the highest demand (node 1) cause the largest response in the network.Additionally, this example shows that closing the switches at the lines 6-7 and 4-11 can reduce the vulnerability to power fluctuations by forming meshes in the grid (see nodes 4, 5, and 6).

FIG. 1 .
FIG. 1. Color plot of the L 2 norm of the frequency response at single nodes.The color corresponds to the L 2 norm of the frequency response at the output nodes (x-axis) given a turbulent fluctuation at the input nodes (y-axis) in the grid depicted in Fig. 2. (a) Simulations of the full nonlinear system for each pair of input and output nodes.(b) Simulation of the linearized system.(c) Analytic prediction calculated with the peak approximation (8).(d) Contribution of the bulk mode to the analytic prediction.

FIG. 2 .
FIG. 2. Vulnerability of nodes in the power flow network for bulk-dominated dynamics.Nodes are colored according to the deviation norm (12) of the system response for fluctuations at this node.The link arrows indicate the direction of the power flow.

FIG. 3 .
FIG. 3. Vulnerability of load nodes to intermittent fluctuations in the IEEE Reliability Test System-1996.Load nodes are depicted as circles and colored according to the L 2 -norm of the system response for fluctuations at this node.Generator nodes are depicted as squares.The link arrows indicate the direction of the power flow.

FIG. 4 .
FIG. 4. Vulnerability of load nodes to intermittent fluctuations in the CIGRE MV grid.Nodes are colored according to the L 2 -norm of the system response for fluctuations at this node.The link arrows indicate the direction of the power flow.Lines 6-7 and 4-11 can be opened by a switch.
Langevin type model developed by Schmietendorf et al.,14 respectively.The power fluctuation δP(t) is a combination of the signals generated with these models for wind and solar power fluctuations δP(t) = 0.5δP W (t) + 0.5δP S (t).(D1)