Optimal placement of inertia and primary control : a matrix perturbation theory approach

The increasing penetration of inertialess new renewable energy sources reduces the overall mechanical inertia available in power grids and accordingly raises a number of issues of grid stability over short to medium time scales. It has been suggested that this reduction of overall inertia can be compensated to some extent by the deployment of substitution inertia - synthetic inertia, flywheels or synchronous condensers. Of particular importance is to optimize the placement of the limited available substitution inertia, to mitigate voltage angle and frequency disturbances following a fault such as an abrupt power loss. Performance measures in the form of H2-norms have been recently introduced to evaluate the overall magnitude of such disturbances on an electric power grid. However, despite the mathematical conveniance of these measures, analytical results can be obtained only under rather restrictive assumptions of uniform damping ratio, or homogeneous distribution of inertia and/or primary control in the system. Here, we introduce matrix perturbation theory to obtain analytical results for optimal inertia and primary control placement where both are heterogeneous. Armed with that efficient tool, we construct two simple algorithms that independently determine the optimal geographical distribution of inertia and primary control. These algorithms are then implemented on a model of the synchronous transmission grid of continental Europe. We find that the optimal distribution of inertia is geographically homogeneous but that primary control should be mainly located on the slow modes of the network, where the intrinsic grid dynamics takes more time to damp frequency disturbances.


I. INTRODUCTION
The penetration of new renewable energy sources (RES) such as photovoltaic panels and wind turbines is increasing in most electric power grids worldwide. In their current configuration, these energy sources are essentially inertialess, and their increased penetration leads to low inertia situations in periods of high RES production [1]. This raises important issues of power grid stability, which is of higher concern to transmission system operators than the volatility of the RES productions [2], [3]. The substitution of traditional productions based on synchronous machines with inertialess RES may in particular lead to geographically inhomogeneous inertia profiles. It has been suggested to deploy substitution inertia -synthetic inertia, flywheels or synchronous condensersto compensate locally or globally for missing inertia. Two related question naturally arise, which are (i) where is it safe to substitute synchronous machines with inertialess RES and (ii) where is it optimal to distribute substitution inertia ? Problem (ii) has been investigated in small power grid models with up to a dozen buses, optimizing the geographical distribution of inertia against cost functions based on eigenvalue damping ratios [4], H p -norms [5], [6] or RoCoF [7], [8] and frequency excursions [8]. Investigations of problem (i) on large power grids emphasized the importance of the geographical extent of the slow network modes [9]. Numerical optimization can certainly be performed for any given network on a case-bycase basis, however it is highly desirable to shed light on the problem with analytical results. So far, such results have been either restricted to small systems or derived under the assumption that damping and inertia parameters, or their ratio, are homogeneous. In this manuscript we go beyond these assumptions and construct an approach applicable to large power grids with inhomogeneous independent damping and inertia parameters.
Inspired by theoretical physics, we introduce matrix perturbation theory [10] as an analytical tool to tackle this problem. That method is widely used in quantum physics, where it delivers approximate solutions to complex, perturbed problems, extrapolated from known, exact solutions of integrable problems [11]. The approximation is valid as long as the difference between the two problems is small and it makes sense to consider the full, complex problem as a perturbation of the exactly soluble, simpler problem. The procedure is spectral in nature. It identifies a small, dimensionless parameter in which eigenvalues and eigenvectors of the perturbed problem can be systematically expanded in a series not unlike a Taylor-expansion about the unperturbed, integrable solution. Depending on the value of that small parameter, the series can be truncated at low orders already and still deliver rather accurate results. In the context of electric power grids, the method was applied for instance in Ref. [12], where quadratic performance measures similar to those discussed below were calculated following a line fault, starting from the eigenvalues and eigenvectors of the network Laplacian before the fault.
In this paper, we apply matrix perturbation theory [10] to calculate performance measures following an abrupt power loss in the case of a transmission power grid with geographically inhomogeneous inertia and damping (primary control) parameters. Our perturbation theory is an expansion in two parameters which are the maximal deviations δm and δd of the rotational inertia and damping parameters from their average values m and d. The approach is valid as long as these local deviations are small, | δm/m | < 1, | δd/d | < 1. These conditions tolerate in principle that inertia and damping parameters vanish or are twice as large as their average on some buses. The main step forward brought about by our approach is that we are able to derive analytical results without relying on the often used homogeneity assumptions that damping, inertia or their ratio is constant -assumptions which are not satisfied in real electric power grids. Our main results are given in Theorems 1 and 2 below, which formulate arXiv:1906.06922v1 [math.OC] 17 Jun 2019 algorithms for optimal placement of local inertia and damping parameters. The spectral decomposition approach used here has recently drawn the attention of a number of groups and has been used to calculate performance measures in power grids and consensus algorithms e.g. in [7], [13], [14], [15].
The article is organized as follows. Section II deals with the case where inertia and primary control are uniformly distributed in the system. The performance measure that quantifies system disturbances is introduced and we calculate its value for abrupt power losses. In Section III we apply matrix perturbation theory to calculate the sensitivities of our measure in local variations of inertia and primary control. Section IV presents the optional placement of inertia and primary control in the case of weak inhomogeneity. In Section V we apply our optimal placements to the continental European grid. Section VI concludes our article.

II. HOMOGENEOUS CASE
We are interested in the dynamical response of an electric power grid to a disturbance such as an abrupt power loss. To that end, we consider the power system dynamics in the lossless line approximation, which is a standard approximation used for high voltage transmission grids [16]. That dynamics is governed by the swing equations, The inertia-weighted matrix L M is real and symmetric, therefore it can be diagonalized with an orthogonal matrix U , the α th row of which gives the components u α,i , i = 1, . . . N of the α th eigenvector u α of L M . The diagonal matrix Λ = diag({λ 1 = 0, λ 2 , · · · , λ N }) contains the eigenvalues of L M with λ α < λ α+1 . For connected networks only the smallest eigenvalue λ 1 vanishes, which follows from the zero row and column sum property of the Laplacian matrix L M . Rewriting (3) in the basis diagonalizing L M gives where δθ M = U ξ. This change of coordinates is nothing but a spectral decomposition of angle deviations δθ M into their components in the basis of eigenvectors of L M . These components are cast in the vector ξ. The formulation (5) of the problem makes it clearer that, if Γ is a multiple of identity, the problem can be recast as a diagonal ordinary differential equation problem that can be exactly integrated. This is done below in (11), and provides an exact solution about which we will construct a matrix perturbation theory in the next sections.
Proof: The proof goes along the lines of the diagonalization procedure proposed in [7], [12], [17]. Equation (5) can be rewritten as where P = U M −1/2 δP and 0 N ×M is the N × M matrix of zeroes. The matrix H 0 is block-diagonal up to a permutation of rows and columns [12], and can easily be diagonalized block by block, where each 2 × 2 block corresponds to one of the eigenvalues λ α of L M . The α th block is diagonalized by the transformation with the eigenvalues µ (0) α± of the α th block, The two rows (columns) of T L α (T R α ) give the nonzero components of the two left (right) eigenvectors t of H 0 . Following this transformation, (7) reads d dt The solutions of (11) are Inserting (12) back into (8), one finally finds (6) which proves the proposition.

B. Performance measure
We want to mitigate disturbances following an abrupt power loss. To that end, we use performance measures which evaluate the overall disturbance magnitude over time and the whole power grid. Performance measures have been proposed, which can be formulated as L 2 and squared H 2 norms of linear systems [6], [7], [12], [17], [18], [19], [20], [21], [22], [23] and are time-integrated quadratic forms in the angle, δθ, or frequency, ω, deviations. Here we focus on frequency deviations and use the following performance measure whereω = (ω sys , ω sys , . . . ω sys ) is the instantaneous average frequency vector with components It is straightforward to see that M reads when rewritten in the eigenbasis of L M , once one notices that the first eigenvector of L M (the one with zero eigenvalue) has Proposition 2. For an abrupt power loss, δP (t) = δP Θ(t) on a single bus labeled b, δP i = δ ib δP , and with an homogeneous damping ratio, Γ = γ 1 with the N ×N identity matrix 1, in terms of the eigenvalues λ α and the components u αb of the eigenvectors u α of L M .
Note that we introduced the subscript b to indicate that the fault is localized on that bus only. The power loss is modeled as which, when summed over α > 1 gives (16). α of the Laplacian L. In that case, the performance measure reads where the superscript (0) refers to inertia homogeneity. This expression has an interesting graph theoretic interpreation. We recall the definitions of the resistance distance Ω ij between two nodes on the network, the associated centrality C j and the generalized Kirchhoff indices Kf p [22], [24], where L † is the Moore-Penrose pseudo inverse of L. With these definitions, one can show that [15], [22], [25] α>1 by using the spectral representation of the resistance distance [26], [27] Because Kf 1 is a global quantity characterizing the network, it follows from (18) with (22) that, when inertia and primary control are homogeneously distributed in the system, the disturbance magnitude as measured by M is larger for disturbances on peripheral nodes [9], [28].

III. MATRIX PERTURBATION
The previous section treats the case where inertia and primary control are uniformly distributed in the system. Our goal is to lift that restriction and to obtain M b when some mild inhomogeneities are present. We parametrize these inhomogeneities by writing with the average m and γ and the maximum deviation amplitudes δm and δγ of inertia and damping ratio. Inhomogeneities are determined by the coefficients −1 ≤ a i , r i ≤ 1 with i r i = i a i = 0 which are determined following a minimization of the performance measure M b of (13). In the following two paragraphs we construct a matrix perturbation theory to linear order in the inhomogeneity parameters δm, and δγ to calculate the performance measure

A. Inhomogeneity in inertia
When inertia is inhomogeneous, but the damping ratios remain homogeneous, the system dynamics and M b are still given by (7) and (16). However, the eigenvectors of the inertiaweighted Laplacian matrix L M differ from those of L and consequently M b is no longer equal to M (0) b . In general there is no simple way to diagonalize L M , but one expects that if the inhomogeneity is weak, then the eigenvalues and eigenvectors of L M only slightly differ from those of m −1 L, which allows to construct a perturbation theory.
Assumption 1 (Weak inhomogeneity in inertia). The deviations δm r i of the local inertias m i are all small compared to their average m. We write 1 is a small, dimensionless parameter.
To linear order in µ, the series expansion of L M reads with In this form, the inertia-weighted Laplacian matrix L M is given by the sum of an easily diagonalizable matrix, m −1 L, and a small perturbation matrix, (µ/m) V 1 . Truncating the expansion of L M at this linear order gives an error of order ∼ µ 2 , which is small under Assumption 1. Matrix perturbation theory gives approximate expressions for the eigenvectors u α and eigenvalues λ α of L M in terms of those (u with From (16), (27) and (28), the first-order approximation of M b in µ reads Proposition 3. For an abrupt power loss, δP (t) = δP Θ(t) on a single bus labeled b, δP i = δ ib δP , and under Assumption 1, the susceptibilites ρ i ≡ ∂M b /∂r i are given by Proof: Taking the derivative of (31) with respect to r i , with λ (1) α and u (1) αb given in (29) and (30), one gets The first term in the square bracket in (33) gives where we used β u This terms therefore exactly cancels out with the last two terms in the square bracket in (33) and one obtains

B. Inhomogeneity in damping ratios
Equation (6) gives exact solutions to the linearized dynamical problem defined in (5), under the assumption of homogeneous damping ratio, m i /d i ≡ γ. In this section we lift that constraint and write γ i = γ + δγ a i . With inhomogeneous damping ratios, (7) becomes which differs from (7) only through the additional term −δγV 2 . Under the assumption that that the dimensionless parameter g ≡ δγ/γ 1, this additional term gives only small corrections to the unperturbed problem of (7), and we use matrix perturbation theory to calculate these corrections in a polynomial expansion in g.
Assumption 2 (Weak inhomogeneity in damping ratios). The deviations δγ a i of the damping ratio γ i from their average γ are all small compared to their average. We write Γ = γ 1 + g diag {a i } , where g ≡ δγ/γ 1 is a small, dimensionless parameter.
We want to integrate (36) using the spectral approach that provided the solutions (12). In principle this requires to know the eigenvalues and eigenvectors of H in (36), which is not possible in general, because V 2 does not commute with Λ. When g is small enough, the eigenvalues and eigenvectors are only slightly altered [10] and can be systematically calculated order by order in a polynomial expansion in g. We therefore follow a perturbative approach which expresses solutions to (36) in such a polynomial expansion in g. Formally, one has, for the eigenvalues µ αs and for the left and right eigenvectors t L,R αs of H where the m = 0 terms are given by the eigenvalues, µ αs , and the left and right eigenvectors, t (0)L,R αs , of the matrix H 0 in (7), corresponding to homogeneous inertia. In order for the sums in (37) and (38) to converge, a necessary condition is that g < 1. The task is to calculate the terms µ 1, one expects that only few, low order terms already give a good estimate of the eigenvalues and eigenvectors of H. In this manuscript, we calculate the first-order corrections, m = 1. They are given by formulas similar to (29) and (30), where indicates that the sum runs over (β, s ) = (α, s). One obtains with V 2;αβ = i a i u αi u βi .
Remark 3. By definition, −1 ≤ V 2;αα ≤ 1. Therefore, (42) indicates, among others, that when the parameters {a i } are correlated (anticorrelated) with the square components {u 2 αi } for some α then that mode is more strongly (more weakly) damped. Accordingly, Theorem 2 will distribute the set {a i } to increase the damping of the slow modes of H. Proposition 4. For an abrupt power loss, δP (t) = δP Θ(t) on a single bus labeled b, δP i = δ ib δP , and under Assumption 2,ξ α (t) reads, to leading order in g, where s α = sin(f α t/2) and c α = cos(f α t/2), and P α and f α are defined below (6).
The proof is based on (42) to (44) and is given in Appendix A.
Proposition 5. For an abrupt power loss, δP (t) = δP Θ(t) on a single bus labeled b, δP i = δ ib δP , and under Assumption 2, the susceptibilities α i ≡ ∂M b /∂a i are given by Proof: From (45), to first order in g, one has Taking the derivative of (47) with respect to a i with the definition of V 2;αβ given below (44), and summing over α > 1 one obtains (46).

Remark 4.
We have found numerically that the second term is generally much smaller than the first one and gives only marginal corrections to our optimized solution.

IV. OPTIMAL PLACEMENT OF INERTIA AND PRIMARY
CONTROL In general it is not possible to obtain closed-form analytical expressions for the parameters a i and r i determining the optimal placement of inertia and primary control. Simple optimization algorithms can however be constructed that determine how to distribute these parameters to minimize M b . Theorems 1 and 2 give two such algorithms for optimization under Assumption 1 and 2 respectively. Additionally, Conjecture 1 proposes an algorithm for optimization under both Assumption 1 and 2.
The proof is in Appendix A.
Proof: With Proposition 5 and M = m1, we get (46). The proof is the same as the one for Theorem 1 given in We next conjecture an algorithmic combined linear optimization treating simultaneously Assumptions 1 and 2. The difficulty is that for fixed total inertia and damping, one must have (25), the second condition requires i a i r i = 0. This is a quadratic, nonconvex constraint, which makes the problem nontrivial to solve. The following conjecture presents an algorithm that starts from the distribution {a i } and {r i } from Theorems 1 and 2 and orthogonalizes them while trying to minimize the related increase in M b .

Conjecture 1 (Combined linear optimization).
For an abrupt power loss, under Assumptions 1 and 2, the optimal placement of a fixed total amount of inertia i m i = mN and primary control i d i = dN that minimize M b is obtained as follows.
1) Compute the parameters r i and a i from Theorems 1 and 2. a) If N is odd, align the zeros of {r i } and {a i } . Let i r0 and i a0 be the indexes of these zeros. Their new common index is Interchange the parameter values r ir0 ↔ r i align and a ir0 ↔ a i align . b) If N is even, do nothing 2) If n ≡ i r i a i = 0, the optimization is done.
3) Find the set I = {i | sgn(r i a i ) = sgn(n)}. To reach i r i a i → 0, our strategy is to set to zero some elements of I. Since however i a i = i r i = 0 must be conserved, this must be accompanied by a simultaneous change of some other parameter.

4)
Find the pair (a i1 , a i2 = −a i1 ) or (r i1 , r i2 = −r i1 ) ∈ I × I which, when sent to (0, 0), induce the smallest increase of the objective function M b . Send it to (0, 0). Because the pair has opposite sign, this does not affect the condition i a i = i r i = 0. 5) go to step # 2.
It is not at all guaranteed that the algorithm presented in Conjecture 1 is optimal, however, numerical results to be presented below indicate that it works well.
The optimization considered so far focused on a single fault on bus labeled b. We are interested, however, in finding the optimal distribution of inertia and/or primary control for all possible faults. To that end we introduce the following global vulnerability measure where the sum runs over all generator buses. The vulnerability measure V gives a weighted average over all possible fault positions, with the weight η b accounting for the probability that a fault occurs at b and δP b accounting for its potential intensity as given, e.g. by the rated power of the generator at bus b. For equiprobable fault locations and for the same power loss everywhere, η b ≡ 1, with Remark 2, it is straightforward to see that ∂V/∂r i = 0 + O(µ 2 ). Therefore, to leading order, there is no benefit in scaling up the inertia anywhere. On the other hand, with Remark 5, we get ∂V/∂a i = −gδP 2 α>1 u The corresponding optimal placement of primary control can be obtained with Theorem 2, from which we observe that the damping ratios are increased for the buses with large squared components u (0)2 αi of the slow modes of L -those with the smallest λ (0) α . These modes are displayed in Fig 1. One concludes that, with a non-weighted vulnerability measure, η b ≡ 1 in (50), an homogeneous inertia location is a local optimum for V, for which damping parameters need to be increased primarily on peripheral buses. (2) (3) (4) u max -u max 0

V. NUMERICAL INVESTIGATIONS
We illustrate our main results on a model of the synchronous power grid of continental Europe. The network has 3809 nodes, among them 618 generators, connected through 4944 lines. For details of the model and its construction we refer the reader to [9], [25]. To connect to the theory presented above, we remove inertialess buses through a Kron reduction [29] and uniformize the distribution of inertia to m i = 29.22MWs 2 , and primary control d i = 12.25MWs. This guarantees that the total amounts of inertia and primary control are kept at their initial levels. At the end of the previous section we argued that an homogeneous distribution of inertia, together with primary control increased on the slowest eigenmodes of the network Laplacian minimize the global vulnerability measure V of (50) for η b ≡ 1. This conclusion is confirmed numerically in Fig. 2 (a) and (b). The optimal placement of primary control displayed in panel (b) decreases V by more than 12% with respect to the homogeneous case.
Setting η b ≡ 1 in (50) is convenient mathematically, however it treats all faults equally, regardless of their impact. One may instead adapt η b to obtain inertia and primary control distributions that reduce the impact of the strongest faults with largest M b . We do this in two different ways, first with and second with The corresponding geographical distributions of inertia and primary control redistribution parameters r i and a i parameters are shown in Fig. 2 (c) and (d) and Fig. 2 (e) and (f) respectively. Compared to the choice η b ≡ 1 [ Fig. 2 (a) and (b)], we see rather small differences. More importantly, the impact of various choices of η b on M b is almost negligible, as can be seen in Fig. 2 (g). In all cases, our optimization algorithm reduces first and foremost the impact of the strongest faults, with little or no influence on the faults that have little impact on grid stability. It is finally interesting to note that our three choices of η b are close to being optimal, especially when considering the strongest faults. This can be seen in Fig. 2 (g), where the purple line shows the maximal obtainable reduction, when inertia and primary control distributions are optimized individually fault by fault, i.e. with a different redistribution for each fault. The inset of Fig. 2 (g) shows in particular that for the strongest fault, the three considered choices of η b lead to reductions in M b that are very close to the maximal one. We conclude that rather generically, inertia is optimally distributed homogeneously, while primary control should be preferentially located on the slow modes of the grid Laplacian.

VI. CONCLUSION
To find the optimal placement of inertia and primary control in electric power grids with limited such resources is a problem of paramount importance. Here, we have made what we think is an important step forward in constructing a perturbative analytical approach to this problem. In this approach, both inertia and primary control are limited resources, as should be. Most importantly, our method goes beyond the usually made assumption of constant inertia to damping ratio. In our approach inertia and primary control can vary spatially independently from one another. Our results suggest that the optimal inertia distribution is close to homogeneous over the whole grid, but that primary control should be reinforced on buses located on the support of the slower modes of the network Laplacian. Further work should try to extend the approach to the next order in perturbation theory. Work along those lines is in progress.

ACKNOWLEDGMENT
This work was supported by the Swiss National Science Foundation under grants PYAPP2 154275 and 200020 182050.

Proof of Proposition 4:
The proof follows the same steps as for Proposition 1. The calculations are rather tedious, though algebraically straightforward. In the following, we sketch the calculational steps. Assuming that H can be diagonalized as t R µ t L , where µ ≡ diag({µ αs }) and t R,L are matrices containing the right and left eigenvectors of H, the problem is resolved by: 1) changing variables χ ≡ t L [ξ ξ ] to diagonalize (36), aṡ 2) solving (52) as 3) Obtainingξ α via the inverse transformation [ξ ξ ] = t R χ. These three steps are carried out with the approximate expressions t R,L = t R,L(0) + gt R,L(1) and µ α± = µ obtained with the first order in g corrections presented in (42)-(44). One gets with where (45) is obtained from (54) by applying trigonometric identities.
Proof of Theorem 1: To leading order in µ = δm/m, this optimization problem is equivalent to the following linear programming problem [30] It is solved by the Lagrange multipliers method, with the Lagrangian function where ε i and ε 0 are Lagrange multipliers. We get The solution must satisfy the Karush-Kuhn-Tucker (KKT) conditions [30], in particular the complementary slackness (CS) condition which imposes that either ε i = 0 or r i = ±1 , ∀i. The former choice leads generally to a contradiction. From (61) and dual feasibility condition, one gets This imposes that r i = −sgn(ε 0 + ρ i ). To ensure that i r i = 0 is satisfied, ε 0 is set to minus the median value of ρ i . If the number of bus N is odd, the r i corresponding to the median value of ρ i is set to zero. Philippe Jacquod (M'16) received the Diplom degree in theoretical physics from the ETHZ, Zürich, Switzerland, in 1992, and the PhD degree in natural sciences from the University of Neuchâtel, Switzerland, in 1997. He is a professor with the engineering department, University of Applied Sciences of Western Switzerland, Sion, Switzerland, with a joint appointment with the Department of Quantum Matter Physics, University of Geneva, Switzerland. From 2003 to 2005 he was an assistant professor with the theoretical physics department, University of Geneva, Switzerland and from 2005 to 2013 he was a professor with the physics department, University of Arizona, Tucson, USA. His main research topics is in power systems and how they evolve as the energy transition unfolds. He has published about 100 papers in international journals, books and conference proceedings.