Bounds and Estimates for the Response to Correlated Fluctuations in Asymmetric Complex Networks
BBounds and Estimates for the Responseto Correlated Fluctuations in Asymmetric Complex Networks
Anton Plietzsch,
1, 2, ∗ Sabine Auer,
1, 3
J¨urgen Kurths,
1, 2 and Frank Hellmann Potsdam-Institute for Climate Impact Research Department of Physics, Humboldt-Universit¨at zu Berlin elena international GmbH (Dated: February 26, 2021)We study the spreading of correlated fluctuations through networks with asymmetric and weightedcoupling. This can be found in many real systems such as renewable power grids. These systemshave so far only been studied numerically. By formulating a network adapted linear response theory,we derive an analytic bound for the response. For colored we find that vulnerability patterns noiseare linked to the left Laplacian eigenvectors of the overdamped modes. We show for a broad classof tree-like flow networks, that fluctuations are enhanced in the opposite direction of the flow.This novel mechanism explains vulnerability patterns that were observed in realistic simulations ofrenewable power grids. Networked systems are ubiquitous in nature and en-gineering. They include neural systems, gene expres-sion networks, and energy transport systems such as ACpower grids and gas networks. Such networks often set-tle into a steady state, in a process typically called self-organization. In all such real systems these steady statesare subject to complex external fluctuations. The inter-action of the network and the structure of fluctuationsshows a rich and barely understood interplay. For ex-ample in [1] the case of oscillator networks was studied,and a resonant regime identified in which the frequencyof the perturbation and the network structure lead to acomplex spreading behaviour through the network. Forthe case of power grids subject to renewable fluctuationscomplex response patterns with a rich phenomenologywere observed in [2].In this work we develop a linear response theorywell suited to disentangling the interplay between net-works and auto-correlated fluctuations. We derive upperbounds of the response that are highly predictive for theactual behaviour of the system in many cases of interest.This response theory is used to explain key features ofthe complex phenomenology reported in [2].Of particular importance is that the linear responsetheory given here is capable of dealing with weighted andasymmetric or directed networks, features that typicallyoccur in real systems. An example for this are energytransporting networks, where losses on links are ubiqui-tous and result in asymmetric coupling of the nodes. Ithas been shown that such line losses can fundamentallyalter both the linear [2] and non-linear behaviour of ACpower grids [3].From previous empirical work we know that both as-pects, the asymmetric coupling due to losses [2, 4] andthe autocorrelation of fluctuations [5] have to be takeninto account to get a qualitatively correct understandingof the systems dynamical behavior, and the complicated ∗ [email protected] response patterns that arise. In contrast all prior analyticworks only considered undirected, lossless networks.Besides the already mentioned work [1] dealing withharmonic perturbations, analytical results were alsogiven for singular perturbations [6], white Gaussian noise[7] and exponentially correlated noise [8–10]. The spread-ing of intermittency from fluctuations to the frequencyresponse throughout a lossless network was calculated in[11]. Of these works the latter is most similar, albeitmore cumbersome and less general to our approach.In [1] three frequency regimes of the network responsewere identified for symmetric networks: a bulk, a res-onant and a local regime. The bulk regime covers lowfrequency perturbations, and the network responds as awhole. In contrast, as already pointed out in [12], highfrequency perturbations stay localized at the fluctuatingnode and decay exponentially. They constitute the lo-cal regime. The resonant regime is where the fluctuationspectrum and the oscillatory dynamics of the networkoverlap, and produces complex resonant response pat-terns.The linear response approach introduced in this pa-per is especially well suited to the study of the bulk andresonant regime. Crucially we allow for heterogeneousparameters, asymmetric coupling (losses on the line) andarbitrary fluctuation spectra.We find that there appear potentially distinct sets ofnodes from which fluctuations spread strongly into thenetwork (vulnerable nodes) and to which fluctuationsspread easily (excitable nodes) [2, 4]. For autocorrelatedfluctuations that are strongest for small frequencies, thenetworks response is dominated by the bulk, with impor-tant corrections from the resonant regime. Whereas inthe symmetric network case the bulk regime exhibits nonetwork structure, our approach shows that the presenceof an asymmetric effective network Laplacian gives riseto a non-trivial structure in the node vulnerability.We consider a network of coupled nonlinear dynamicalsystems. The dynamical state of each node i can be de- a r X i v : . [ n li n . AO ] F e b scribed by a vector x ( i ) ( t ) containing a certain numberof dynamical variables. The state of the entire dynami-cal system x ( t ) is the vector containing the states of allcomponents. We assume the system has a stable fixpoint x ( t ) = x ◦ . Our goal is now to analyse how the systemresponds to an additive perturbation η ( t ).˙ x = f ( x ) + η ( t ) . (1)If the amplitude of η ( t ) is sufficiently small, the systemstays in the vicinity of the fixpoint x ◦ and we can lin-earize the dynamical system. The linear response of thedeviation δx i ( t ) := x i ( t ) − x ◦ i is given by δx i ( t ) = (cid:90) ∞−∞ χ ( t − t (cid:48) ) η ( t (cid:48) ) dt (cid:48) . (2)Here, χ ( t ) = θ ( t ) e Jt is the network response matrix withthe systems Jacobian J = Df | x = x ◦ at the fixpoint andthe Heaviside function θ ( t ). We quantify the responsesignal by applying the L norm (cid:107) δx ( t ) (cid:107) := (cid:115)(cid:88) i (cid:90) ∞−∞ | δx i ( t ) | dt. (3)Using Parseval’s theorem (cid:107) δx i ( t ) (cid:107) = √ π (cid:107) δ ˆ x i ( ν ) (cid:107) andinserting the Fourier transformed response equation (2)yields (cid:107) δx ( t ) (cid:107) = (cid:115) π (cid:90) ∞−∞ Tr ( ˆ χ ( ν ) S ηη ( ν ) ˆ χ † ( ν )) dν, (4)where S ηη ( ν ) := ˆ η ( ν )ˆ η † ( ν ) is the cross-spectral den-sity matrix. The network response matrix can be de-composed into the response of the single eigenmodesˆ χ ( ν ) = (cid:80) n ˆ χ ( n ) ( ν ), whereˆ χ ( n ) ik ( ν ) = v ( n ) ri v ( n ) lk jν − σ n . Here, v ( n ) l , v ( n ) r are the left and right eigenvectors and σ n the eigenvalues of the systems Jacobian J . An importantclass of systems studied for example in the Master Sta-bility Function approach to networks, are those that arehomogeneous at the linear level, and diffusively coupledon the network. For these systems the Jacobian takeson the following form in terms of the (effective) networkLaplacian: J = A ⊗ B ⊗ L In this case, we can show that the left and right eigen-vectors of the system are given explicitly in terms of thenetwork eigenmodes: v ( n ) = v ( n ) int ⊗ v ( λ ) net λ v ( λ ) net = L v ( λ ) net (5) σ n v ( n ) int = ( A + λB )v ( n ) int . In the following we will focus on the case of a singlenode fluctuation, i.e. there is only one non-zero entry η k ( t ) that drives the system. We denote the L of the re-sponse to a perturbation at node k by (cid:107) δx ( t ) (cid:107) ,k . Sincethe L -norm is a quadratic measure, the mode expan-sion yields not only single mode terms ( n = m ) but alsoterms with cross-mode coupling ( n (cid:54) = m ). For the singlemode terms, we define a spectral excitation factor thatexpresses how much a certain eigenmode of the systemis excited when the power spectral density of the exter-nal fluctuation S η k η k ( ν ) is strong at the eigenfrequency ν = (cid:61) ( σ n ) of this mode: SEF ( n ) = 12 π (cid:90) ∞−∞ S η k η k ( ν ) (cid:60) ( σ n ) + ( ν − (cid:61) ( σ n )) dν . (6)We apply the triangle inequality to the mode expansionand find that neglecting the cross-mode terms ( n (cid:54) = m )defines an upper bound for the L norm of the response (cid:107) δx ( t ) (cid:107) ,k (cid:46) (cid:115)(cid:88) n,i | v ( n ) ri | | v ( n ) lk | SEF ( n ) . (7)Now the core observation is that the L norm is oftenvery close to that bound, which therefore can also serveas a good approximation. As shown in the Supplemen-tal Material, the cross-mode terms will be small if theeigenvalues corresponding to the modes n and m are wellseparated. However, the terms will also be negligible ifthe factors v ( n ) ri v ( m ) ri or ¯ v ( n ) lk ¯ v ( m ) lk are small. For systemswith a Jacobian of the form (5) where eigenvectors fac-torize into network and internal parts, this implies inparticular that these factors are small if the eigenvec-tors have little overlap in their support on the network.Thus cross-mode terms only contribute if the eigenfre-quencies of the two modes are close to each other, andthe corresponding eigenvectors have support at node i and k respectively. While it is hard to give general ruleswhen this will be the case, simulations suggest that thisis typically the case for modes with high frequencies (seeFigure 1). Thus, for auto-correlated noise that has alower power spectral density at higher frequencies, themode-decoupling approximation is more exact.For auto-correlated perturbation signals the integral(6) is generally hard to solve, particularly when the powerspectral density is not known analytically. However, inthat case we could determine the power spectral den-sity from measurements and compute the L -norm semi-analytically. Another possibility is to approximate theintegral by making some general assumptions. In thelow-damping-limit the harmonic modes have a small linewidth compared to the variation of the power spectrum S ηη ( ν ) in that specific frequency range. This meansthat the integral has only significant contributions inthe vicinity of the eigenfrequency ν = (cid:61) ( σ n ). For the n th mode, we can therefore make the approximation S ηη ( ν ) ≈ S ηη ( (cid:61) ( σ n )), i.e. the power spectral densitydoes not depend on the frequency and by integrating the Frequency Frequency responsedecoup-approxpeak-approx FIG. 1.
Power spectral density of fluctuation series, network kernel and network response.
Upper Left: Turbulentpower spectrum of a power fluctuation time series generated with stochastic models of solar and wind power [5, 13]. Lowerleft: Network response spectrum for a certain pair of input and output nodes and the dynamical parameters D = 0 . s and M = 0 . s . The vertical lines correspond to the harmonic eigenfrequencies of the network. Right: Power spectral density ofthe network response as well as its mode-decoupling and peak approximation. Lorentzian function we get
SEF ( n ) ≈ S ηη ( (cid:61) ( σ n ))2 |(cid:60) ( σ n ) | . (8)Fig. 1 shows the spectral response and its differentapproximations for the example system described. It canbe clearly seen that the mode-decoupling approximationand even the peak approximation are relatively closeto the real response in the bulk and in the resonantregime. Only for the local regime (high frequencies)these approximations overestimate the power spectraldensity. However, this regime is only weakly excited bycoloured noise.We now will turn towards an example system todemonstrate the the power of our theoretical approach.We analyse the impact of turbulent fluctuations on powergrids by characterising the vulnerability and excitability of the different nodes in the network. Here, we definethe vulnerability of a node as the strength of networkresponse to power fluctuations at this particular node.Conversely, the excitability of a node is the strength ofthe response at this node given a power fluctuations atanother node.Using our estimates of the L norm we are able to ex-plain three previously observed properties of power grids[2]: I) There is a pronounced fine structure in both vul-nerability and excitability of nodes. II) Losses on thelines lead to a pronounced network structure in the vul-nerability of the nodes, but not in which nodes are ex-cited. III) The vulnerability appears high in parts of the network that are consumer heavy, and within these ar-eas tends to rise the further away from the center of thenetwork the node is.In the following we will provide an analytical expla-nation for these observations for an islanded microgridwith fluctuation renewable infeed. Following [14] droopcontrolled inverters with virtual inertia in their simplestform can be modelled by the swing equation. We includeresistive losses of the power flow on the lines via the con-ductance matrix G ij and power fluctuations δP ( t ) of thepower infeed. This model reads then:˙ φ i = ω i ,M i ˙ ω i = P i + δP i ( t ) − D i ω i − N (cid:88) k =1 P ik ( φ i , φ k ) , (9) P ik = | V i || V k | [ G ik cos( φ i − φ k ) + B ik sin( φ i − φ k )] . We simulate this system on a 100 node network gener-ated by a random growth model for power grids [15].The power fluctuation signal is generated by a combina-tion of stochastic wind and solar power fluctuation mod-els [5, 16]. The power spectrum of the resulting signalis power-lawed with the Kolmogorov exponent of turbu-lence. Further details on the parametrization and thefluctuation modelling can be found in the SupplementalMaterial.We analyse this by looking at pairs of perturbed nodes k and measured nodes i and compute the L norm of thefrequency fluctuations at the measured nodes (cid:107) δω i ( t ) (cid:107) ,k .Following [2], we simulate single node fluctuations inthe full nonlinear system (9) for every pair of input and Output node I n p u t n o d e (a) Nonlinear System Output node (b) Linear System
Output node (c) Peak approximation
Output node (d) Bulk contribution L n o r m o f t h e f r e q u e n c y FIG. 2.
Colour plot of the L norm of the frequency response at single nodes. (a) L norm of the frequencysignals from simulations of the full nonlinear system (9) for each pair of input and output nodes. The dynamic parametersare M = 0 . s and D = 0 . s . (b) L norm of the linearized system. (c) Approximation of the L norm calculated with thepeak approximation of the spectral excitation factor (8). (d) Contribution of the bulk mode to the L norm of the frequencyresponse. output nodes and depict the response strength as a colorcoded matrix plot (Figure 2). There, the horizontal andvertical lines correspond to nodes with large vulnerabil-ity or excitability, respectively. In the following we willanalyse this response fine structure in more detail by ap-plying our linear response theory.For this purpose we first have to linearize the system atthe operation point. Comparing the simulation of the lin-earized system with the full nonlinear system shows thatthe main response pattern remains also in the linearizeddynamics. However, there are some artefacts that are notcaptured by the linear model. These nonlinear effects canobviously not be analysed with our linear theory.For a given time series of a power fluctuation δP ( t ),we can numerically determine its power spectral densityand thereby semi-analytically compute the approxima-tions (7) and (8) for the L norm of the frequency re-sponses. In Figure 2 it can clearly be seen that eventhe very rough peak approximation is able to reproducethe response pattern of the linearized system. The rea-son for this is that the approximations of the spectrumare very good at low frequencies and that for turbulentpower spectra of wind and solar, the deviations in thehigher frequencies are suppressed. This means that byonly knowing the systems eigenvectors and the spectraldensity of the power fluctuations at the systems eigenfre-quencies, we are able to analytically predict the responsestrength at every node in the network. The full networkresponse is a superposition of the responses of every sin-gle mode. The contribution of each eigenmode is deter-mined by the spectral excitation factor and the responsepattern by the left and rights eigenvectors of that mode.When a single mode is dominating the response of thenetwork, the vulnerability and excitability of the nodescan be linked to the left and right eigenvectors of thismode.In the Supplemental Material we show that for homo-geneous and weakly heterogeneous damping-to-inertia-ratios D i /M i in the network there always exists an over- damped eigenmode with an eigenvalue given by σ bulk = − (cid:80) i D i (cid:80) i M i . This mode corresponds to the eigenvalue λ = 0 of theeffective network Laplacian. For homogeneous damping-to-inertia-ratio γ = D i /M i ∀ i , we can show that allthe other modes of the system are harmonic if the alge-braic connectivity of the network is larger than a certainthreshold: | λ | > γ . In this regime the overdamped mode fully determines thebehaviour of the system in the bulk regime (Figure 1) andwe therefore refer to it as the bulk mode. When the alge-braic connectivity is much larger than the threshold wefind that the network response is entirely dominated bythis single mode. For this specific mode we can show,that the right eigenvector is homogeneous but the lefteigenvector has heterogeneous entries. This means thatall nodes are equally excited but certain nodes have muchhigher vulnerability to external fluctuations, even for ho-mogeneous inertia parameters in the network. In theSupplemental Material we show that this effect is due tothe losses of the lines. The resulting dynamical asymme-try corresponds to the continuous horizontal lines in thebulk mode plot (Figure 2). Moreover it can be seen thateven if we are not directly in the bulk regime, the bulkmode can form a substantial contribution to the eventual L norm. This is due to the turbulent power spectra ofwind and solar power fluctuations which are strongest forvery low frequencies (see Figure 1)When the bulk mode is dominating the response ofthe system, the difference of the vulnerability of nodesto power fluctuations is almost entirely determined bythe left eigenvector of this mode. For the case of tree-networks with homogeneous inertia ( M i = M k ) and smallline resistances r ik (cid:28) x ik we derive the following relationfor the entries of the corresponding eigenvector of thenetwork Laplacian:v (0) l,i = (cid:18) − φ ◦ ik ) r ik x ik (cid:19) v (0) l,k . From this it follows that the entries of the left bulk modeeigenvector are monotonously increasing with the powerflow. Consequently, when this mode is dominating, the L norm is increasing along the power flow. This meansthat the nodes located at the sinks of the flow are themost vulnerable to power fluctuations. In Figure 3 it canclearly be seen that the network branch where the poweris flowing from the center towards the outlying nodes ismuch more vulnerable than the network branches wherethe power is flowing towards the center. This explainsFigure 5 in [2] showing the relation between the vulnera-bility of nodes and the closeness centrality of the network. L norm of the frequency response0.250 0.275 0.300 0.325 0.350 0.375 0.400 0.425 L norm of the frequency response FIG. 3.
Vulnerability of nodes in the power flow net-work for bulk dominated dynamics.
The dynamic pa-rameters are D = 0 . s and M = 0 . s . Nodes are colouredaccording to the L -norm of the system response for fluctua-tions at this node. The link arrows indicate the direction ofthe power flow and the width of the links is scaled accordingto the amount of flow. We have seen that the analysis of the power spec-trum of frequency fluctuations provides a powerful toolto understand how fluctuating infeeds spread in com-plex power grids. This is especially true for the reso-nant regime in which no good analytic approximationswere previously known. Crucially we did not require theassumption of lossless power lines. This allowed us to unveil novel patterns in the vulnerability of nodes, andreveal network patterns in the bulk response of coupledoscillators.This work focuses on single node fluctuations. How-ever, the formulation of the response in terms of powerspectra also provides an elegant starting point for futureinvestigations of multi-node fluctuations. In that case notonly the auto-correlation but also the cross-correlation offluctuations will play a crucial role for understanding thedynamical interactions.With some straightforward assumptions on the powerspectrum of the power fluctuations, we were able to de-rive reasonable approximations for the network responsespectrum and thereby for the L norm of frequencyfluctuations. It should be noted that the approxima-tions we used here, while successful at describing the ob-served phenomena in the particular regime of interest, aremostly uncontrolled and rather rough. Assuming specificanalytic models of the power spectrum of input fluctua-tions would enable us to derive much more sophisticatedapproximations from the mode expansion of the networkresponse matrix.The approximations were numerically validated for theimportant example case of an islanded microgrid withhigh renewable penetration. With these we were ableto explain structures in the node vulnerability that havebeen observed in numerical simulations [2, 4] and to showthat these vulnerability patterns originate from the losseson power lines.Considering the generality of our theoretical approach,it should be mentioned that the application to powergrids is not limited to fluctuation of RES but might alsobe the basis for studying the impact of demand fluctua-tions on grids.For tree-like grids the location of vulnerable nodes isrelated to the power flow throughout the network. Inparticular, fluctuations at net power flow sinks resultsin strong frequency fluctuations at all nodes in the net-work. For microgrids grids with a a very unbalancedpower production, we therefore expect the consumerheavy branches to be much more vulnerable to turbulentinfeed of renewable power. In the Supplemental Materialwe show that the distinct relation between the flow di-rection and vulnerability is true for a much broader classof dynamical flow networks. For second order node dy-namics we are able to derive a specific criterion for theparameter regime where this is the case. For general sys-tems, the relation is also true if the flow is stabilizing thedynamics at the fixpoint. Deriving a more solid criterionfor this effect to occur is a topic for future research. ACKNOWLEDGMENTS
The work presented was partially funded by the BmBF(Grant No. 03EK3055A) and the Deutsche Forschungs-gemeinschaft (Grant No. KU 837/39-1 / RA 516/13-1). [1] X. Zhang, S. Hallerberg, M. Matthiae, D. Witthaut, andM. Timme, Science Advances , eaav1027 (2019).[2] S. Auer, F. Hellmann, M. Krause, and J. Kurths, Chaos , 10.1063/1.5001818 (2017).[3] F. Hellmann, P. Schultz, P. Jaros, R. Levchenko, T. Kap-itaniak, J. Kurths, and Y. Maistrenko, Nature Commu-nications , 1 (2020).[4] S. Auer, The Stability and Control of Power Grids withHigh Renewable Energy Share , Ph.D. thesis, Humboldt-Universit¨at zu Berlin, Humboldt-Universit¨at (insgesamt)(2018).[5] K. Schmietendorf, J. Peinke, and O. Kamps, The Euro-pean Physical Journal B , 222 (2017).[6] S. Tamrakar, M. Conrath, and S. Kettemann, Scientificreports (2018).[7] B. Sch¨afer, M. Matthiae, X. Zhang, M. Rohden,M. Timme, and D. Witthaut, Phys. Rev. E , 060203(2017).[8] M. Tyloo, T. Coletta, and P. Jacquod, Phys. Rev. Lett. , 084101 (2018).[9] T. Coletta, B. Bamieh, and P. Jacquod, in (IEEE,2018) pp. 6163–6167.[10] M. Tyloo, L. Pagnier, and P. Jacquod, Science Advances , eaaw8359 (2019).[11] H. Haehne, K. Schmietendorf, S. Tamrakar, J. Peinke,and S. Kettemann, Physical Review E , 050301 (2019).[12] S. Kettemann, Phys. Rev. E , 062311 (2016).[13] M. Anvari, G. Lohmann, M. W¨achter, P. Milan,E. Lorenz, D. Heinemann, M. R. R. Tabar, and J. Peinke,New Journal of Physics , 063027 (2016).[14] J. Schiffer, D. Goldin, J. Raisch, and T. Sezi, in Decisionand Control (CDC), 2013 IEEE 52nd Annual Conferenceon (IEEE, 2013) pp. 2334–2339.[15] P. Schultz, J. Heitzig, and J. Kurths, The EuropeanPhysical Journal Special Topics , 2593 (2014).[16] M. Anvari, B. Werther, G. Lohmann, M. W¨achter,J. Peinke, and H.-P. Beck, Solar Energy , 735 (2017).[17] T. Boßmann and I. Staffell, Energy , 1317 (2015).[18] S. Auer, F. Steinke, W. Chunsen, A. Szabo, and R. Sol-lacher, in PES Innovative Smart Grid Technologies Con-ference Europe (ISGT-Europe), 2016 IEEE (IEEE, 2016)pp. 1–6.[19] P. C. Sen,
Principles of electric machines and power elec-tronics (John Wiley & Sons, 2007).[20] J. Schiffer, D. Zonetti, R. Ortega, A. M. Stankovi´c,T. Sezi, and J. Raisch, Automatica , 135 (2016).[21] E. A. A. Coelho, P. C. Cortizo, and P. F. D. Garcia, IEEETransactions on Industry Applications , 533 (2002).[22] S. Asmussen, Applied probability and queues , Vol. 51(Springer Science & Business Media, 2008).[23] P. Milan, M. W¨achter, and J. Peinke, Physical reviewletters , 138701 (2013).[24] H. Hurst, Hydrological Sciences Journal , 13 (1956).[25] E. Alessio, A. Carbone, G. Castelli, and V. Frappietro,The European Physical Journal B-Condensed Matter andComplex Systems , 197 (2002).[26] T. Preis, P. Virnau, W. Paul, and J. J. Schneider, NewJournal of Physics , 093024 (2009).[27] A. Carbone, G. Castelli, and H. Stanley, Physical ReviewE , 026105 (2004). SUPPLEMENTAL MATERIALParametrization of the Microgrid Model
The microgrid model case is kept at a conceptual levelto study the effect of local fluctuations on dynamic gridstability and isolate the influence of the network struc-ture.Germany has 4,500 MV distribution networks thatconnect 500,000 LV distribution networks [17]. Thus themicrogrid is chosen as a network of 100 nodes to representan average German grid at medium-voltage (MV) level.The MV level is a good testing case for modelling powergrids with a high renewable energy share, since most PVpower plants are connected to low-voltage- (LV) or MVlevels. An islanded microgrid must be internally power-balanced and not connected to a higher grid level.Being power balanced, we assume that there are 50net producers and 50 net consumers with P i = ± . M W power infeed before losses. The power infeeds are chosenhomogeneously to focus on topology and network effectsin the model. As there is no connection to upper gridlevels, losses are compensated locally at each node, andthe net power infeed is given by ˜ P i = ( P i + P loss /N ).Mathematically this is equivalent to switching to the co-rotating frame.The network topology is generated by a randomgrowth model for power grids [15]. We have chosena parametrization such that we get radial grids whichis a typical structure for distribution grids. The gridparametrization follows from the voltage level. The lineimpedance for typical MV grids with 20kV base voltageequals Z = R + jX ≈ (0 . . j ) Ω /km [18]. For sim-plicity all power, voltage and impedance values are trans-formed into per unit with a base voltage of 20kV and abase power of 1MW, which are typical values for MVgrids [18, 19]. The absolute impedance of each line scaleswith the geographic distance l between linked nodes andconsequently differs per link. The average line length,according to [18], is 23 . τ p and the droop control parameter k p fromgrid-forming inverters: M = τ p /k p , D = 1 /k p , ∀ i with i = 1 , .., N . Standard parameters for the droop and time constants of grid-forming inverters are in the range k p = 0 . . . . s − and τ p = 0 . . . . s [14, 21]. In thelow inertia case, with only a few low powered grid form-ing inverters at each node, we assume a weakly reacting,strongly smoothed system.In the simulations intermittent time series for solar andwind power fluctuations were generated by a clear skyindex model, based on a combination of a Langevin anda jump process, developed in [13], and a Non-MarkovianLangevin type model developed in [5], respectively. Thepower fluctuation δP ( t ) is a combination of the signalsgenerated with these models for wind and solar powerfluctuations ∆ d ( t ) = 0 . δP W ( t ) + 0 . δP S ( t ) . (10)Both the PDF of the fluctuations and their incrementtime series are fat tailed (the tails are not exponen-tially bounded [22]). The power generation from windand solar power plants has a power spectrum that ispower-lawed with the Kolmogorov exponent of turbu-lence [13, 23]. Thus, these time series show long-termtemporal correlations [24–27]. Mode Expansion of the Response Function
In general, the Jacobian has distinct left and righteigenvectors v ( n ) l , v ( n ) r for every eigenvalue σ n . The leftand right eigenvectors are orthonormal to each other butnot orthonormal within themselves v ( n ) l · v ( m ) r = δ nm . When the eigenvectors of J are linearly independent, itcan be decomposed to J = Q · diag( σ n ) · Q − , where Q = ( v (1) r , ..., v ( N ) r ) , Q − = v (1) l ... v ( N ) l . In an alternative representation, the matrix decomposi-tion can be written as J = (cid:88) n σ n v ( n ) r v ( n ) l , where v ( n ) r v ( n ) l is the outer product of right and left eigen-vectors. Using the fact that the matrix exponential e J has the same eigenvectors as the Jacobian itself, the re-sponse function can be decomposed to χ ( t ) = (cid:88) n χ ( n ) ( t ) = (cid:88) n θ ( t ) e σ n t v ( n ) r v ( n ) l . With this, we decomposed the response of the system intothe response of its single eigenmodes. The Fourier trans-form of a single element of these mode response matricesis given by ˆ χ ( n ) ik ( ν ) = v ( n ) ri v ( n ) lk jν − σ n . Dynamics in Flow Networks
In many network dynamical systems, the nodes arecoupled by a flow: The flow function F i → j ( x ij ) representsthe flow going from node i into the link between i and j .It depends on the state variable difference x ij := x i − x j .We separate the flow function in an antisymmetric anda symmetric part: F i → j ( x ij ) = F ai → j ( x ij ) + F si → j ( x ij ) . The antisymmetric part F ai → j corresponds to the flowthat gets transmitted to j and the symmetric part F si → j corresponds to the loss on the link.The flow is positive in the direction from higher tolower state variable: F ai → j ( x ij ) > x ij > . When the states at both nodes are equal, there is no flowon the link: F ai → j (0) = 0 . For the derivative of the antisymmetric part it followsthat: (cid:0) F ai → j (cid:1) (cid:48) ( x ij ) ≥ . In order to have only losses and no gains on the links,the symmetric part has to be positive: F si → j ( x ij ) ≥ . When the flow is zero, there are no losses: F si → j (0) = 0 . For the derivative of the symmetric part it follows that: (cid:0) F si → j (cid:1) (cid:48) ( x ij ) ≥ x ij > , (cid:0) F si → j (cid:1) (cid:48) ( x ij ) ≤ x ij < . In the linearization of a dynamical flow network at afixpoint x ◦ , we get an effective network Laplacian of theform L = L a + L s ,L a = − δ ij (cid:88) l w ail + w aij ,L s = − δ ij (cid:88) l w sil + w sij , where we introduced the asymmetic and symmetricweights w sij = (cid:0) F ai → j (cid:1) (cid:48) ( x ◦ ij ) ,w aij = (cid:0) F si → j (cid:1) (cid:48) ( x ◦ ij ) . Note that the weights for the antisymmetric flow are sym-metric whereas the weights for the symmetric flow areantisymmetric: w sij = w sji ,w aij = − w aji . Thus, the symmetric part of the flow i.e. the losses onthe link, correspond to an asymmetry in the effectiveLaplacian.The Laplacian always has one zero eigenvalue as can beseen by multiplying with the homogeneous vector v (0) r,j =(1 , ..., T : L v (0) r,j = L a v (0) r,j + L s v (0) r,j = 0 . The left eigenvector is inhomogeneous for L s (cid:54) = 0. Wehave the condition:0 = (cid:88) i v (0) l,i L ij = (cid:88) i (cid:104) w aij (v (0) l,i + v (0) l,j ) + w sij (v (0) l,i − v (0) l,j ) (cid:105) =: (cid:88) i S ij . In tree-like networks this equation gives a strict relationfor the eigenvector entries of two neighbouring nodes. Wecan see this by starting at the nodes with degree one. Atthese nodes there is only one summand S ij which there-fore has to be zero. Now the summands are antisym-metric S ij = − S ji and therefore we know that the cor-responding summand S ji in the condition for the neigh-bouring node is also zero. By going up the tree structureto nodes of higher degree we can eliminate the summandscorresponding to all previously visited nodes. Doing this,we see that in the above equations not only the sum isequal to zero but every single summand has to be zeroitself and therefore:v (0) l,i = v (0) l,j · w sij − w aij w sij + w aij We always have w sij >
0. The antisymmetric weights arepositive w aij > x ◦ ij > w aij < x ◦ ij <
0. It follows thatv (0) i < v (0) j for x ij > , v (0) i > v (0) j for x ij < . Therfore, the entries of the left eigenvector increase alongthe direction of the flow.For dynamical flow networks with homogeneous nodedynamics the Jacobian takes the following form in termsof the (effective) network Laplacian: J = A ⊗ B ⊗ L In this case, the left and right eigenvectors of the systemare given explicitly in terms of the network eigenmodes: v ( λ,µ ) = v ( λ,µ ) int ⊗ v ( λ ) net λ v ( λ ) net = L v ( λ ) net σ ( λ,µ ) v ( λ,µ ) int = ( A + λB )v ( λ,µ ) int . This is an eigenvector as: Jv ( λ,µ ) = ( A ⊗ B ⊗ L )v ( λ,µ ) int ⊗ v ( λ ) net = A v ( λ,µ ) int ⊗ v ( λ ) net + B v ( λ,µ ) int ⊗ L v ( λ ) net = A v ( λ,µ ) int ⊗ v ( λ ) net + λB v ( λ,µ ) int ⊗ v ( λ ) net = ( A + λB )v ( λ,µ ) int ⊗ v ( λ ) net = σ ( λ,µ ) v ( λ,µ ) For the zero eigenmode of the Laplacian, the eigenvalueof the Jacobian is given by the eigenvalues of the internaldynamics: σ (0 ,µ ) v (0 ,µ ) int = A v (0 ,µ ) int From a control theory perspective, it is desirable thatthe internal dynamics of an uncoupled system is non-oscillatory and therefore we can assume that (cid:61) ( σ (0 ,µ ) ).Now for a certain class of flow network dynamics thenetwork coupling introduces an oscillatory behavior andfor certain parameter regimes, all modes except the bulkmode oscillate. For systems with second order internaldynamics, this is the case iftr( A + λB ) < A + λB ) ∀ λ (cid:54) = 0 . For correlated perturbations where the excitation of thehigh frequency is suppressed, the zero eigenmode willdominate the response of the system.More generally, the zero eigenmode of the Laplacianis also dominating the response for systems where thecoupling is stabilizing the system at the fixpoint (cid:60) ( σ ( λ,µ ) ) ≤ (cid:60) ( σ (0 ,µ ) ) . Mode-Decoupling Approximation
In the following we will focus on fluctuations η k ( t ) at asingle dynamical variable. The L norm of the responseat another variable δx i ( t ) is given by (cid:107) δx i ( t ) (cid:107) ,k = 12 π (cid:90) ∞−∞ | ˆ χ ik ( ν ) | S η k η k ( ν ) dν. Inserting the mode decomposition of the spectral re-sponse function yields (cid:107) δx i ( t ) (cid:107) ,k = 12 π (cid:88) n,m (cid:90) ∞−∞ ˆ χ ( n ) ik ( ν ) ¯ˆ χ ( m ) ik ( ν ) S η k η k ( ν ) dν. We will now make an approximation where we discardall cross mode terms (cid:107) δx i ( t ) (cid:107) ,k (cid:47) π (cid:88) n (cid:90) ∞−∞ | ˆ χ ( n ) ik ( ν ) | S η k η k ( ν ) dν, where the squared mode response functions have theshape of a Lorentzian function | ˆ χ ( n ) ik ( ν ) | = | v ( n ) ri | | v ( n ) lk | (cid:60) ( σ n ) + ( ν − (cid:61) ( σ n )) . First of all, it follows from the triangle inequality thatthe approximation above is an upper bound of the fullexpression. Moreover, we will show that the cross modeterms are negligible under certain conditions.First, we will show that the cross mode terms are sup-pressed if the spectral gap between the modes is largecompared to the damping of the modes. We introducethe shorthands γ n = (cid:60) ( σ n ) and ν n = (cid:61) ( σ n ). The maxi-mum of the diagonal mode contributions to the networktransfer function is given by ν = ν n max ν | ˆ χ ( n ) ik ( ν ) | = | v ( n ) ri | | v ( n ) lk | γ n . For the cross mode contributions we getˆ χ ( n ) ik ( ν ) ˆ χ ( m ) ik ( ν ) = v ( n ) ri v ( n ) lk v ( m ) ri v ( m ) lk [ γ n − j ( ν − ν n )][ γ m + j ( ν − ν m )] . For simplicity we assume that both modes have similardamping γ := γ n = γ m . Then the magnitude of thedenominator reaches its extrema at0 = (2 ν − ν n − ν m ) (cid:0) γ + ( ν − ν n )( ν − ν m ) (cid:1) . If γ ≥ ( ν n − ν m ) , there is only one minimum of thedenominator at ν min = ν n + ν m . If we instead have γ < ( ν n − ν m ) , this minimum turns into a local maxi-mum and the new local minima are given by ν min = ν n + ν m ± (cid:114) ( ν n − ν m ) − γ . Therefore, in this case, the maximum of the magnitudeof the off-diagonal term is given byˆ χ ( n ) ik ( ν min ) ˆ χ ( m ) ik ( ν min ) = v ( n ) ri v ( n ) lk v ( m ) ri v ( m ) lk | γ || ν n − ν m | We see that if | ν n − ν m | (cid:29) | ρ | the off-diagonal modeterms are suppressed relative to the diagonal terms. Ifthe network modes are not overdamped, this is true formost pairs of modes.0 Swing Equation with Line Losses
The swing equation with power fluctuations and lossesis given by˙ φ i = ω i ,M i ˙ ω i = P i + δP i ( t ) − D i ω i − N (cid:88) k =1 P ik ( φ i , φ k ) ,P ik ( φ i , φ k ) = | V i || V k | [ G ik cos( φ i − φ k ) + B ik sin( φ i − φ k )] . Here, G ik and B ik are the entries of the conductance andsusceptance matrices which are given by the real andimaginary part of the complex admittance matrix Y ik = δ il (cid:88) l y il − y ik ,y ik = ( r ik + jx ik ) − . Here r ik and x ik are resistance and reactance of thelines. In the swing equation the voltage amplitudes areconstant. We assume all buses in the network to be atthe same voltage level | V i | = | V k | = V ref .The system is measured in the co-rotating frame of thesynchronous state ω i = ω k ∀ i, k of the unperturbed sys-tem δP = 0. In the co-rotating frame, the synchronousstate is a fixpoint and we denote the constant phase an-gles as φ ◦ i . This fixpoint is assumed to be linearly stable,i.e. the real parts of all eigenvalues are smaller or equalto zero. Linearising around the fixpoint yields (cid:20) δ ˙ φ ˙ ω (cid:21) = (cid:20) N × N N × N M − L − M − D (cid:21) (cid:20) δφω (cid:21) + (cid:20) N M − δP ( t (cid:48) ) (cid:21) , were M = diag( M i ), D = diag( D i ) are the parametermatrices and L is the weighted Laplacian matrix at thefixpoint. The Laplacian can be decomposed into a sym-metric and an antisymmetric part L ik = L sik + L aik whichare given by L sik = − δ ik (cid:88) l V ref B il cos( φ ◦ il ) + V ref B ik cos( φ ◦ ik ) ,L aik = δ ik (cid:88) l V ref G il sin( φ ◦ il ) − V ref G ik sin( φ ◦ ik ) , Here we used the shorthand φ ◦ ik := φ ◦ i − φ ◦ k . Let σ n bethe Jacobian eigenvalues and v ( n ) l = (cid:104) u ( n ) l w ( n ) l (cid:105) , v ( n ) r = (cid:34) u ( n ) r w ( n ) r (cid:35) , the corresponding left and right eigenvectors. Let us firstassume a constant ratio between the droop and inertiaconstants in the network D i = γM i ∀ i . The eigenequa- tions for the Jacobian matrix are then given by (cid:20) N × N N × N M − L − γ · N × N (cid:21) (cid:34) u ( n ) r w ( n ) r (cid:35) = σ n (cid:34) u ( n ) r w ( n ) r (cid:35) , (cid:104) u ( n ) l w ( n ) l (cid:105) (cid:20) N × N N × N M − L − γ · N × N (cid:21) = σ n (cid:104) u ( n ) l w ( n ) l (cid:105) . From this we can derive the following equations w ( n ) l · M − L = σ n ( σ n + γ ) w ( n ) l M − L · w ( n ) r = σ n ( σ n + γ ) w ( n ) r . These equations imply that in case of a homogeneousratio of droop and inertia the Jacobian eigenvalues canbe calculated from the eigenvalues λ of M − L by σ ( λ, ± ) = − γ ± (cid:114) γ λ. Further, the lower halves of the Jacobian eigenvectors areproportional to the eigenvectors of the inertia weightedLaplacian M − L v ( λ ) r ∝ w ( λ, ± ) r , v ( λ ) l ∝ w ( λ, ± ) l Using this the linear frequency response to the powerfluctuation δP ( t ) is given by ω ( t ) = (cid:88) n v ( n ) r v ( n ) l (cid:90) t −∞ e σ n ( t − t (cid:48) ) M − δP ( t (cid:48) ) dt (cid:48) . Bulk Mode
When the network is connected and if the systemhas no further symmetry the inertia weighted Lapla-cian has exactly one zero eigenvalue λ = 0. The cor-responding eigenvalues of the Jacobian are σ (0 , +) = 0and σ (0 , − ) = − γ . Here, σ (0 , +) corresponds to the sym-metry of homogeneous phase shifts that have no impacton the dynamics, whereas σ (0 , − ) corresponds to homo-geneous frequency shifts leading to an exponentially de-caying response of the nodes’ frequencies with rate γ .We can therefore conclude that in the case of homoge-neous droop-to-inertia-ratio, the system always has oneoverdamped eigenmode. Moreover, if the algebraic con-nectivity λ of the network fulfils the condition | λ | > γ , then mode of the eigenvalue σ (0 , − ) is the only over-damped mode in the system.The right eigenvector corresponding to this mode ishomogeneous, i.e. v (0) r,i = v (0) r,k ∀ i, k , and consequentiallythe linear response to this mode is the equal for allfrequencies ω i . Because of this bulk-like behaviour the1bulk mode σ bulk := σ (0 , − ) .The left eigenvector of the bulk mode is generally het-erogeneous and fulfils the condition0 = (cid:88) i v (0) l,i M − i L ik = (cid:88) i v (0) l,i M − i L sik + (cid:88) i v (0) l,i M − i L aik ∀ k. With the definitions of the symmetric and antisymmetricpart we can calculate the two terms and get0 = (cid:88) i V ref (cid:34) B ik cos( φ ◦ ik ) (cid:32) v (0) l,i M i − v (0) l,k M k (cid:33) − G ik sin( φ ◦ ik ) (cid:32) v (0) l,i M i + v (0) l,k M k (cid:33)(cid:35) ∀ k. In tree-like networks these equations give a strict condi-tion for the eigenvector entries of the neighbouring nodes.Since there are no loops in these networks, there exits noadditional degree of freedom for choosing v (0) l,k once a sin-gle entry is fixed to a specific value. This can be easilyseen when we set the eigenvector entry at a node withdegree one, calculate the entry of its neighbour and con-tinue going up the tree structure. Mathematically thismeans that in the above equations not only the sum isequal to zero but every single summand has to be zeroitself. From this condition, we can derive the followingequation for the ratio of the eigenvector entries of twoneighbouring nodesv (0) l,i v (0) l,k = M i M k (cid:32) G ik B ik tan( φ ◦ ik )1 − G ik B ik tan( φ ◦ ik ) (cid:33) . For a small ohmic resistance are small compared to thereactance r ik (cid:28) x ik we can Taylor expand the expressionand get v (0) l,i ≈ M i M k (cid:18) − φ ◦ ik ) r ik x ik (cid:19) v (0) l,k . From this, it follows that for tree-like networks evenin the case of homogeneous inertia parameters M i = M k ∀ i, k the left eigenvector corresponding to the bulkmode has non-identical entries if there are losses on theline ( r ik > Parameter Heterogeneity
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii M ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii M ii = − (cid:80) ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii M ii = − (cid:80) ii D ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii M ii = − (cid:80) ii D ii (cid:80) ii
We have shown above that in a grid with heterogeneousdroop-to-inertia-ratio D i = γM i there always exists anoverdamped mode that is directly related to the powerlosses on the lines. In the following we will show thatif we allow for small heterogeneities in the parameters D i = ( γ + δγ i ) M i this bulk mode will still persist. Inthe heterogeneous case it is not possible to analyticallyderive the Jacobian eigenvalues since the sub-matricesin the Jacobian do not commute any more. However,for small heterogeneities we can use matrix perturbationtheory to calculate the change δσ n of the eigenvalues dueto the heterogeneity. The Jacobian of the heterogeneouscase can be separated into a homogeneous part J hom = (cid:20) N × N N × N M − L − γ · N × N (cid:21) and an inhomogeneous part δJ = (cid:20) N × N N × N N × N − diag( δγ i ) (cid:21) . To linear order, the change of the eigenvalues is given by δσ n = v ( n ) l δJ v ( n ) r = − (cid:88) i w li δγ i w ri . For the homogeneous case, the eigenvalue of the bulkmode is σ bulk = − γ . For the sake of simplicity we assumethat the line losses are very small G ik B ik (cid:28) w li ∝ M i . A possible representation of theleft and right eigenvectors of the bulk mode is then givenby u li = 0 , w li = M i (cid:80) i M i , u li = − γ , w li = 1 ∀ i. Therefore the linear order approximation of the bulkmode eigenvalue is δσ bulk = − (cid:80) i δγ i M i (cid:80) i M i . This means that for small heterogeneities, the bulk modeeigenvalue is real and thus the mode is still overdamped.Moreover, if the sum of heterogeneities vanishes, theeigenvalue is the same as in the homogeneous case. Thiscan be generalized to the statement that the dampingcoefficient of the bulk mode is given by the ratio of totaldroop to total inertia in the grid σ bulk = − γ − (cid:80) i δγ ii M ii (cid:80) ii M ii = − (cid:80) ii D ii (cid:80) ii M ii