Global analysis of synchronization performance for power systems: bridging the theory-practice gap

TL;DR

This paper extends the analysis of power grid synchronization to realistic scenarios with heterogeneous generators, providing new metrics and insights that bridge theoretical models and practical engineering applications.

Contribution

It introduces a decomposition of frequency responses, defines new metrics for synchronization, and validates these through simulations, connecting theory with real-world power system practice.

Findings

System frequency captures grid behavior and steady-state deviations.

New metrics like Nadir, RoCoF, and synchronization cost quantify dynamic performance.

System frequency remains accurate as network connectivity increases.

Abstract

The issue of synchronization in the power grid is receiving renewed attention, as new energy sources with different dynamics enter the picture. Global metrics have been proposed to evaluate performance and analyzed under highly simplified assumptions. In this paper, we extend this approach to more realistic network scenarios and more closely connect it with metrics used in power engineering practice. In particular, our analysis covers networks with generators of heterogeneous ratings and richer dynamic models of machines. Under a suitable proportionality assumption in the parameters, we show that the step response of bus frequencies can be decomposed in two components. The first component is a {system-wide frequency} that captures the aggregate grid behavior, and the residual component represents the individual bus frequency deviations from the aggregate. Using this decomposition, we define --and compute in closed form-- several metrics that capture dynamic behaviors that are of relevance for power engineers. In particular, using the \emph{system frequency}, we define industry-style metrics (Nadir, RoCoF) that are evaluated through a representative machine. We further use the norm of the residual component to define a \emph{synchronization cost} that can appropriately quantify inter-area oscillations. Finally, we employ robustness analysis tools to evaluate deviations from our proportionality assumption. We show that the system frequency still captures the grid steady-state deviation, and becomes an accurate reduced-order model of the grid as the network connectivity grows. Simulation studies with practically relevant data are included to validate the theory and further illustrate the impact of network structure and parameters on synchronization. Our analysis gives conclusions of practical interest, sometimes challenging the conventional wisdom in the field.

Equations

$$w_i(s)$$

$$p^e(s)$$

$$[L]_{ij}=\frac{\partial}{\partial{}\theta_{j}}{\sum_{j=1}^{n}|V_{i}||V_{j}|b_{ij}\sin(\theta_{i}-\theta_{j})}\Bigr{|}_{\theta=\theta_{0}}.$$

$$m_i\dot{w_i}$$

$$g_i(s)=\frac{1}{m_{i}s+d_i}.$$

$$m_i{\dot{w}}_i$$

$${\tau }_i{\dot{q}}_i$$

$$g_i(s)=\frac{{\tau }_{i}s+1}{m_i{\tau }_i{s}^2+(m_i+d_i{\tau }_i)s+d_i+r_i^{-1}}.$$

$$g_i(s)=\frac{1}{f_i}{g}_0(s).$$

$$G(s)=\mathrm{diag}(g_i(s))=F^{-\frac{1}{2}}[g_0(s)I]F^{-\frac{1}{2}},$$

$$L_F:=F^{-\frac{1}{2}}L{F}^{-\frac{1}{2}},$$

$$L_F=V\Lambda V^T,$$

$$V$$

$$V_{\perp }\in {\mathbb{R}}^{n\times (n-1)},V_{\perp }^T{V}_{\perp }=I_{n-1},V_{\perp }^T{v}_0=0.$$

$$\displaystyle v_{0}=\alpha_{F}F^{\frac{1}{2}}\mathbf{1},\quad\mbox{ with }\alpha_{F}:=\Big{(}\sum_{i}f_{i}\Big{)}^{-\frac{1}{2}}.$$

$$H(s)$$

$$h_k(s)$$

$$T_{wu}(s)=F^{-\frac{1}{2}}VH(s)V^T{F}^{-\frac{1}{2}}.$$

$$\overline{w}(t):=\frac{{\sum }_i{m}_i{w}_i(t)}{{\sum }_i{m}_i},$$

$$w(s)=T_{wu}(s)\frac{1}{s}{u}_0=F^{-\frac{1}{2}}V\frac{H(s)}{s}{V}^T{F}^{-\frac{1}{2}}{u}_0.$$

$${\tilde{h}}_k(s):=\frac{h_k(s)}{s}=\frac{g_0(s)}{s+{\lambda }_k{g}_0(s)},k=0,\dots ,n-1;$$

$$V\frac{H(s)}{s}{V}^T={\tilde{h}}_0(s)v_0{v}_0^T+V_{\perp }\tilde{H}(s)V_{\perp }^T,$$

$$F^{-\frac{1}{2}}{v}_0{v}_0^T{F}^{-\frac{1}{2}}{u}_0={\alpha }_F^{2}11^T{u}_0=\frac{{\sum }_i{u}_{0i}}{{\sum }_i{f}_i}1.$$

$$\displaystyle w(s):=\overbrace{\frac{\sum_{i}u_{0i}}{\sum_{i}f_{i}}\tilde{h}_{0}(s)}^{\mbox{\normalsize$\overline{w}(s)$}}\cdot\mathbf{1}+\overbrace{F^{-\frac{1}{2}}V_{\perp}\tilde{H}(s)V_{\perp}^{T}F^{-\frac{1}{2}}u_{0}}^{\mbox{\normalsize$\tilde{w}(s)$}};$$

$$w(t)=\overline{w}(t)1+\tilde{w}(t),$$

$$(1^{T}F)\tilde{w}(s)\equiv 0.$$

$$\overline{w}(t)=\frac{{\sum }_i{f}_i{w}_i(t)}{{\sum }_i{f}_i};$$

$$\overline{w}(t)=\frac{{\sum }_i{u}_{0i}}{{\sum }_i{f}_i}\cdot {\tilde{h}}_0(t).$$

$$\Vert \tilde{w}{\Vert }_2^2={\int }_0^{\infty }|\tilde{w}(t){|}^{2}dt.$$

$$y_{kl}={\gamma }_{kl}⟨{\tilde{h}}_k,{\tilde{h}}_l⟩$$

$$\text{mbox}where\Gamma =({\gamma }_{kl})$$

$$\tilde{w}(t)=F^{-\frac{1}{2}}{V}_{\perp }\tilde{H}(t)z_0,$$

Full text

Using this decomposition, we define –and compute in closed form– several metrics that capture dynamic behaviors that are of relevance for power engineers. In particular, using the system frequency, we define industry-style metrics (Nadir, RoCoF) that are evaluated through a representative machine. We further use the norm of the residual component to define a synchronization cost that can appropriately quantify inter-area oscillations. Finally, we employ robustness analysis tools to evaluate deviations from our proportionality assumption. We show that the system frequency still captures the grid steady-state deviation, and becomes an accurate reduced-order model of the grid as the network connectivity grows.

Simulation studies with practically relevant data are included to validate the theory and further illustrate the impact of network structure and parameters on synchronization. Our analysis gives conclusions of practical interest, sometimes challenging the conventional wisdom in the field.

I Introduction

The electric power grid, constructed a century ago under the alternating current model, relies strongly on the synchronization of the multiple oscillatory variables at a common nominal frequency (e.g. 60Hz in the U.S., 50Hz in Europe). Automatic control mechanisms are deployed to keep deviations from these set points very small (up to the order of 100mHz). If larger fluctuations are observed, drastic measures are taken: machine protections or load shedding  [1] automatically disconnect network elements, with the potential cost of cascading failures and ultimately blackouts [2].

One challenge for the design of such controls is that frequency measurements are node dependent, and individual nodes see the superposition of two different phenomena: system-wide frequency changes which reflect a global supply-demand imbalance, and inter-area frequency oscillations  [3, 4, 5, 6] due to weak dynamic coupling. In this regard, frequency regulation controllers for normal operation include a fast, decentralized response by local generators, and a slower set-point correction coordinated by the system operator [7]. To identify when protections must be triggered due to a sudden fault, power engineers traditionally rely on step response notions such as the maximum frequency deviation (Nadir) or the maximum rate of change of frequency (RoCoF)[8]; again these are typically based on node specific measurements.

The analysis of grid synchronization from a global perspective has been pursued in academia, going back to eigenvalue methods for reduced order modelling of inter-area oscillations (slow coherency [9], participation factors [10]). More recently, researchers from the control field have proposed to use $\mathcal{H}_{2}$ or $\mathcal{H}_{\infty}$ norms of global transfer functions [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] as figures of merit for synchronization performance. Ideally, these measures should isolate the effect of different aspects of the network structure (topology, damping, inertia, etc.). Of particular importance is to determine whether the decreased inertia in new systems with renewable generation causes, as is feared [21], a deterioration of frequency regulation.

Some model simplifications (mainly, linearization) are naturally required for such analysis. However, other common assumptions (homogeneous machines modeled by swing equations) appear too narrow from the practitioners’ perspective. To bridge the gap between this theory and power engineering practice, we pursue two objectives in this paper. On one hand, to extend the “system norm” approach to the more realistic situation of machines with heterogeneous ratings and where the turbine dynamics is accounted for; secondly, to incorporate step-response metrics (Nadir, RoCoF) for an appropriately defined global system frequency, and to separately characterize inter-area oscillations. Our main contributions are as follows:

We identify a family of heterogeneous machines for which the previous objectives are attainable exactly in closed form. This occurs when dynamic parameters are proportional to machine rating, a mild restriction compared with homogeneity. In this case, we show how to decompose the system step response into a system frequency $\overline{w}(t)$ (motion of the center of inertia) and a transient component of inter-area oscillations. 2. 2.

In the above coordinates, we provide natural definitions of performance metrics relevant to power engineers: Nadir and RoCoF are defined as the $\mathcal{L}_{\infty}$ norms of, respectively, $\overline{w}(t)$ and $\dot{\overline{w}}(t)$. A synchronization cost measuring oscillations is defined as the $\mathcal{L}_{2}$ norm of the vector of deviations. We give a general result on how the latter cost is computed, separating the topological and dynamic components. 3. 3.

We apply these methods to different machine models, in particular beyond the traditional swing equation model to include turbine control. The theory helps reveal the impact on relevant parameters; in particular we find the rather limited role of system inertia on synchronization performance, as compared to damping. 4. 4.

We analyze the situation where our proportionality assumption does not hold, invoking robust control tools. Although in this case it is not possible to perfectly decouple system-wide and inter-area dynamics, we provide bounds on the deviation from the idealized case. We show in particular how the desired decomposition emerges in the limit of high network connectivity. 5. 5.

We provide simulation studies to validate the approach, based on models of the Icelandic power grid [22]. We first modify the machine data to impose paramter proportionality, but respecting the heterogeneous ratings. Our results provide insightful conclusions on the role of the different parameters. Finally we study the real (non-proportional) network, and also explore the influence of connectivity.

The rest of the paper is organized as follows. Section II contains some background material on the network model. In Section III we introduce our proportionality assumption, carry out the relevant decomposition and introduce the system norms, expressed in terms of a representative machine and the network structure. In Sections IV and V we apply the results to different machine models, respectively to the second-order swing dynamics, and a third-order model that incorporates the turbine control. Section VI addresses the case where the proportionality assumption is removed. In Section VII we present our simulation studies, and conclusions are presented in Section VIII. Some of the more extensive derivations are deferred to the appendices. Partial versions of the results in this paper were presented in [23].

II Preliminaries

We consider a set of $n$ generator buses, indexed by $i\in\{1,\dots,n\}$, dynamically coupled through an AC network. Assuming operation around an equilibrium, the linearized dynamics are represented by the block diagram in Fig. 1.

$G(s)=\mathrm{diag}(g_{i}(s))$ is the diagonal transfer function of generators –with local controllers in closed loop– at each bus. Each $g_{i}(s)$ has as output the angular velocity $w_{i}=\dot{\theta}_{i}$, and as input the net power at its generator axis, relative to its equilibrium value. This includes an outside disturbance $u_{i}$, reflecting variations in mechanical power or local load, minus the variation $p^{e}_{i}$ in electrical power drawn from the network:

[TABLE]

The network power fluctuations $p^{e}$ are given by a linearized model of the power flow equations:

[TABLE]

where $L$ is a undirected weighted Laplacian matrix of the network with elements

[TABLE]

Here $\theta_{0}$ are the equilibrium angles, $|V_{i}|$ is the (constant) voltage magnitude at bus $i$, and $b_{ij}$ is the line susceptance.

Remark 1 (Model Assumptions).

The linearized network model (1)-(2) implicitly makes the following assumptions which are standard and well-justified for frequency control on transmission networks [24]:

• •

Bus voltage magnitudes $|V_{i}|$ are constant; we are not modeling the dynamics of exciters used for voltage control; these are assumed to operate at a much faster time-scale.

• •

Lines $ij$ are lossless.

• •

Reactive power flows do not affect bus voltage phase angles and frequencies.

• •

Without loss of generality, the equilibrium angle difference ($\theta_{0,i}-\theta_{0,j}$) accross each line is less than $\frac{\pi}{2}$.

Also implicit in our model is the so-called Kron reduction, a procedure by which algebraic constraints resulting from buses with no generators are eliminated; see e.g. [25, 26]. For a first principle derivation of the model we refer to [27, Section VII]. For applications of similar models for frequency control within the control literature, see, e.g., [28, 29, 30].

Two examples of generator dynamics, to be considered explicitly in this paper, are:

Example 1.

The swing equation dynamics

[TABLE]

which takes the form of the rotational Newton’s law, with an inertia $m_{i}$, and a damping $d_{i}$; forcing is the power imbalance. The typical swing equations (see e.g. [26]) are of second order in the machine angles; here we are not making these explicit so the model is of first order.

The corresponding transfer function of Fig. 1 is

[TABLE]

Example 2.

The swing equation with a first-order model of turbine control:

[TABLE]

Here $q_{i}$ is the (variation of) turbine power, $\tau_{i}$ the turbine time constant and $r_{i}$ the droop coefficient; the governor dynamics are considered to be faster and neglected. The corresponding transfer function is

[TABLE]

Of course, other models are possible within this framework. We make the following assumption:

Assumption 1.

The transfer function $g_{i}(s)$ is stable (has all its poles in $Re(s)<0$) and strictly positive real, i.e. $Re[g_{i}(j\omega)]>0$ for all $\omega\in{\mathbb{R}}$.

This assumption implies that stability of the system, whichever the network.

Proposition 1.

The feedback loop of Figure 1 is internally stable.

Proof.

The transfer function $\frac{1}{s}L$ is always positive real (see [31]) since the Laplacian is a positive semidefinite matrix. Invoking the passivity theorem (see e.g. [31]) implies the result. ∎

It is straightforward that the swing model (3) satisfies Assumption 1. Checking it for the model (4) requires some calculations, or to note that this model is obtained by feedback of (3) with a passive first-order controller.

For more relaxed conditions for internal stability that are robust to network uncertanties we refer the reader to [32].

III Modal decomposition for proportionally heterogeneous machines

A popular research topic in recent years [11, 12, 13, 14, 15, 16, 17] has been the application of global metrics from robust control to this kind of synchronization dynamics, as a tool to shed light on the role of various parameters, e.g. system inertia. Most of the analytical results, however, consider a homogeneous network where all machines are identical (i.e., common $m_{i}$, $d_{i}$, etc.), a very restrictive scenario.111Some bounds on heterogenous systems are given in [11, 14]. Numerical studies with heterogeneity are given in [12].

In a real network, where generators have different power ratings, it is natural for parameters to scale accordingly: for instance, the inertia $m_{i}$ of a machine will grow with its rating, and it is clear that “heavier” machines will have a more significant impact in the overall dynamics.

While in principle one would like to cover general parameters, we will develop most of our theory for parameters satisfying a certain proportionality.

Assumption 2 (Proportionality).

There exists a fixed transfer function $g_{0}(s)$, termed the representative machine, and a rating parameter $f_{i}>0$ for each bus $i$, such that

[TABLE]

To interpret this, consider first the swing dynamics of Example 1. Then the assumption is satisfied provided that inertia and damping are both proportional to $f_{i}$, i.e. $m_{i}=f_{i}m$, $d_{i}=f_{i}d$, where $m$, $d$ are those of a representative machine. Equivalently, the ratios $m_{i}/d_{i}$ are uniform over $i$; this kind of proportionality is termed “uniform damping” in [25].

Going to the case of Example 2 with the turbine dynamics, we find that the assumption is satisfied provided that $m_{i}=f_{i}m$, $d_{i}=f_{i}d$, $r_{i}^{-1}=f_{i}r^{-1}$, $\tau_{i}=\tau$; here the inverse droop coefficient is assumed proportional to rating, but the turbine time-constant is taken to be homogeneous.

Regarding the practical relevance of our simplifying assumption: empirical values reported in [33] indicate that at least in regard to orders of magnitude, proportionality is a reasonable first-cut approximation to heterogeneity, substantially more realistic than the homogeneous counterpart. In Section VI we will discuss deviations from this assumption.

III-A Diagonalization

We will now exploit Assumption 2 to decompose the dynamics of Fig. 1 in a manner that allows for a suitable decoupling in the analysis. In what follows, $F=\mathrm{diag}(f_{i})$ denotes the diagonal matrix of rating parameters. Writing

[TABLE]

we transform the feedback loop into the equivalent form of Fig. 2.

We introduce a notation for the scaled Laplacian matrix222This scaling already appears in the classical paper [9] and, more recently, in [34, 35].

[TABLE]

which is positive semidefinite and of rank $n-1$. Applying the spectral theorem we diagonalize it as

[TABLE]

where $\Lambda=\mathrm{diag}(\lambda_{k}),\quad 0=\lambda_{0}<\lambda_{1}\leq\cdots\leq\lambda_{n-1},$ and $V$ is unitary. Distinguishing the eigenvector $v_{0}$ that corresponds to the zero eigenvalue, we write

[TABLE]

In fact $v_{0}$ can be made explicit by recalling that $\ker(L)=\mathrm{span}\{\mathbf{1}\}$, so $\ker(L_{F})=\mathrm{span}\{F^{\frac{1}{2}}\mathbf{1}\}$, from where

[TABLE]

Substitution of (6) into Fig. 2 and some block manipulations leads to the equivalent representation of Fig. 3.

Noting finally that the $V$ block commutes with $g_{0}(s)I$ and thus cancels out with $V^{T}$, the internal loop is now fully diagonalized, yielding the closed-loop transfer function

[TABLE]

Note that by virtue of Assumption 1, $H(s)$ is a stable transfer function, as a consequence of the passivity of the integrator feedback. The overall transfer function between the vector of external power disturbances and the machine frequency outputs is

[TABLE]

III-B Step response decomposition

Global metrics for synchronization performance in e.g. [11, 14, 15], are system norms ($\mathcal{H}_{2}$, $\mathcal{H}_{\infty}$) applied to $T_{wu}$ (frequency output), or to a “phase coherency” output based on differences in output angles. The choice of metric carries an implicit assumption on the power disturbances considered (white noise, or a worst-case $\mathcal{L}_{2}$ signal).

In this paper, we wish to analyze metrics which more closely reflect industry practice. The standard [36] provides a detailed description of the Frequency Response to a step change in generation or load, and the associated time-domain performance metrics. The steady-state values of the response are part of the Balancing Authority’s requirements; other transient features (Nadir, RoCoF, see below) are also mentioned as important.

Now, since buses do not respond in unison, some globally representative system frequency must be considered. A candidate is the weighted average

[TABLE]

which corresponds to the motion of the center of inertia (COI), a classical notion [24, 25].

We now show that for our family of heterogeneous systems, the behavior of $\overline{w}(t)$ decouples nicely from the individual bus deviations $w_{i}(t)-\overline{w}(t)$, opening the door for a separate analysis of both aspects of the problem.

Remark 2.

We are not considering other disturbances such as line outages which alter the structure of the network itself, in particular the Laplacian matrix. An extreme case would be a system split, which would lead to different COIs in each portion of the disturbed network. These alternatives are a topic for future research.

Our input will be a step function $u(t)=u_{0}\mathds{1}_{t\geq 0}$; here $u_{0}$ is a given vector direction. In Laplace transforms:

[TABLE]

Denote by the diagonal elements of $\frac{H(s)}{s}$ by

[TABLE]

we note they are all stable except for $\tilde{h}_{0}(s)=g_{0}(s)/s$. This suggests the following decomposition:

[TABLE]

where $\tilde{H}(s)=\mathrm{diag}\left(\tilde{h}_{1}(s),\ldots,\tilde{h}_{n-1}(s)\right)$, and we have invoked the submatrices of $V$ from (7). Furthermore, recalling from (8) that $v_{0}=\alpha_{F}F^{\frac{1}{2}}\mathbf{1}$, we have

[TABLE]

Substitution into (12) gives the decomposition

[TABLE]

or in the time domain,

[TABLE]

interpreted as follows:

• •

$\overline{w}(t)$ is a system frequency term, applied to all nodes;

• •

the transient term $\tilde{w}(t)$ represents the individual node deviations from the synchronous response.

III-C System frequency

We can obtain more information on the system frequency by observing that since $\mathbf{1}^{T}F^{\frac{1}{2}}V_{\perp}=\alpha_{F}^{-1}v_{0}^{T}V_{\perp}=0$, we have

[TABLE]

Therefore $\mathbf{1}^{T}Fw(t)=\overline{w}(t)\mathbf{1}^{T}F\mathbf{1}$ by (15), which gives

[TABLE]

the system frequency is a weighted mean of bus frequencies, in proportion to their rating. Noting that $m_{i}=mf_{i}$, it follows that $\overline{w}(t)$ is exactly the COI frequency from (11).

Also, returning to (14) we have

[TABLE]

Since $\tilde{h}_{0}(s)=g_{0}(s)/s$, then $\tilde{h}_{0}(t)$ is the step response of the representative machine. Thus $\overline{w}(t)$ corresponds to the angular frequency observed when exciting the representative machine (in open loop) with the total system disturbance normalized by the total scale.

Remark 3.

Note that this result is independent of $L$, i.e. the electrical network does not affect the time response of the system frequency, only the machine ratings themselves. Thus, when the network dependent term ($\tilde{w}$) converges fast to zero, (16) is a natural candidate for a reduced order model similar to the ones recently considered in [37, 38].

In the following sections we analyze its behavior for the previously discussed examples.

III-D Quantifying the deviation from synchrony

We now turn our attention to the term $\tilde{w}(t)$ which represents individual bus deviations from a synchronous response. A natural way of quantifying the size of this transient term is through the $\mathcal{L}_{2}$ norm

[TABLE]

We now show how this norm can be computed in terms of the parameters of the scaled network Laplacian, and the impulse response matrix $\tilde{H}(t)=\mathrm{diag}\left(\tilde{h}_{1}(t),\ldots,\tilde{h}_{n-1}(t)\right)$, Laplace inverse of $\tilde{H}(s)$, which encapsulates all information on the machine model.

Proposition 2.

$\|\tilde{w}\|_{2}^{2}=z_{0}^{T}Yz_{0}$, where:

• •

${Y}\in{\mathbb{R}}^{(n-1)\times(n-1)}$* is the matrix with elements*

[TABLE]

• •

$z_{0}:=V_{\perp}^{T}F^{-\frac{1}{2}}u_{0}\in{\mathbb{R}}^{n-1}$.

Proof:

With the introduced notation we have

[TABLE]

therefore $\tilde{w}(t)^{T}\tilde{w}(t)=z_{0}^{T}\tilde{H}(t)\Gamma\tilde{H}(t)z_{0}.$ The matrix in the above quadratic form has elements $\tilde{h}_{k}(t)\gamma_{kl}\tilde{h}_{l}(t)$, therefore integration in time yields the result. ∎

Remark 4.

The metric $\|\tilde{w}\|_{2}^{2}$ does depend on the electrical network, through the eigenvalues and eigenvectors of $L_{F}$.

III-E *Mean synchronization cost for random disturbance

step*

Since the cost discussed above is a function of the disturbance step $u_{0}$, it may be useful to find its average over a random choice of this excitation. Recalling that the components $u_{0i}$ correspond to different buses, it is natural to assume them to be independent, and thus $E[u_{0}u_{0}^{T}]=\Sigma^{u}$, a diagonal matrix.

Using the fact that $z_{0}=V_{\perp}^{T}F^{-\frac{1}{2}}u_{0}$ we get

[TABLE]

and the expectation for the cost in Proposition 2 is

[TABLE]

Therefore, given $\Sigma^{u}$ we define the mean synchronization cost as

[TABLE]

We look at some special cases:

• •

$\Sigma^{u}=I$ (uniform disturbances). Then $\Sigma^{z}=\Gamma$, and

[TABLE]

• •

$\Sigma^{u}=F$. This means disturbance size follows the square root of the bus rating. Here $\Sigma^{z}=I$, and

[TABLE]

• •

$\Sigma^{u}=F^{2}$. This is probably most natural, with disturbances proportional to bus rating. Here $\Sigma^{z}=V_{\perp}^{T}FV_{\perp}=\Gamma^{\dagger}$ (pseudoinverse); $E\left[\|\tilde{w}\|_{2}^{2}\right]=\mbox{\rm Tr}(Y\Gamma^{\dagger})$.

III-F The homogeneous case

If all machines have the same response $g_{0}(s)$, setting $F=I$ we can obtain some simplifications:

• •

$\{\lambda_{k}\}$ are the eigenvalues of the original Laplacian $L$.

• •

The system frequency is the average $\overline{w}(t)=\frac{1}{n}\sum_{i}w_{i}(t)$, and satisfies

[TABLE]

• •

$z_{0}=V_{\perp}^{T}u_{0}$, and $\Gamma=V_{\perp}^{T}V_{\perp}=I$. Therefore the matrix $Y$ in Proposition 2 is diagonal, $Y=\mathrm{diag}(\|\tilde{h}_{k}\|^{2})$, and

[TABLE]

where $z_{0k}=v_{k}^{T}u_{0}$ is the projection of the excitation vector $u_{0}$ in the direction of the $k$-th eigenvector of the Laplacian $L$.

• •

The mean synchronization cost for $E[u_{0}u_{0}^{T}]=I$ (here all the preceding cases coincide) is

[TABLE]

the $\mathcal{H}_{2}$ norm of the transfer function $\tilde{H}(s)$. We recall that this was obtained by isolating the portion $\tilde{h}_{0}(s)$ corresponding to the synchronized response (in this case, projecting onto $\mathbf{1}^{\perp}$). In this form, the cost resembles other proposals [11, 15], for the price of synchrony, and [16] for the evaluation of the synchronization cost under step changes in homogeneous systems.

IV Application to the swing dynamics

In this section we assume we are in the situation of Example 1, i.e., the representative machine is

[TABLE]

and $m_{i}=f_{i}m$, $d_{i}=f_{i}d$ are the individual bus parameters. The diagonal closed loop transfer functions in (III-A) are

[TABLE]

and the corresponding step response elements are

[TABLE]

IV-A System frequency

Inverting the transform for the case $\lambda_{0}=0$ we have

[TABLE]

and therefore invoking (16) we find that

[TABLE]

We remark the following:

• •

Again, $\overline{w}(t)$ does not depend on the electrical network.

• •

The first-order evolution of $\overline{w}(t)$ implies there is no overshoot; system frequency never deviates to a “Nadir” further from equilibrium than its steady-state value.

• •

The asymptotic frequency is the ratio $w_{\infty}=\frac{\sum_{i}u_{0i}}{\sum_{i}d_{i}}$ of total disturbance to total damping. As it is well known [24], it does not depend on the inertia $m$, which only affects the time constant in which this asymptote is achieved.

• •

The maximum RoCoF (rate-of-change-of frequency) occurs at $t\to 0+$, and is given by

[TABLE]

here the total inertia appears, which is natural in the response to a step in forcing. RoCoF increases for low inertia, however it need not have a detrimental impact: system frequency initially varies quickly but never deviates more than $w_{\infty}$, independent of $m$.

IV-B Synchronization cost

The synchronization cost $\|\tilde{w}\|_{2}$ for this case, can be computed by particularizing the result in Proposition 2. The following result is proved in Appendix A.

Proposition 3.

Let $\tilde{h}_{k}(s)$ be given in (23), and $\tilde{h}_{k}(t)$ its inverse transform, for $k=1,\ldots,n-1$. Then:

[TABLE]

It follows that the matrix $Y$ in (17) will depend on both inertia $m$ and damping $d$, so in general both have an impact on the “price of synchrony”. Note however that inertia only appears in off-diagonal terms, and the matrix remains bounded as $m\to 0$ or $m\to\infty$; we thus argue that inertia has limited impact. Let us look at this issue in more detail.

IV-B1 Homogeneous case

In the case of homogeneous machines, we saw above that $\Gamma=I$ and $Y$ is diagonal, so inertia disappears completely: indeed using (20) we have

[TABLE]

The cost is inversely proportional to damping, and the direction of the disturbance $u_{0}$ also matters. Recalling that $v_{k}$ is the $k$-th Laplacian eigenvector, the worst-case for a given magnitude $|u_{0}|$ is when it is aligned to $v_{1}$, the Fiedler eigenvector.

If the disturbance direction is chosen randomly as in Section III-F, then (21) gives

[TABLE]

again a similar result to those in [11] for homogeneous systems.

IV-B2 Heterogeneous, high inertia case

Assume for this discussion that all the $\lambda_{k}$ are distinct; then as $m\to\infty$ we have $y_{kl}\to 0$ for $k\neq l$, so $Y$ again becomes diagonal, and the cost has the limiting expression

[TABLE]

So the high inertia behavior is of a similar structure to the homogeneous case in (28). Comparisons are not straightforward, though, since the scaling factor $F$ affects $z_{0k},\gamma_{kk}$ and $\lambda_{k}$ in each of the above terms.

IV-B3 Heterogeneous, low inertia case

If $m\to 0$, then the limiting $Y$ matrix is not diagonal. The corresponding limiting cost is

[TABLE]

Note, however, that the diagonal terms are the same as in the high inertia case. This suggests that inertia plays a limited role in the $\mathcal{L}_{2}$ price of synchrony, even in the heterogeneous machine case. The simulations of Section VII are consistent with this observation.

Reviewing all results for the swing equation model, one may object that it is perhaps not expressive enough to capture the behaviors of interest. In particular, the lack of an overshoot in the COI step response departs from practical observations [36]; one of the main missing features is the presence of a governor with delayed response. To address this objection we consider the more detailed model of the following section.

V Model with turbine dynamics

In this section we use the model of Example 2, where the representative machine is

[TABLE]

The transfer functions for the step response in (13) are:

[TABLE]

It can be checked (e.g. by applying the Routh-Hurwitz criterion) that they are stable whenever $\lambda_{k}>0$.

V-A System frequency

We can again use (16) and (32) to compute the system frequency

[TABLE]

but now finding the inverse transform of $\tilde{h}_{0}(s)$ is more involved. Using partial fractions we first express

[TABLE]

The first term provides the steady-state response, which is

[TABLE]

this is analogous to the swing equation case, except than the droop control has been added to the damping; the regulation standards [36] apply to this quantity. Again, as expected, inertia plays no role at all in this steady-state deviation.

The transient term is a second-order transfer function, which we proceed to analyze now. Its behavior critically depends on whether its poles are real or complex conjugate. In particular, whenever

[TABLE]

the system is under-damped with poles $\eta\pm j\omega_{d}$, and

[TABLE]

where

[TABLE]

The system frequency time evolution is thus given by

[TABLE]

A few observations are in order:

• •

Including the turbine model has a nontrivial effect on the system frequency $\bar{w}(t)$. It is the presence of the turbine dynamics that provides the characteristic under-damped behavior that produces a Nadir.

• •

We have only provided here the solution of $\bar{w}(t)$ for the (practically more relevant) under-damped case.

• •

Interestingly, (35) shows that the system may become over-damped by either increasing $m$, or decreasing $m$! However, the behavior is different for each case: in the very high inertia case the Nadir disappears; whereas when $m$ goes to zero, there is an overshoot in the overdamped response. Since in practice this occurs only for very low inertia and already way beyond the acceptable deviation, we are justified in our focus on the under-damped case.

We now proceed to compute the Nadir and RoCoF for this situation. The proofs of the following propositions are found in Appendix B.

Proposition 4 (Nadir).

Given a power system under Assumption 2 with generators containing first order turbine dynamics ($g_{i}(s)$ given by (4)). Then under the under-damped condition (35), the Nadir of the system frequency $\overline{w}(t)$ is given by

[TABLE]

where the phase $\phi\in(-\frac{\pi}{2},\frac{\pi}{2})$ is uniquely determined by

[TABLE]

The dependence of (39) on $m$ is not straightforward, as $\phi$, $\eta$, and $\omega_{d}$ depend on it. The next proposition shows that the dependence is as expected by conventional power engineering wisdom.

Proposition 5.

Given a power system under Assumption 2 with generators containing first order turbine dynamics ($g_{i}(s)$ given by (4)). Then under the under-damped condition (35), the maximum frequency deviation $||\bar{w}||_{\infty}$ is a decreasing function of $m$, i.e., $\frac{\partial}{\partial m}||\overline{w}||_{\infty}<0$.

We now turn to the maximum rate of change of frequency.

Proposition 6 (RoCoF).

Given a power system under Assumption 2 with generators containing first order turbine dynamics ($g_{i}(s)$ given by (4)). Then under the under-damped condition (35), the RoCoF is given by

[TABLE]

The main difference with Proposition 4 is that, while not trivial to establish, here the maximum is always achieved at $t=0+$, exactly as in the second order case of (26).

Again we find, as expected, that the RoCoF decreases with $m$. But this initial slope need not be consequential.333 If there are system protections activated by RoCoF, these may need to be recalibrated to operate with low inertia, see [21].

V-B Synchronization cost

The synchronization cost $\|\tilde{w}\|^{2}$ can once again be computed through Proposition 2, which requires finding the inner products $\langle{\tilde{h}_{k}},{\tilde{h}_{l}}\rangle$, in this case for the functions in (33).

Since the corresponding expression is in general rather unwieldy (see Appendix A), we will present some simpler cases, beginning with $k=l$; the norm from Appendix A is:

[TABLE]

V-B1 Homogeneous case

The above expression suffices to analyze the case of homogeneous machines, where $\Gamma=I$ and $Y$ is diagonal. We have from (20) that

[TABLE]

from (42) we see that, in contrast to the second order machine model, the inertia $m$ does affect the synchronization cost. A closer look at $\|\tilde{h}_{k}\|^{2}$ as a (linear fractional) function of $m$ shows that it is decreasing in $m\in(0,\infty)$, going from

[TABLE]

to

[TABLE]

So higher inertia is beneficial here. Recalling that the corresponding cost for the swing dynamics is $\frac{1}{2\lambda_{k}d}$, we see that it has been reduced. In the high inertia case, the main change is the increased damping through the droop coefficient $r^{-1}$.

V-B2 Heterogeneous, high inertia case

As mentioned, the formula for $\langle{\tilde{h}_{k}},{\tilde{h}_{l}}\rangle$ for $k\neq l$ is quite formidable, but we can give its approximation in the limit of large $m$:

[TABLE]

This assumes $\lambda_{k}\neq\lambda_{l}$. So if the eigenvalues of the scaled Laplacian $L_{F}$ are distinct, we see that again the matrix $Y$ becomes diagonal as $m\to\infty$. The limiting cost is

[TABLE]

This expression amounts to reducing to the cost (30) for the swing dynamics, by the fraction $\frac{d}{r^{-1}+d}$. So the role of the turbine in a high inertia system is again mainly a change in the droop coefficient.

V-B3 Heterogeneous, low inertia case

In the low inertia limit, we find that $\langle{\tilde{h}_{k}},{\tilde{h}_{l}}\rangle\stackrel{{\scriptstyle m\to 0}}{{\longrightarrow}}\frac{N}{D}$, where

[TABLE]

So the limiting matrix $Y$ is not diagonal, as in the second order case; an expression analogous to (31) can be written. Comparisons between the two are not straightforward here, and must be pursued by numerical experimentation; some of these are reported in Section VII.

Despite the more intricate formulas that preclude from simple conclusions, one observation can be made: system inertia still does not play a major role in the synchronization cost. In particular it does not explode in either end of the inertia spectrum, reflecting a certain degree of indifference; we will give further comments after our experimental results.

VI Beyond the proportionality assumption

The analysis of the preceding sections was predicated on the proportionality of Assumption 2. Without it, it is not possible to isolate the center of inertia response from the remaining oscillatory fluctuations. Nevertheless, some conclusions can be drawn by analyzing perturbations from the proportionality assumption. We start with the swing equation case.

VI-A Non-proportionality in swing dynamics

We consider the situation of Example 1, but where the mass to damping ratios are not uniform across buses. Choosing the inertia as indication of scale, we will continue to write $m_{i}=mf_{i}$ for a representative mass $m$. Since the ratios $d_{i}/f_{i}$ are not uniform, we will select as representative damping

[TABLE]

and define for each bus the perturbation parameter

[TABLE]

A key property of the chosen definition is that

[TABLE]

We now find an expression for the transfer function $g_{i}(s)$, in terms of the representative machine $g_{0}(s)$ in (22) and the perturbation:

[TABLE]

This can be verified through standard calculations, and is equivalent to the feedback interconnection of Fig. 5.

Including such feedback uncertainty in each machine $g_{i}(s)$ in Fig. 1, and commuting it with the network feedback, leads to the diagram of Fig 6. Here $\Delta=\mathrm{diag}(\delta_{i})$ is a diagonal perturbation.

In this commuted representation, the inner loop corresponds exactly to the proportional situation analyzed in Section III, yielding the transfer function $T_{wu}(s)$ in (10). The overall perturbed transfer function is

[TABLE]

Replacing with (10) leads after some manipulations to

[TABLE]

for $H(s)$ in (III-A) and $\tilde{\Delta}:=V^{T}\Delta V$.

We observe that the above expression reverts back to the one in (10) if $\Delta=0$ (nominal case with proportionality). This suggests a robustness analysis for the perturbed performance. A first comment is that robust stability is not an issue: since our non-proportional plant model still satisfies Assumption 1, the passivity argument for closed loop stability still holds.

Now when considering the step response, the elegant decoupling of Section III breaks down, due to the term $(I-H(s)\tilde{\Delta})^{-1}$ which does not have diagonal structure, since $\tilde{\Delta}$ is not diagonal. There is no natural way to isolate a “system frequency”; one may still define it by the motion of the COI, but its dynamics is of high order depending on all synchronization modes.

Some partial analysis is nevertheless possible.

VI-A1 Steady-state step response

This quantity can be found through the simple structure of $H(s)$ at $s=0$, namely

[TABLE]

where $e_{1}$ is the first coordinate vector.

Lemma 7.

$H(0)\tilde{\Delta}H(0)=0$, and consequently:

• (i)

$(I-H(0)\tilde{\Delta})^{-1}=I+H(0)\tilde{\Delta}$,

• (ii)

$(I-H(0)\tilde{\Delta})^{-1}H(0)=H(0)$.

Proof.

For the first identity it suffices to note that

[TABLE]

invoking (44). Then (i) and (ii) follow directly. ∎

The preceding lemma implies that the nominal and perturbed systems coincide at $s=0$:

[TABLE]

Here we have used that $g_{0}(0)=1/d$ for the swing model.

Under a step perturbation of direction $u_{0}$, the residual at $s=0$ of $T^{\Delta}_{wu}(s)\frac{1}{s}u_{0}$ is

[TABLE]

so all buses converge to the same steady-state frequency, with the same value as the one in (25); but here we cannot claim the same transient behavior for the center of inertia frequency.

VI-A2 Behavior as connectivity grows

While one cannot give a clean expression for the transfer function $T^{\Delta}_{wu}(s)$ outside $s=0$, it is worth considering what happens when the network connectivity grows, reflected in the growth of the nonzero eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{n-1}$. Indeed, for any fixed $s$ in the right half-plane $Re(s)\geq 0$, we have

[TABLE]

Therefore we can emulate the argument above to conclude that

[TABLE]

This amounts to saying that all buses respond in the same way as the representative machine, appropriately scaled, excited by the aggregate of all bus disturbances. The use of $g_{0}(s)$ as the natural reduced-order model of the system is thus formally justified (recall Remark 3).

The limit can also be extended from a single point $s$ to a compact region of the plane444Since the transfer functions $h_{k}(s)$ have a resonant peak at $s=j{\sqrt{\lambda_{k}}}/{m}$ with value $\frac{1}{d}$, the limit in (46) cannot be taken uniformly over right-half plane, i.e. convergence in $\mathcal{H}_{\infty}$.; for instance one can claim that the Fourier response $T_{wu}^{\Delta}(j\omega)$ for $|\omega|\leq B$ converges uniformly as connectivity grows to a response of the form (47). So the low bandwidth behavior of the generators becomes highly synchronized, irrespective of non-proportionality. We will revisit this issue in our simulation section.

VI-B Non-proportionality in the turbine model

As mentioned in Section II, the turbine model can be viewed as the negative feedback interconnection of the swing dynamics with a controller $K_{i}(s)$, that adds a power injection as a function of measured bus frequency.

Our development was for the first order controller $K_{i}(s)=\frac{r_{i}^{-1}}{\tau_{i}s+1}$. Proportionality in this context required, in addition to the proportional natural damping $d_{i}=df_{i}$, the use of a proportional droop coefficient $r_{i}^{-1}=r^{-1}f_{i}$ and a uniform $\tau_{i}=\tau$. In particular in this case we have $K_{i}(s)=f_{i}K_{0}(s)$, where $K_{0}(s)=\frac{r^{-1}}{\tau s+1}$ is a fixed transfer function.

Non-proportionality is now modeled as a dynamic perturbation of the nominal, normalized controller $K_{0}(s)$:

[TABLE]

Choosing its gain appropriately one can always ensure that

[TABLE]

is zero at $s=0$. In this way, the steady-state analysis carried out in the swing model extends.

In particular, the feedback structure is analogous to that in Figure 5, where the static $\delta_{i}$ is replaced by $\delta_{i}+\Delta^{K}_{i}(s)$; with this change, the rest of the analysis can be carried through, to an expression (45) where now $\tilde{\Delta}(s)=V^{T}(\Delta+\Delta^{K}(s))V$.

Focusing on the behavior at $s=0$, the perturbation will still satisfy the conditions of Lemma 7; therefore the steady-state response is again unchanged from the proportional case. Similar statements can be made about the low-bandwidth behavior as the connectivity grows.

VII Numerical Illustrations

In this section we present numerical simulations to validate our analysis, as well as to provide insight on how generator and network parameters affect the synchronization performance.

Our starting point is a set of real network data available for the Icelandic power grid [22]. Since the parameters $m_{i}$, $d_{i}$, $r_{i}$, $\tau_{i}$ in this data do not satisfy the proportionality of Assumption 2, to illustrate the theory we first produce a set of synthetic parameters as follows. First, define $f_{i}$ and the representative mass using

[TABLE]

then define the remaining parameters of the representative machine:

[TABLE]

In this way we obtain machines with the same heterogeneity in ratings as the real ones, but satisfying Assumption 2; this network is simulated in subsections VII-A and VII-B. Finally in subsection VII-C we will simulate the real network data to assess the effect of the lack of proportionality.

We use the per unit system, standard in the power engineering community, to run our simulations and present our numerical results. We refer to [24, 39, 40] for more details.

VII-A Swing Dynamics

We first illustrate the behavior of a network obeying the swing dynamics and satisfying the proportionality assumption.

Figure 7 corresponds to the parameters obtained from using (48) and (49) on the Icelandic grid data [22]. It shows the response of the bus frequency $w$, the system frequency $\bar{w}$ and the synchronization error $\tilde{w}$, after a disturbance of $-3$p.u. is introduced at bus number 2. A synchronization cost $||\tilde{w}||_{2}=4.77$ is incurred.

The system frequency evolution $\overline{w}(t)$ is consistent with a first order response, as predicted by equation (25). In particular it presents no overshoot, no aggregate oscillatory behavior.

Synchronization performance is reflected in the norm of deviations $\tilde{w}$ from system frequency. To study its dependence with parameters we use $||\tilde{w}||_{2,\Sigma_{u}}$ (averaging in the disturbance direction with $\Sigma_{u}=F^{2}$). Figure 8 presents values for this metric as a function of $m$ and $d$, with fixed heterogeneity scalings $f_{i}$. The red dot represents the parameters used in Figure 7. It can be seen that while changes on the inertia $m$ have nearly no effect on the synchronization cost, changes on the damping $d$ can significantly affect it.

This conclusion is corroborated by the time simulations in Fig. 9, which correspond to increasing by a factor of 2.5 either the system damping $d$ or the inertia $m$. As expected, increasing $d$ not only reduces the synchronization cost $||\tilde{w}||_{2}$ (from 4.77 to 2.96), but also renders a smaller steady-state frequency deviation. Increasing inertia has essentially no effect on our cost, and none on the steady state; it is mainly reflected in a slower rate of convergence.

VII-B System with Turbine Dynamics

We now turn to simulations with the model of Section V, which includes the lagged response of the turbine control.

Figure 10 shows the effect of a $-3$p.u. step change on bus 2. When compared with Figure 7, the system frequency $\overline{w}$ now is seen to have the characteristic second order response, including an overshoot. The steady-state frequency is closer to the nominal due to the action of the turbine droop ($r^{-1}$).

As for the synchronization cost, it is noteworthy that it is essentially unchanged with respect to the swing model. Indeed, if we explore in Fig. 11 the effect on this cost of varying parameters $m$ and $d$, we find very similar behavior to Fig. 8.

We supplement this in Fig. 12 with a study of the dependence of the mean synchronization cost on parameters $(r,\tau)$; again these are not very significant around the nominal (red) point. In fact we only see an impact for very low values of $\tau$, which would correspond to very fast droop control, essentially equivalent to increasing the damping $d$. This prescription for an improved synchronization is consistent with results in [41].

The dependence of the Nadir $||\bar{w}||_{\infty}$ on the parameters is illustrated in Figs. 13 and 14. In particular, Figure 13 shows again the limited impact of an increase of $m$ on the Nadir in comparison to an increase of $d$. Similarly, increasing $r^{-1}$ seems to have a more significant impact than a reduction in $\tau$. Moreover, while both $d$ and $r^{-1}$ have a direct impact on the final steady-state (which constitutes a bound on the Nadir), neither $m$ nor $\tau$ affect it, as expected.

We further illustrate in Fig. 15 the effect of individually varying parameters $d$, $r$ and $m$ in the scenario of Fig. 10. The most significant impact both in time characteristics of $\overline{w}(t)$ (Nadir, steady state) and on synchronization cost is increasing the damping $d$. The droop coefficient $r^{-1}$ has some positive influence in the time response, less so the inertia $m$; neither of the two reduces the synchronization cost.

As to the role of inertia, it is illustrative to compare the bottom plot of Fig. 15 with Fig. 10. For the system frequency, we see that the increased inertia produces a very slight reduction in the deviation at the Nadir. The oscillatory components have been reduced in amplitude (the largest almost in half) but now they last longer in time; our $\mathcal{L}_{2}$ cost is largely insensitive. Intuitively speaking: a heavier system initially deviates less, but it takes more effort to subsequently control it to settle down.

VII-C Non-proportionality and Impact of Connectivity

So far we have simulated a network that corresponds to the real-world Icelandic grid, both electrically and in the machine inertias, but where other machine parameters were altered to satisfy Assumption 2. What happens if we simulate with the real, unmodified data? Figure 16 represents the response under the same scenario as in Figure 10, but using instead the true machine parameters. We observe the close similarity of both sets of plots, indicating that our modeling technique is capturing the essential dynamics. We note that there is a certain reduction in our synchronization cost; real data behave somewhat better than the proportional version.

Another issue to explore is the effect of network connectivity; in Section VI we showed that increased connectivity (reflected in larger eigenvalues for $L_{F}$) aligns the response around the system frequency, and reduces the impact of uncertainty due to non-proportionality. To investigate this we randomly add lines to the Icelandic power network with impedances uniformly distributed between its minimum and maximum values. Figure 18 shows some of the resulting network topologies; here $K$ denotes the number of lines added, $K=0$ in Figure 18(a) being the original topology.

Figure 17 shows the mean synchronization cost for proportionally correlated disturbances ($\Sigma_{u}=F^{2}$) for the true grid parameters as well as the proportional parameters computed using (48)-(49) for the swing dynamics with turbines. We can see how not only the synchronization error goes to zero as the connectivity increases (as expected) but also the relative error between the synchronization cost with proportional parameters and real parameters remains within 10%.

VIII Conclusions

In this paper we have analyzed power grid synchronization, with the explicit aim of reconciling system theoretic tools with power engineering insights and practice. With this motivation, we favored performance metrics based on the step response, and attempted to isolate two components: system wide frequency deviations characterized by time domain specifications, and transient oscillatory deviations measured by an $\mathcal{L}_{2}$ norm. Our analysis obtains such a decomposition for networks with heterogeneous machine ratings, provided a certain proportionality condition holds. Under this assumption, closed-form performance metrics are given in terms of the network, the ratings, and the parameters of a representative machine. We also applied robustness tools to analyze deviations from the proportionality assumption, and presented simulation evidence on a real set of data to validate the relevance of our analysis.

A first observation is that machine models matter: to replicate observed behaviors such as overshoots in mean frequency response, requires enhancing classical swing equation model with a some representation of turbine controls. A less obvious conclusion is the rather minor role of system inertia in these performance metrics: while inertia clearly influences the speed of initial ramping (RoCoF), its impact on the maximal deviation (Nadir) of the mean system frequency, and on the oscillation cost is very limited. Details depend on the machine model employed, but on the whole it appears that inertia is not the main parameter on which synchronization performance rests. This should assuage some concerns (e.g. [42]) regarding the decreased inertia of power grids with high penetration of renewable energy sources.555Of course, this analysis assumes some inertia is still present; in the extreme case of a grid where all generation is connected by power electronics, other issues like voltage dynamics are likely to dominate.

What is the key parameter, then? Both our theoretical results and our simulations point to the damping $d$ in the swing equation model as the main “knob” to control synchronization 666Consistently, the recent reference [41] that ostensibly focuses on low inertia ends up advocating for an ”augmented frequency damping” as a solution.. The damping $d$ is mainly due to the response of the local load to variations in frequency, which suggests that active load control is the best frequency regulation tool to mitigate step disturbances. Indeed, recent literature [28, 30, 37, 43, 44] has explored the potential of this method, also enabled by smarter grids. Our results lend support for this control strategy for improving synchronization in future power networks.

If more active control is enabled, one could choose to go beyond the emulation of physical parameters such as damping or inertia, and use a more general dynamic controller; indeed such approaches are also being considered [13, 17].

Appendix A Inner product computation

We show here how to evaluate the inner product in (17) using state-space methods. We start with a minimal state-space realization of the step response of the representative machine,

[TABLE]

If follows easily that each $\tilde{h}_{k}$ in (13) has realization

[TABLE]

Note that the state matrix $A_{k}$ is the only one that depends on the eigenvalue $\lambda_{k}$ under consideration. Under Assumption 1, $A_{k}$ is a Hurwitz matrix for any $\lambda_{k}>0$.

Writing $\tilde{h}_{k}(t)=Ce^{A_{k}t}B$ for the impulse response, we compute the inner product between two such functions:

[TABLE]

here T denotes matrix transpose. A standard calculation shows that $Q_{kl}$ satisfies the Sylvester equation

[TABLE]

Furthermore since the eigenvalues of $A_{k},A_{l}$ never add up to zero it follows (see [45]) that (51) has a unique solution $Q_{kl}$. Thus, the relevant inner product can be found by solving the above linear equation and substituting into (50).

Second order machine model

In the situation of Section IV, it is easily checked that

[TABLE]

In this case the solution to the Sylvester equation is

[TABLE]

Substitution into (50) for the given $C$ proves Proposition 3.

Third order machine model

Here the relevant matrices are

[TABLE]

The Sylvester equations for $Q_{kl}$ in this case (9 linear equations, 9 unknowns) give unwieldy expressions. We report first the case $k=l$, which remains tractable; here we have a Lyapunov equation with symmetric solution

[TABLE]

where $\Delta=2[m(r^{-1}+d)+\tau d(r^{-1}+d+\lambda_{k}\tau)].$ By looking at the $(1,1)$ element of this matrix we find the norm $\|\tilde{h}_{k}\|^{2}=\langle{\tilde{h}_{k}},{\tilde{h}_{k}}\rangle$, which coincides with the expression given in (42).

Going now to the general case $k\neq l$, we report here only the inner product obtained from the (1,1) element of $Q_{kl}$, itself found by solving the Sylvester equation using the Matlab symbolic toolbox:

[TABLE]

where

[TABLE]

The limiting cases $m\to 0$ and $m\to\infty$, presented in Section V were obtained from this general formula.

Appendix B Time-domain results with turbine control

Proof:

Since the second order system $g_{0}(s)$ is stable, the maximum of the impulse response is either at $t=0$ or the first time $\dot{g_{0}}(t)=0$. Since $\bar{g}_{0}(0)=0$, then it must be the latter.

Therefore we need to find the time instance $t_{\text{Nadir}}$ such that $|g_{0}(t_{\text{Nadir}})|=\underset{t\geq 0}{\sup}|g_{0}(t)|$.

The time derivative of (36) is given by

[TABLE]

where $\phi$ is defined as in the (40).

Setting now $\dot{g}_{0}(t)=0$ in (54) gives

[TABLE]

and since $\phi+\frac{\pi}{2}>0$, the first maximum is for $k=0$. Therefore

[TABLE]

which after substituting in (36) gives

[TABLE]

where the last step follows from (35), (37), (40) and

[TABLE]

as well as several trigonometric identities.

Therefore the Nadir of the system frequency is given by

[TABLE]

∎

Proof of Proposition 5

The proof is based on two previous lemmas.

Lemma 8.

Given a power system under Assumption 2 with $g_{i}(s)$ given by (4), the derivative of $\phi$ with respect to $m$ is given by

[TABLE]

Proof:

Notice that while it is not possible to derive a closed form condition for $\phi$ in terms of the system parameters without the aid of a trigonometric function, it is possible to achieve such expression for $\frac{\partial\phi}{\partial m}$. This is achieved by first computing

[TABLE]

Now using the fact that $\phi\in(-\frac{\pi}{2},\frac{\pi}{2})$ we get

[TABLE]

∎

Lemma 9.

Given a power system under Assumption 2 with $g_{i}(s)$ given by (4), the derivative of $\frac{\eta}{\omega_{d}}$ with respect to $m$ is given by

[TABLE]

Proof:

The proof is just by direct computation

[TABLE]

∎

Proof:

Using lemmas 8 and 9, and $\alpha=\frac{\pi}{2}+\phi$ we can compute

[TABLE]

It is easy to see that whenever $m\geq d\tau$, then $\frac{\partial}{\partial m}\bar{w}<0$. However it is possible to show that even when $m<d\tau$, the remaining negative terms are still dominant making $\frac{\partial}{\partial m}\bar{w}<0$ always true. ∎

Proof of Proposition 6

Again we first state a lemma with some of the calculations.

Lemma 10.

Given $g_{0}(s)$ as in (4)

[TABLE]

where the angle $\beta\in(-\frac{\pi}{2},\frac{\pi}{2})$, $\beta>\phi$ for $\phi\in(-\frac{\pi}{2},\frac{\pi}{2})$ defined by (55) and (40), and $\beta$ is uniquely defined by

[TABLE]

and

[TABLE]

Proof:

We first compute

[TABLE]

where

[TABLE]

Thus we can compute

[TABLE]

Therefore, it is enough to compute

[TABLE]

∎

Proof:

When $t=0^{+}$, $\ddot{g}_{0}(0)<0$. Therefore, the initial trend the $\dot{g_{0}}(t)$ is decreasing and therefore

[TABLE]

is a local maximum.

We will now show that this is indeed a global maximum. The function $\ddot{g}_{0}(t)$ crosses zero for the first time when $\omega_{d}t-\beta=\frac{\pi}{2}$ or equivalently

[TABLE]

By substituting (60) into (54) and defining

[TABLE]

we get

[TABLE]

Finally by simplifying

[TABLE]

we obtain

[TABLE]

Therefore, since for any additional instant of time that $\dot{g}_{0}(t)$ achieves a local extremum $t_{k}^{*}=\frac{\beta+\frac{\pi}{2}+k\pi}{\omega_{d}}$, the factor $|\cos(\omega_{d}t_{k}^{*}-\phi)|$ does not change, then the maximum is achieved at $t=0$.

∎