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

TL;DR

This paper introduces a matrix perturbation theory approach to optimize the placement of inertia and primary control in power grids with heterogeneous parameters, improving stability and disturbance mitigation.

Contribution

It develops analytical tools and algorithms for optimal placement of inertia and primary control considering heterogeneity, advancing beyond previous uniform assumptions.

Findings

Optimal inertia distribution is geographically homogeneous.

Primary control should be placed on slow network modes.

Algorithms effectively determine optimal control placement.

Abstract

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.

Figures

Equations

$$m_i{\dot{\omega }}_i+d_i{\omega }_i=P_i-\sum _j{B}_{ij}\sin ({\theta }_i-{\theta }_j),$$

$$\mathbf{M}\dot{\mathbf{\omega }}+\mathbf{D}\mathbf{\omega }=\mathbf{\delta P}-\mathbf{L\delta \theta },$$

$${\dot{\mathbf{\omega }}}_{\mathbf{M}}+\underset{\mathbf{\Gamma }}{\underbrace{{\mathbf{M}}^{-1}\mathbf{D}}}{\mathbf{\omega }}_{\mathbf{M}}+\underset{{\mathbf{L}}_{\mathbf{M}}}{\underbrace{{\mathbf{M}}^{-1/2}\mathbf{L}{\mathbf{M}}^{-1/2}}}\mathbf{\delta }{\mathbf{\theta }}_{\mathbf{M}}={\mathbf{M}}^{-1/2}\mathbf{\delta P},$$

$${\mathbf{L}}_{\mathbf{M}}={\mathbf{U}}^{\top }\mathbf{\Lambda }\mathbf{U}$$

$$\ddot{\mathbf{\xi }}+\mathbf{U}\mathbf{\Gamma }{\mathbf{U}}^{\top }\dot{\mathbf{\xi }}+\mathbf{\Lambda }\mathbf{\xi }=\mathbf{U}{\mathbf{M}}^{-1/2}\mathbf{\delta P},$$

$${\dot{\xi }}_{\alpha }(t)$$

$$\frac{{\rm d}}{{\rm d}t}\left[\!\begin{array}[]{c}\bm{\xi}\\ \bm{\dot{\xi}}\end{array}\!\right]=\underbrace{\left[\!\begin{array}[]{cc}\mathbb{0}_{N\times N}&\!\!\!\!\mathbb{1}\\ \!\!\!\!-\bm{\Lambda}&\!\!\!\!\!\!-\gamma\,\mathbb{1}\end{array}\!\right]}_{\bm{H}_{0}}\!\left[\!\begin{array}[]{c}\bm{\xi}\\ \bm{\dot{\xi}}\end{array}\!\right]+\bigg{[}\!\begin{array}[]{c}\mathbb{0}_{N\times 1}\\ \bm{\mathcal{P}}\end{array}\!\bigg{]}\,,$$

$$\displaystyle\left[\!\begin{array}[]{c}\chi^{(0)}_{\alpha+}\\ \chi^{(0)}_{\alpha-}\end{array}\!\right]$$

$$\displaystyle\left[\!\begin{array}[]{c}\xi_{\alpha}\\ \dot{\xi}_{\alpha}\end{array}\!\right]$$

$${\mu }_{\alpha \pm }^{(0)}=-\frac{1}{2}(\gamma \mp i{f}_{\alpha }).$$

$$\frac{{\rm d}}{{\rm d}t}\left[\!\!\begin{array}[]{c}\chi^{(0)}_{\alpha+}\\ \chi^{(0)}_{\alpha-}\end{array}\!\!\right]=\left[\!\!\begin{array}[]{cc}\mu_{\alpha+}^{(0)}&0\\ 0&\mu_{\alpha-}^{(0)}\end{array}\!\!\right]\!\left[\!\!\begin{array}[]{c}\chi^{(0)}_{\alpha+}\\ \chi^{(0)}_{\alpha-}\end{array}\!\!\right]+\frac{i}{f_{\alpha}}\left[\!\!\begin{array}[]{c}-\mathcal{P}_{\alpha}\\ \phantom{-}\mathcal{P}_{\alpha}\end{array}\!\!\right]\,.$$

$${\chi }_{\alpha \pm }^{(0)}$$

$$\mathcal{M}=\intop_{0}^{\infty}\big{(}\bm{\omega}^{\top}-\bm{\bar{\omega}}^{\top}\big{)}\bm{M}\big{(}\bm{\omega}-\bm{\bar{\omega}}\big{)}{\rm d}t\,,$$

$$\omega_{\rm sys}(t)=\sum_{i}m_{i}\omega_{i}(t)\Big{/}\sum_{i}m_{i}\,.$$

$$\mathcal{M}=\underset{0}{\overset{\infty }{\int }}\sum _{\alpha >1}{\dot{\xi }}_{\alpha }^2(t)\mathrm{d}t,$$

$${\mathcal{M}}_b=\frac{\delta P^2}{2\gamma m_b}\sum _{\alpha >1}\frac{u_{\alpha b}^2}{{\lambda }_{\alpha }},$$

$$\underset{0}{\overset{\infty }{\int }}{\dot{\xi }}_{\alpha }^2(t)\mathrm{d}t=\frac{\delta P^2{u}_{\alpha b}^2}{2\gamma m_b{\lambda }_{\alpha }},\alpha >1,$$

$${\mathcal{M}}_b^{(0)}=\frac{\delta P^2}{2\gamma }\sum _{\alpha >1}\frac{u_{\alpha b}^{(0)2}}{{\lambda }_{\alpha }^{(0)}},$$

$${\Omega }_{ij}$$

$$C_j$$

$$K{f}_p$$

$$\sum _{\alpha >1}\frac{u_{\alpha b}^{(0)2}}{{\lambda }_{\alpha }^{(0)}}=C_b^{-1}-N^{-2}K{f}_1,$$

$$\Omega_{ib}=\sum_{\alpha>1}\big{(}u_{\alpha i}^{(0)}-u_{\alpha b}^{(0)}\big{)}^{2}\big{/}\lambda_{\alpha}^{(0)}\,.$$

$$m_i$$

$$d_i$$

$${\mathbf{L}}_{\mathbf{M}}$$

$${\lambda }_{\alpha }$$

$${\mathbf{u}}_{\alpha }$$

$${\lambda }_{\alpha }^{(1)}$$

$${\mathbf{u}}_{\alpha }^{(1)}$$

$${\mathcal{M}}_b$$

$${\rho }_i=-\frac{\mu \delta P^2}{\gamma N}\sum _{\alpha >1}\frac{u_{\alpha b}^{(0)}{u}_{\alpha i}^{(0)}}{{\lambda }_{\alpha }^{(0)}}.$$

$$\frac{\partial {\mathcal{M}}_b}{\partial r_i}$$

Full text

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 condensers – to 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], $\mathcal{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-by-case 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 $\delta m$ and $\delta 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, $|\,\delta m/m\,|<1$, $|\,\delta 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 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,

[TABLE]

which determine the time-evolution of the voltage angles $\theta_{i}$ and frequencies $\omega_{i}=\dot{\theta}_{i}$ at each of the $N$ buses labelled $i$ in the power grid in a rotating frame such that $\omega_{i}$ measures the angle frequency deviation to the rated grid frequency of 50 or 60 Hz. Each bus is characterized by inertia, $m_{i}$, and damping, $d_{i}$, parameters and $P_{i}$ is the active power injected ($P_{i}>0$) or extracted ($P_{i}<0$) at bus $i$. We introduce the damping ratio $\gamma_{i}\equiv d_{i}/m_{i}$. Buses are connected to one another via lines with susceptances $B_{ij}$. Stationary solutions $\{\theta_{i}^{(0)}\}$ are power flow solutions determined by $P_{i}=\sum_{j}B_{ij}\sin(\theta_{i}^{(0)}-\theta_{j}^{(0)})$. Under a change in active power $P_{i}\rightarrow P_{i}+\delta P_{i}$, linearizing the dynamics about such a solution with $\theta_{i}(t)=\theta_{i}^{(0)}+\delta\theta_{i}(t)$ gives, in matrix form,

[TABLE]

where $\bm{M}={\rm diag}(\{m_{i}\})$, $\bm{D}={\rm diag}(\{d_{i}\})$ and voltage angles and frequencies are cast into vectors $\bm{\delta\theta}$ and $\bm{\omega}\equiv\bm{\delta\dot{\theta}}$. The Laplacian matrix $\bm{L}$ has matrix elements $L_{ij}=-B_{ij}\cos(\theta_{i}^{(0)}-\theta_{j}^{(0)})$, for $i\neq j$ and $L_{ii}=\sum_{k}B_{ik}\cos(\theta_{i}^{(0)}-\theta_{k}^{(0)})$.

II-A Exact solution for homogeneous damping ratio

When the damping ratio is constant, $d_{i}/m_{i}=\gamma_{i}=\gamma$, $\forall i$, (2) can be integrated exactly [7, 12]. To see this we first transform angle coordinates as $\bm{\delta\theta}=\bm{M}^{-1/2}\bm{\delta\theta_{M}}$ to obtain

[TABLE]

where we introduced the diagonal matrix $\bm{\Gamma}={\rm diag}(\{d_{i}/m_{i}\})\equiv{\rm diag}(\{\gamma_{i}\})$. The inertia-weighted matrix $\bm{L_{M}}$ is real and symmetric, therefore it can be diagonalized

[TABLE]

with an orthogonal matrix $\bm{U}$, the $\alpha^{\rm th}$ row of which gives the components $u_{\alpha,i}$, $i=1,\ldots N$ of the $\alpha^{\rm th}$ eigenvector $\bm{u}_{\alpha}$ of $\bm{L_{M}}$. The diagonal matrix $\bm{\Lambda}={\rm diag}(\{\lambda_{1}=0,\lambda_{2},\cdots,\lambda_{N}\})$ contains the eigenvalues of $\bm{L}_{M}$ with $\lambda_{\alpha}<\lambda_{\alpha+1}$. For connected networks only the smallest eigenvalue $\lambda_{1}$ vanishes, which follows from the zero row and column sum property of the Laplacian matrix $\bm{L_{M}}$. Rewriting (3) in the basis diagonalizing $\bm{L}_{M}$ gives

[TABLE]

where $\bm{\delta\theta_{M}}=\bm{U}^{\top}\bm{\xi}$. This change of coordinates is nothing but a spectral decomposition of angle deviations $\bm{\delta\theta_{M}}$ into their components in the basis of eigenvectors of $\bm{L}_{\bm{M}}$. These components are cast in the vector $\bm{\xi}$. The formulation (5) of the problem makes it clearer that, if $\bm{\Gamma}$ 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 (23), and provides an exact solution about which we will construct a matrix perturbation theory in the next sections.

Proposition 1** (Unperturbed evolution).**

For an abrupt power loss, $\bm{\delta P}(t)=\bm{\delta P}\,\Theta(t)$, with the Heaviside step function defined by $\Theta(t>0)=1$, $\Theta(t<0)=0$, and with homogeneous damping ratio, $\bm{\Gamma}=\gamma\,\mathbb{1}$ with the $N\times N$ identity matrix $\mathbb{1}$, the frequency coordinates $\dot{\xi}_{\alpha}$ evolve independently as

[TABLE]

where $f_{\alpha}=\sqrt{4\lambda_{\alpha}-\gamma^{2}}$ and $\mathcal{P}_{\alpha}=\sum_{i}u_{\alpha i}\,\delta P_{i}/m_{i}^{1/2}$.

This result generalizes Theorem III.3 of [14].

Proof:

The proof goes along the lines of the diagonalization procedure proposed in [7, 12, 17]. Equation (5) can be rewritten as

[TABLE]

where $\bm{\mathcal{P}}=\bm{U}\bm{M}^{-1/2}\bm{\delta P}$ and $\mathbb{0}_{N\times M}$ is the $N\times M$ matrix of zeroes. The matrix $\bm{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\times 2$ block corresponds to one of the eigenvalues $\lambda_{\alpha}$ of $\bm{L_{M}}$. The $\alpha^{\rm th}$ block is diagonalized by the transformation

[TABLE]

with the eigenvalues $\mu_{\alpha\pm}^{(0)}$ of the $\alpha^{\rm th}$ block,

[TABLE]

The two rows (columns) of $\bm{T}_{\alpha}^{L}$ ($\bm{T}_{\alpha}^{R}$) give the nonzero components of the two left (right) eigenvectors $\bm{t}_{\alpha\pm}^{(0)L}$ ($\bm{t}_{\alpha\pm}^{(0)R}$) of $\bm{H}_{0}$. Following this transformation, (7) reads

[TABLE]

The solutions of (23) are

[TABLE]

Inserting (24) back into (14), one finally finds (6) which proves the proposition. ∎

II-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 $\mathcal{L}_{2}$ and squared $\mathcal{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, $\bm{\delta\theta}$, or frequency, $\bm{\omega}$, deviations. Here we focus on frequency deviations and use the following performance measure

[TABLE]

where $\bm{\bar{\omega}}=(\omega_{\rm sys},\omega_{\rm sys},\ldots\omega_{\rm sys})^{\top}$ is the instantaneous average frequency vector with components

[TABLE]

It is straightforward to see that $\mathcal{M}$ reads

[TABLE]

when rewritten in the eigenbasis of $\bm{L}_{\bm{M}}$, once one notices that the first eigenvector of $\bm{L}_{\bm{M}}$ (the one with zero eigenvalue) has components $u_{1i}=\sqrt{m_{i}}/\sqrt{\sum_{j}m_{j}}$.

Proposition 2**.**

For an abrupt power loss, $\bm{\delta P}(t)=\bm{\delta P}\,\Theta(t)$ on a single bus labeled $b$, $\delta P_{i}=\delta_{ib}\,\delta P$, and with an homogeneous damping ratio, $\bm{\Gamma}=\gamma\,\mathbb{1}$ with the $N\times N$ identity matrix $\mathbb{1}$,

[TABLE]

in terms of the eigenvalues $\lambda_{\alpha}$ and the components $u_{\alpha b}$ of the eigenvectors $\bm{u}_{\alpha}$ of $\bm{L}_{\bm{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 $P_{i}=P_{i}^{(0)}-\delta P_{i}\,\Theta(t)$ with $\delta P_{i}=\delta_{ib}\,\delta P$ with the Kronecker symbol $\delta_{ib}=1$ if $i=b$ and 0 otherwise.

Proof:

Equation (6) straightforwardly gives

[TABLE]

which, when summed over $\alpha>1$ gives (28). ∎

Remark 1**.**

For homogeneous inertia coefficients, $\bm{M}=m\mathbb{1}$, the eigenvectors and eigenvalues of the inertia-weighted Laplacian $\bm{L}_{\bm{M}}$ defined in (3) are given by $\bm{u}_{\alpha}=\bm{u}_{\alpha}^{(0)}$, and $\lambda_{\alpha}=m^{-1}\lambda_{\alpha}^{(0)}$, in terms of the eigenvectors $\bm{u}_{\alpha}^{(0)}$ and eigenvalues $\lambda_{\alpha}^{(0)}$ of the Laplacian $\bm{L}$. In that case, the performance measure reads

[TABLE]

where the superscript (0) refers to inertia homogeneity. This expression has an interesting graph theoretic interpreation. We recall the definitions of the resistance distance $\Omega_{ij}$ between two nodes on the network, the associated centrality $C_{j}$ and the generalized Kirchhoff indices $Kf_{p}$ [22, 24],

[TABLE]

where $\bm{L}^{\dagger}$ is the Moore–Penrose pseudo inverse of $\bm{L}$. With these definitions, one can show that [15, 22, 25]

[TABLE]

by using the spectral representation of the resistance distance [26, 27]

[TABLE]

Because $Kf_{1}$ is a global quantity characterizing the network, it follows from (30) with (34) that, when inertia and primary control are homogeneously distributed in the system, the disturbance magnitude as measured by $\mathcal{M}_{b}^{(0)}$ 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 $\mathcal{M}_{b}$ when some mild inhomogeneities are present. We parametrize these inhomogeneities by writing

[TABLE]

with the average $m$ and $\gamma$ and the maximum deviation amplitudes $\delta m$ and $\delta\gamma$ of inertia and damping ratio. Inhomogeneities are determined by the coefficients $-1\leq a_{i},r_{i}\leq 1$ with $\sum_{i}r_{i}=\sum_{i}a_{i}=0$ which are determined following a minimization of the performance measure $\mathcal{M}_{b}$ of (25). In the following two paragraphs we construct a matrix perturbation theory to linear order in the inhomogeneity parameters $\delta m$, and $\delta\gamma$ to calculate the performance measure $\mathcal{M}_{b}=\mathcal{M}_{b}^{(0)}+\sum_{i}r_{i}\rho_{i}+\sum_{i}a_{i}\alpha_{i}+\mathcal{O}(\delta m^{2},\delta\gamma^{2})$. This requires to calculate the susceptibilities $\rho_{i}\equiv\partial\mathcal{M}_{b}/\partial r_{i}$ and $\alpha_{i}\equiv\partial\mathcal{M}_{b}/\partial a_{i}$.

III-A Inhomogeneity in inertia

When inertia is inhomogeneous, but the damping ratios remain homogeneous, the system dynamics and $\mathcal{M}_{b}$ are still given by (7) and (28). However, the eigenvectors of the inertia-weighted Laplacian matrix $\bm{L_{M}}$ differ from those of $\bm{L}$ and consequently $\mathcal{M}_{b}$ is no longer equal to $\mathcal{M}_{b}^{(0)}$. In general there is no simple way to diagonalize $\bm{L_{M}}$, but one expects that if the inhomogeneity is weak, then the eigenvalues and eigenvectors of $\bm{L_{M}}$ only slightly differ from those of $m^{-1}\bm{L}$, which allows to construct a perturbation theory.

Assumption 1** (Weak inhomogeneity in inertia).**

The deviations $\delta m\,r_{i}$ of the local inertias $m_{i}$ are all small compared to their average $m$. We write \bm{M}=m\big{[}\mathbb{1}+\mu\,{\rm diag}\big{(}\{r_{i}\}\big{)}\big{]}, where $\mu\equiv\delta m/m\ll 1$ is a small, dimensionless parameter.

To linear order in $\mu$, the series expansion of $\bm{L_{M}}$ reads

[TABLE]

with \bm{V}_{1}=-\big{(}\bm{R}\bm{L}+\bm{L}\bm{R}\big{)}/2 and $\bm{R}={\rm diag}(\{r_{i}\})$. In this form, the inertia-weighted Laplacian matrix $\bm{L_{M}}$ is given by the sum of an easily diagonalizable matrix, $m^{-1}\bm{L}$, and a small perturbation matrix, $(\mu/m)\,\bm{V}_{1}$. Truncating the expansion of $\bm{L_{M}}$ at this linear order gives an error of order $\sim\mu^{2}$, which is small under Assumption 1.

Matrix perturbation theory gives approximate expressions for the eigenvectors $\bm{u}_{\alpha}$ and eigenvalues $\lambda_{\alpha}$ of $\bm{L_{M}}$ in terms of those ($\bm{u}_{\alpha}^{(0)}$ and $\lambda_{\alpha}^{(0)}$) of $\bm{L}$ [10]. To leading order in $\mu$ one has

[TABLE]

with

[TABLE]

From (28), (39) and (40), the first-order approximation of $\mathcal{M}_{b}$ in $\mu$ reads

[TABLE]

Proposition 3**.**

For an abrupt power loss, $\bm{\delta P}(t)=\bm{\delta P}\,\Theta(t)$ on a single bus labeled $b$, $\delta P_{i}=\delta_{ib}\,\delta P$, and under Assumption 1, the susceptibilites $\rho_{i}\equiv\partial\mathcal{M}_{b}/\partial r_{i}$ are given by

[TABLE]

Proof:

Taking the derivative of (43) with respect to $r_{i}$, with $\lambda_{\alpha}^{(1)}$ and $u_{\alpha b}^{(1)}$ given in (41) and (42), one gets

[TABLE]

The first term in the square bracket in (45) gives

[TABLE]

where we used $\sum_{\beta}u_{\beta i}^{(0)}u_{\beta b}^{(0)}=\delta_{ib}$. This terms therefore exactly cancels out with the last two terms in the square bracket in (45) and one obtains

[TABLE]

The argument of the double sum in (47) is odd under permutation of $\alpha$ and $\beta$, therefore only terms with $\beta=1$ survive. With $u_{1i}^{(0)}=1/\sqrt{N}$, one finally obtains (44). ∎

Remark 2**.**

By summing over every fault locations $b$, one gets $\sum_{b}\rho_{i}(b)=0$. This follows from the properties of the eigenvector $\bm{u}_{\alpha}^{(0)}$.

III-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}\equiv\gamma$. In this section we lift that constraint and write $\gamma_{i}=\gamma+\delta\gamma\,a_{i}$. With inhomogeneous damping ratios, (7) becomes

[TABLE]

which differs from (7) only through the additional term $-\delta\gamma\bm{V}_{2}$ with $\bm{V}_{2}=\bm{U}\bm{A}\,\bm{U}^{\top}$, $\bm{A}={\rm diag}(\{a_{i}\})$. Under the assumption that that the dimensionless parameter $g\equiv\delta\gamma/\gamma\ll 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 $\delta\gamma\,a_{i}$ of the damping ratio $\gamma_{i}$ from their average $\gamma$ are all small compared to their average. We write \bm{\Gamma}=\gamma\big{[}\mathbb{1}+g\,{\rm diag}\big{(}\{a_{i}\}\big{)}\big{]}, where $g\equiv\delta\gamma/\gamma\ll 1$ is a small, dimensionless parameter.

We want to integrate (48) using the spectral approach that provided the solutions (24). In principle this requires to know the eigenvalues and eigenvectors of $\bm{H}$ in (48), which is not possible in general, because $\bm{V}_{2}$ does not commute with $\bm{\Lambda}$. 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 (48) in such a polynomial expansion in $g$. Formally, one has, for the eigenvalues $\mu_{\alpha s}$ and for the left and right eigenvectors $\bm{t}_{\alpha s}^{L,R}$ of $\bm{H}$

[TABLE]

where the $m=0$ terms are given by the eigenvalues, $\mu_{\alpha s}^{(0)}$, and the left and right eigenvectors, $\bm{t}_{\alpha s}^{(0)L,R}$, of the matrix $\bm{H}_{0}$ in (7), corresponding to homogeneous inertia. In order for the sums in (49) and (50) to converge, a necessary condition is that $g<1$. The task is to calculate the terms $\mu_{\alpha s}^{(m)}$ and $\bm{t}_{\alpha s}^{(m)L,R}$ with $m=1,2,...$. When $g\ll 1$, one expects that only few, low order terms already give a good estimate of the eigenvalues and eigenvectors of $\bm{H}$. In this manuscript, we calculate the first-order corrections, $m=1$. They are given by formulas similar to (41) and (42),

[TABLE]

where ${\overline{\sum}}$ indicates that the sum runs over $(\beta,s^{\prime})\neq(\alpha,s)$. One obtains

[TABLE]

with $\bm{V}_{2;\alpha\beta}=\sum_{i}a_{i}\,u_{\alpha i}\,u_{\beta i}$.

Remark 3**.**

By definition, $-1\leq\bm{V}_{2;\alpha\alpha}\leq 1$. Therefore, (60) indicates, among others, that when the parameters $\{a_{i}\}$ are correlated (anticorrelated) with the square components $\{u_{\alpha i}^{2}\}$ for some $\alpha$ 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 $\bm{H}$.

Proposition 4**.**

For an abrupt power loss, $\bm{\delta P}(t)=\bm{\delta P}\,\Theta(t)$ on a single bus labeled $b$, $\delta P_{i}=\delta_{ib}\,\delta P$, and under Assumption 2, $\dot{\xi}_{\alpha}(t)$ reads, to leading order in $g$,

[TABLE]

where $s_{\alpha}=\sin(f_{\alpha}t/2)$ and $c_{\alpha}=\cos(f_{\alpha}t/2)$, and $\mathcal{P}_{\alpha}$ and $f_{\alpha}$ are defined below (6).

The proof is based on (60) to (62) and is given in Appendix A.

Proposition 5**.**

For an abrupt power loss, $\bm{\delta P}(t)=\bm{\delta P}\,\Theta(t)$ on a single bus labeled $b$, $\delta P_{i}=\delta_{ib}\,\delta P$, and under Assumption 2, the susceptibilities $\alpha_{i}\equiv\partial\mathcal{M}_{b}/\partial a_{i}$ are given by

[TABLE]

Proof:

From (63), to first order in $g$, one has

[TABLE]

Taking the derivative of (65) with respect to $a_{i}$ with the definition of $\bm{V}_{2;\alpha\beta}$ given below (62), and summing over $\alpha>1$ one obtains (64). ∎

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.

Remark 5**.**

Close to the homogeneous case $\bm{M}=m\mathbb{1}$ and $\bm{\Gamma}=\gamma\mathbb{1}$, summing over every fault locations $b$ makes the second term in (64) vanish. One gets $\sum_{b}\alpha_{i}=-g\delta P^{2}\sum_{\alpha>1}u_{\alpha i}^{(0)2}/(2\gamma\lambda_{\alpha}^{(0)})$. This follows from the properties of the eigenvector of $\bm{u}_{\alpha}^{(0)}$.

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 $\mathcal{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.

Theorem 1**.**

For an abrupt power loss, under Assumption 1 and with $\bm{\Gamma}=\gamma\mathbb{1}$, the optimal distribution of parameters $\{r_{i}\}$ that minimizes $\mathcal{M}_{b}$ is obtained as follows.

Compute the sensitivities $\rho_{i}=\partial\mathcal{M}_{b}/\partial r_{i}$ from (44) 2. 2.

Sort the set $\{\rho_{i}\}_{i=1,...N}$ in ascending order 3. 3.

Set $r_{i}=1$ for $i=1,...{\rm Int}[N/2]$ and $r_{i}=-1$ for $i=N-{\rm Int}[N/2]+1,...N$

The optimal placement of inertia and primary control is given by

[TABLE]

The proof is in Appendix A.

Theorem 2**.**

For an abrupt power loss, under Assumption 2 and with $\bm{M}=m\mathbb{1}$, the optimal distribution of parameters $\{a_{i}\}$ that minimizes $\mathcal{M}_{b}$ is obtained as follows.

Compute the sensitivities $\alpha_{i}=\partial\mathcal{M}_{b}/\partial r_{i}$ from (64). 2. 2.

Sort the set $\{\alpha_{i}\}$ in ascending order, 3. 3.

Set $a_{i}=1$ for $i=1,...{\rm Int}[N/2]$ and for $i=N-{\rm Int}[N/2]+1,...N$

The optimal placement of primary control is given by

[TABLE]

Proof:

With Proposition 5 and $\bm{M}=m\mathbb{1}$, we get (64). The proof is the same as the one for Theorem 1 given in Appendix A, but with $\{\alpha_{i}\}$ instead of $\{\rho_{i}\}$ . ∎

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 $\sum_{i}m_{i}=N\,m$, $\sum_{i}d_{i}=N\,d$. From (37), the second condition requires $\sum_{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 $\mathcal{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 $\sum_{i}m_{i}=mN$ and primary control $\sum_{i}d_{i}=dN$ that minimize $\mathcal{M}_{b}$ is obtained as follows.

Compute the parameters $r_{i}$ and $a_{i}$ from Theorems 1 and 2.

1. (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

[TABLE]

Interchange the parameter values $r_{i_{r0}}\leftrightarrow r_{i_{\rm align}}$ and $a_{i_{r0}}\leftrightarrow a_{i_{\rm align}}$. 2. (b)

If $N$ is even, do nothing 2. 2.

If $n\equiv\sum_{i}r_{i}a_{i}=0$, the optimization is done. 3. 3.

Find the set $\mathcal{I}=\{i\,|\,{\rm sgn}(r_{i}a_{i})={\rm sgn}(n)\}$. To reach $\sum_{i}r_{i}a_{i}\rightarrow 0$, our strategy is to set to zero some elements of $\mathcal{I}$. Since however $\sum_{i}a_{i}=\sum_{i}r_{i}=0$ must be conserved, this must be accompanied by a simultaneous change of some other parameter. 4. 4.

Find the pair $(a_{i1},a_{i2}=-a_{i1})$ or $(r_{i1},r_{i2}=-r_{i1})\in\mathcal{I}\times\mathcal{I}$ which, when sent to $(0,0)$, induce the smallest increase of the objective function $\mathcal{M}_{b}$. Send it to $(0,0)$. Because the pair has opposite sign, this does not affect the condition $\sum_{i}a_{i}=\sum_{i}r_{i}=0$. 5. 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

[TABLE]

where the sum runs over all generator buses. The vulnerability measure $\mathcal{V}$ gives a weighted average over all possible fault positions, with the weight $\eta_{b}$ accounting for the probability that a fault occurs at $b$ and $\delta 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, $\eta_{b}\equiv 1$, with Remark 2, it is straightforward to see that $\partial\mathcal{V}/\partial r_{i}=0+\mathcal{O}(\mu^{2})$. Therefore, to leading order, there is no benefit in scaling up the inertia anywhere. On the other hand, with Remark 5, we get $\partial\mathcal{V}/\partial a_{i}=-g\delta P^{2}\sum_{\alpha>1}u_{\alpha i}^{(0)2}/(2\gamma\lambda_{\alpha}^{(0)})+\mathcal{O}(g^{2})$. 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_{\alpha i}^{(0)2}$ of the slow modes of $\bm{L}$ – those with the smallest $\lambda_{\alpha}^{(0)}$. These modes are displayed in Fig 1. One concludes that, with a non-weighted vulnerability measure, $\eta_{b}\equiv 1$ in (68), an homogeneous inertia location is a local optimum for $\mathcal{V}$, for which damping parameters need to be increased primarily on peripheral buses.

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.22$MWs2, and primary control $d_{i}=12.25$MWs. 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 $\mathcal{V}$ of (68) for $\eta_{b}\equiv 1$. This conclusion is confirmed numerically in Fig. 2 (a) and (b). The optimal placement of primary control displayed in panel (b) decreases $\mathcal{V}$ by more than 12% with respect to the homogeneous case.

Setting $\eta_{b}\equiv 1$ in (68) is convenient mathematically, however it treats all faults equally, regardless of their impact. One may instead adapt $\eta_{b}$ to obtain inertia and primary control distributions that reduce the impact of the strongest faults with largest $\mathcal{M}_{b}$. We do this in two different ways, first with $\eta_{b}=\mathcal{M}_{b}^{(0)2}$ and second with

[TABLE]

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 $\eta_{b}\equiv 1$ [Fig. 2 (a) and (b)], we see rather small differences. More importantly, the impact of various choices of $\eta_{b}$ on $\mathcal{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 $\eta_{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 $\eta_{b}$ lead to reductions in $\mathcal{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.

Appendix A

Proof:

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 $\bm{H}$ can be diagonalized as $\bm{t}^{R}\,\bm{\mu}\,\bm{t}^{L}$, where $\bm{\mu}\equiv{\rm diag}(\{\mu_{\alpha s}\})$ and $\bm{t}^{R,L}$ are matrices containing the right and left eigenvectors of $\bm{H}$, the problem is resolved by:

1. changing variables $\bm{\chi}\equiv\bm{t}^{L}[\bm{\xi}^{\top}\dot{\bm{\xi}}^{\top}]^{\top}$ to diagonalize (48), as

[TABLE]

2. solving (70) as

[TABLE]

3. Obtaining $\dot{\xi}_{\alpha}$ via the inverse transformation $[\bm{\xi}^{\top}\bm{\dot{\xi}}^{\top}]^{\top}=\bm{t}^{R}\bm{\chi}$.

These three steps are carried out with the approximate expressions $\bm{t}^{R,L}=\bm{t}^{R,L(0)}+g\bm{t}^{R,L(1)}$ and $\mu_{\alpha\pm}=\mu_{\alpha\pm}^{(0)}+g\mu_{\alpha\pm}^{(1)}$ obtained with the first order in $g$ corrections presented in (60)–(62). One gets

[TABLE]

with

[TABLE]

where

[TABLE]

(63) is obtained from (87) by applying trigonometric identities. ∎

Proof:

To leading order in $\mu=\delta m/m$, this optimization problem is equivalent to the following linear programming problem [30]

[TABLE]

It is solved by the Lagrange multipliers method, with the Lagrangian function

[TABLE]

where $\varepsilon_{i}$ and $\varepsilon_{0}$ are Lagrange multipliers. We get

[TABLE]

The solution must satisfy the Karush-Kuhn-Tucker (KKT) conditions [30], in particular the complementary slackness (CS) condition which imposes that either $\varepsilon_{i}=0$ or $r_{i}=\pm 1\,,\;\forall i$. The former choice leads generally to a contradiction. From (112) and dual feasibility condition, one gets

[TABLE]

This imposes that $r_{i}=-{\rm sgn}(\varepsilon_{0}+\rho_{i})$. To ensure that $\sum_{i}r_{i}=0$ is satisfied, $\varepsilon_{0}$ is set to minus the median value of $\rho_{i}$. If the number of bus $N$ is odd, the $r_{i}$ corresponding to the median value of $\rho_{i}$ is set to zero. ∎