Global performance metrics for synchronization of heterogeneously
rated power systems: The role of machine models and inertia

arXiv:1710.07195v2 [math.OC] 24 Oct 2017

Fernando Paganini, Fellow, IEEE, and Enrique Mallada, Member, IEEE

Abstract— A recent trend in control of power systems has
sought to quantify the synchronization dynamics in terms of
a global performance metric, compute it under very simplified
assumptions, and use it to gain insight on the role of system
parameters, in particular, inertia. In this paper, we wish to extend this approach to more realistic scenarios, by incorporating
the heterogeneity of machine ratings, more complete machine
models, and also to more closely map it to classical power
engineering notions such as Nadir, Rate of Change of Frequency
(RoCoF), and inter-area oscillations.
We consider the system response to a step change in power
excitation, and define the system frequency as a weighted
average of generator frequencies (with weights proportional
to each machine’s rating); we characterize Nadir and RoCoF
by the L∞ norm of the system frequency and its derivative,
respectively, and inter-areas oscillations by the L2 norm of the
error of the vector of bus frequencies w.r.t. the system frequency.
For machine models where the dynamic parameters (inertia,
damping, etc.) are proportional to rating, we analytically
compute these norms and use them to show that the role of
inertia is more nuanced than in the conventional wisdom. With
the classical swing dynamics, inertia constant plays a secondary
role in performance. It is only when the turbine dynamics are
introduced that the benefits of inertia become more prominent.

I. I NTRODUCTION
The synchronization performance of the power grid has
been a major concern of system operators since the early
days [1], [2]. Most generators and loads are designed on the
assumption that the grid frequency is tightly regulated around
a nominal value (e.g. 60Hz in the U.S., 50Hz in Europe).
When the frequency deviates significantly (a few hundred
mHz) due to some network fault, several mechanisms, such
as machine protections or under frequency load shedding
(ULFS) [3], automatically disconnect critical network elements, potentially causing cascading failures and ultimately
blackouts [4].
The gradual substitution of conventional electromechanical with renewable generation has raised new concerns about
synchronization performance. The former provide a natural
response to power imbalances, which is not present in the
inverter-based interfaces of renewable sources; in particular
the lack of inertia in the latter is often seen as a threat
to frequency regulation. An objective analysis of this issue
requires identifying appropriate performance metrics.
Two separate mechanisms are at play in frequency fluctuations. On one hand, abrupt changes on the global
Fernando Paganini is with Universidad ORT Uruguay, Montevideo. email: paganini@ort.edu.uy. Enrique Mallada is with Johns Hopkins University, Baltimore. e-mail: mallada@jhu.edu. This work was supported in part
by IDB/MIEM-Uruguay, Project ATN/KF 13883 UR, ANII-Uruguay, grant
FSE 1 2016 1 131605, and NSF, grants CNS 1544771, EPCN 1711188, and
AMPS 1736448.

supply-demand balance induce system-wide frequency
changes, which may persist in steady state. On the
other hand, geographically-distributed frequency oscillations
(a.k.a. inter-area oscillations) [5], [6], [7] may be observed
due to weak global coupling. Good performance metrics
should be able to identify and discriminate between these
two phenomena.
To address these problems, power engineers have traditionally relied on classical control metrics. For example,
they use the system response to a step input to measure the
maximum frequency deviation to imbalances (Nadir) as well
as the maximum rate of change of frequency (RoCoF)[8];
their main limitation is that these quantities are node dependent. To evaluate inter-area oscillations, eigenvalue methods
(slow coherency [9], participation factors [10]) have been
employed.
More recently, a trend in control of power systems has
aimed to quantify synchronization performance in terms of
global system metrics such as H2 or H∞ [11], [12], [13],
[14], [15], [16], [17]. These norms capture the effect of system parameters, such as inertia, damping, and the eigenvalues
of the network (Laplacian) matrix, on system performance.
However, closed-form analytic results depend on oversimplified assumptions –homogeneous machines modeled by
swing equations– and do not directly represent step-response
information which is most important for network operators
to analyze disturbances.
In this paper, we wish to bridge the gap between the two
approaches. Firstly, to extend the latter approach to cover
more realistic scenarios by incorporating heterogeneity of
machine ratings, and adding turbine dynamics. Secondly,
to incorporate step-response metrics (Nadir, RoCoF) for
an appropriately defined global system frequency, and to
separately characterize inter-area oscillations.
To obtain analytically tractable results, we focus on a
specifically family of heterogeneous machines, in which
key dynamic parameters are proportional to rating. This
restriction is mild in comparison with homogeneity, and
enables a diagonalization procedure, generalizing traditional
eigen-analysis. From it, a system frequency suitable for step
response analysis appears naturally, and turns out to be
P
mi wi (t)
,
(1)
w̄(t) := iP
i mi

the weighted average of node frequencies in proportion to
their inertia. This is identical to the frequency of the center
of inertia (COI) a classical notion [18], [19]. Nadir and
RoCoF are defined as the L∞ norms of, respectively, w̄(t)

˙
and w̄(t).
A synchronization cost measuring transient interarea oscillations is defined as the L2 norm of the vector of
deviations w̃i (t) = wi (t) − w̄(t).
The rest of the paper is organized as follows. In Section
II we formulate the model and carry out the diagonalization
of the dynamics, making the above decomposition precise.
The system frequency step response, and a closed form
expression for the synchronization cost, are both expressed in
terms of a representative machine and the network structure.
In Sections III and IV 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. We find that some important aspects, in
particular the importance of inertia, depend crucially on the
chosen model. The paper concludes with a discussion in
Section V, and some derivations are covered in the Appendix.
II. DYNAMIC

MODEL WITH MACHINE HETEROGENEITY

This corresponds to the transfer function
gi (s) =

mi

1
.
+ di s

s2

(2)

Example 2: The swing equation with a first-order model
of the turbine control:
θ̇i = wi ,
mi ẇi = −di wi + qi + ui − pei ,
τi q̇i = −ri−1 ω − qi .

Here qi is the (variation of) turbine power, τi the turbine time
constant and ri the droop coefficient; the governor dynamics
are considered to be faster and neglected. The corresponding
transfer function is
τi s + 1
 . (3)
gi (s) =
s mi τi s2 + (mi + di τi )s + di + ri−1

Of course, other models are possible within this framework
(e.g. a 4th order system including a state for the governors).

G = diag(gi )

A. A family of heterogeneous machines

θ

u
+

-

pe
Fig. 1.

L

Block Diagram of Linearized Power Network

We consider a set of n generator buses, indexed by i ∈
{1, . . . , n}, dynamically coupled through an AC network.
Assuming operation around an equilibrium, the linearized
dynamics are represented by the block diagram in Fig. 1,
where:
• G(s) = diag(gi (s)) is the diagonal transfer function
of generators at each bus. Each gi (s) has as output
the phase angle θi , and as input the net power at its
generator axis, relative to its equilibrium value. This
includes an outside disturbance ui , reflecting variations
in mechanical power or local load, minus the variation
pei in electrical power drawn from the network.
• Using a linear DC model for the network, the vector
of drawn power is written as pe = Lθ, where L is the
weighted Laplacian of the graph defined by the line
susceptances. Thus all the coupling between the bus
subsystems is through this feedback term. L is a rank
n − 1 matrix with kernel spanned by 1, the vector of
all ones.
Two examples of generator dynamics, to be considered
explicitly in this paper, are:
Example 1: The swing equation dynamics
θ˙i = wi ,
mi ẇi = −di wi + ui − pei .

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 mi , di , etc.), a very
restrictive scenario.1
In a real network, where generators have different power
ratings, it is natural for parameters to scale accordingly: for
instance, the inertia mi 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 show here that a compact analysis can be given for
the case where parameters satisfy a certain proportionality.
We formalize this by introducing a rating parameter 0 <
fi ≤ 1, defined in relation to the largest machine which has
fi = 1, and imposing the following:
Assumption 1: There exists a fixed transfer function g0 (s),
termed the representative machine, such that
gi (s) =

1
g0 (s)
fi

for each i, where fi > 0 is the rating parameter of bus i.
To interpret this, consider first the swing dynamics of
Example 1. Then the assumption is satisfied provided that
inertia and damping are both proportional2 to fi , i.e. mi =
fi m, di = fi d, where m, d are those of the largest
(representative) machine,
g0 (s) =

1
.
+ ds

ms2

1 Some bounds on heterogenous systems are given in [11], [14]. Numerical
studies with heterogeneity are given in [12].
2 This kind of proportionality is termed “uniform damping” in [19].

Going to the case of Example 2 with the turbine dynamics,
we find that the above assumption is satisfied provided that
mi = fi m, di = fi d, ri−1 = fi r−1 , τi = τ ; here the inverse
droop coefficient is assumed proportional to rating, but the
turbine time-constant is taken to be homogeneous.
The corresponding representative machine is
τs + 1
.
g0 (s) =
s (mτ s2 + (m + dτ )s + d + r−1 )
Regarding the practical relevance of our simplifying assumption: empirical values reported in [20] 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.
B. Diagonalization
We will now exploit the above assumption to transform
the dynamics of Fig. 1 in a manner that allows for a suitable
decoupling in the analysis. In what follows, F = diag(fi )
denotes the diagonal matrix of rating parameters. Writing
1

1

G(s) = diag(gi (s)) = F − 2 [g0 (s)I]F − 2 ,
we transform the feedback loop into the equivalent form of
Fig. 2.
u0

1

F -2

1

g0(s)I

+

-

w

F -2

u0

1

F -2

VT

+

-

V

g0(s)I

VT

F
Fig. 2.

L

F

We introduce a notation for the scaled Laplacian matrix3
LF := F

LF

− 12

,

(4)

which is positive semidefinite and of rank n − 1. Applying
the spectral theorem we diagonalize it as
LF = V ΛV T ,

(5)

where Λ = diag(λk ), 0 = λ0 < λ1 ≤ · · · ≤ λn−1 , and V
is unitary. Distinguishing the eigenvector v0 that corresponds
to the zero eigenvalue, we write V = [v0 V⊥ ] , where
V⊥ ∈ Rn×(n−1) , V⊥T V⊥ = In−1 , V⊥T v0 = 0.
In fact v0 can be made explicit by recalling that ker(L) =
1
span{1}, so ker(LF ) = span{F 2 1}, from where
 X − 21
1
v0 = αF F 2 1,
.
(6)
fi
with αF :=
i

Substitution of (5) into Fig. 2 and some block manipulations leads to the equivalent representation of Fig. 3.
Noting finally that the V block commutes with g0 (s)I
and thus cancels out with V T , the internal loop is now fully
diagonalized, yielding the closed-loop transfer function
3 This

w

Fig. 3. Equivalent block diagram for heterogeneously rated machines with
diagonalized closed loop

H(s) = diag(hk (s)), with
g0 (s)
hk (s) =
, k = 0, 1, . . . , n − 1.
1 + λk g0 (s)

(7)

Assumption 2: The proportional feedback λk > 0, k =
1, . . . , n − 1 is stabilizing for g0 (s), i.e. hk (s) has all its
poles in Re[s] < 0 for k = 1, . . . , n − 1.
Later on we will verify that the assumption always holds
for the two examples considered.
As for h0 (s), we see that here there is no feedback:
h0 (s) = g0 (s). In the examples above, g0 (s) has a pole at
s = 0, which corresponds to the integration from frequency
to angle. This is the typical situation in which the local
control of the machine has no angle feedback, the latter only
appears when considering coupling through the network.
Returning to Fig. 3 we arrive at the transfer function
1

Tθu (s) = F − 2 V H(s)V T F − 2

-12

Equivalent block diagram for heterogeneously rated machines

− 12

1

F -2

Λ

1

-12

V

scaling already appears in the classical paper [9].

between the vector of external power disturbances and the
machine angle outputs. If the output of interest is chosen to
be the vector of frequencies, the relevant transfer function is
Twu (s) = sTθu (s).
C. Step response characterization
Global metrics for synchronization performance in e.g.
[11], [14], [15], are system norms (H2 , H∞ ) applied to Twu
(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 L2 signal).
In this paper, we wish to bring metrics closer to industry
practice, by considering a step input disturbance (e.g. due
a fault), and analyzing the response of the vector of bus
frequencies wi (t). Practitioners are interested in some timedomain performance metrics (Nadir, RoCoF, see below) for
the response, but transient frequency is bus-dependent.
A candidate global notion of system frequency is the
weighted average w̄(t) in (1). We now show that for our
family of heterogenous systems, the behavior of w̄(t) decouples nicely from the individual bus deviations wi (t) − w̄(t),
opening the door for a separate analysis of both aspects.
Our input will be a step function u(t) = u0 1t≥0 ; here u0
is a given vector direction. In Laplace transforms:
1
1
1
w(s) = Twu (s) u0 = F − 2 V H(s)V T F − 2 u0 .
s

(8)

We now isolate the dynamics corresponding to the eigenvalue
λ0 = 0 and its eigenvector v0 from the rest:
w(s) = F

− 12

v0 h0 (s)v0T F

− 21

u0 + F

− 12

V⊥ H̃(s)V⊥T F

− 12

u0 ;

here H̃(s) = diagk=1,..,n−1 (hk (s))).
1
Noting that F − 2 v0 = αF 1 from (6), and h0 (s) = g0 (s),
the first term above is of the form w̄(s)1, where
P
u0i
2 T
g0 (s).
(9)
w̄(s) := αF 1 u0 g0 (s) = Pi
i fi

The second term will be denoted by w̃(s); by Assumption
2 it has left half-plane poles. So we have obtained the
decomposition
w(t) = w̄(t)1 + w̃(t),

(10)

•

We now turn our attention to the term 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 L2 norm
Z ∞
kw̃k22 =
|w̃(t)|2 dt.
0

We now show how this norm can be computed in terms
of the parameters of the scaled network Laplacian, and
the impulse response matrix H̃(t) = diagk=1,...,n−1 (hk (t)),
Laplace inverse of H̃(s), which encapsulates all information
on the machine model.
Proposition 1: kw̃k22 = z0T Y z0 , where:
(n−1)×(n−1)
• Y ∈R
is the matrix with elements
Z ∞
ykl = γkl hhk , hl i = γkl
hk (t)hl (t)dt,
(12)
0

intepreted as follows:
•

E. Quantifying the deviation from synchrony

w̄(t) is a system frequency term, applied to all nodes;
the transient term w̃(t) represents the individual node
deviations from the synchronous response.

where Γ = (γkl ) := V⊥T F −1 V⊥ ;

(13)

1

•

z0 := V⊥T F − 2 u0 ∈ Rn−1 .
Proof: With the introduced notation we have
1

w̃(t) = F − 2 V⊥ H̃(t)z0 ,

D. System frequency
We can obtain more information on the system frequency
1
T
by observing that since 1T F 2 V⊥ = α−1
F v0 V⊥ = 0, we have
1

1

(1T F ) F − 2 V⊥ H̃(s)V⊥T F − 2 u0 ≡ 0.
{z
}
|

therefore w̃(t)T w̃(t) = z0T H̃(t)ΓH̃(t)z0 . The matrix in the
above quadratic form has elements hk (t)γkl hl (t), therefore
integration in time yields the result.
Remark 2: The metric kw̃k22 does depend on the electrical
network, through the eigenvalues and eigenvectors of LF .

w̃(s)

Therefore 1T F w(t) = w̄(t)1T F 1 by (10), which gives
P
fi wi (t)
;
w̄(t) = iP
i fi

the system frequency is a weighted mean of bus frequencies,
in proportion to their rating. Noting that mi = mfi , it follows
that w̄(t) is exactly the COI frequency from (1).
Also, returning to (9) we have
P
u0i
g0 (t).
(11)
w̄(t) = Pi
i fi

Recall that g0 (t) is the angle impulse response of the representative machine, or equivalently its angular frequency step
response. Thus w̄(t) corresponds to the frequency observed
when exciting the representative machine (in open loop) with
the total system disturbance normalized by the total scale.
Remark 1: 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 (w̃) converges fast
to zero, (11) is a natural candidate for a reduced order model
similar to the ones recently considered in [21], [22].
In the following sections we analyze its behavior for the
previously discussed examples.

F. Mean synchronization cost for random disturbance step
Since the cost discussed above is a function of the
disturbance step u0 , it may be useful to find its average
over a random choice of this excitation. Recalling that the
components u0i correspond to different buses, it is natural
to assume them to be independent, and thus E[u0 uT0 ] = Σu ,
a diagonal matrix. Therefore
E[z0 z0T ] = V⊥T F −1 Σu V⊥ =: Σz ,
and the expectation for the cost in Proposition 1 is


E kw̃k22 = E[z0T Y z0 ] = E[Tr(Y z0 z0T )] = Tr(Y Σz ).

We look at some special cases:
z
• Σu = I (uniform disturbances). Then Σ = Γ, and
X


2
hhk , hl i.
γkl
E kw̃k22 = Tr(Y Γ) =
k,l

•

Σu = F . This means disturbance size follows the square
root of the bus rating. Here Σz = I, and
X


E kw̃k22 = Tr(Y ) =
γkk khk k22 .
k

•

Σu = F 2 . This is probably most natural, with disz
turbances proportional to bus rating.
Here Σ =

V⊥T F V⊥ = Γ† (pseudoinverse); E kw̃k22 = Tr(Y Γ† ).

G. The homogeneous case

•

If all machines have the same response g0 (s), setting F =
I we can obtain some simplifications:
• {λk } are the eigenvalues of the original Laplacian L.
• The system frequency is the average w̄(t)
=
1 P
w
(t),
and
satisfies
i i
n
!
1 X
u0i g0 (t).
w̄(t) =
n
i
•

z0 = V⊥T u0 , and Γ = V⊥T V⊥ = I. Therefore the matrix
Y in Proposition 1 is diagonal, Y = diag(khk k2 ), and
kw̃k22 =

•

n−1
X
k=1

(z0k )2 khk k2 ,

(14)

where z0k = vkT u0 is the projection of the excitation
vector u0 in the direction of the k-th eigenvector of the
Laplacian L.
The mean synchronization cost for E[u0 uT0 ] = I (here
all the preceding cases coincide) is


E kw̃k22 =

n−1
X
k=1

khk k22 = kH̃k2H2 ,

(15)

the H2 norm of the transfer function H̃(s). We recall
that this was obtained by isolating the portion g0 (s) corresponding to the synchronized response (in this case,
projecting onto 1⊥ ). 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.
III. A PPLICATION

TO THE SWING DYNAMICS

In this section we assume we are in the situation of
Example 1, i.e., the representative machine is
g0 (s) =

1
.
s(ms + d)

and mi = fi m, di = fi d are the individual bus parameters.
The corresponding closed loop transfer functions in (7) are
1
,
hk (s) =
ms2 + ds + λk

k = 0, 1, . . . , n − 1.

A. System frequency
Inverting the transform and invoking (11) we find that
g0 (t)

z
}|
{

d
1
−m
t
1−e
d


d
1 − e− m t , t > 0.

•

here the total inertia appears, which is natural in a
second-order 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∞ , independent
of m.
B. Synchronization cost
We now turn to the synchronization cost kw̃k2 , which can
be computed by particularizing the result in Proposition 1.
The following result is proved in Appendix I.
Proposition 2: Let hk (s) be given in (16), and hk (t) its
inverse transform, for k = 1, . . . , n − 1. Then:
hhk , hl i =

(17)
(18)

Some comments are in order:
• Again, we recall w̄(t) does not depend on the electrical
network.

2d
.
m(λk − λl )2 + 2(λk + λl )d2

(20)

It follows that the matrix Y in (12) 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 → 0 or m → ∞; we thus argue that inertia
has limited impact. We look at this issue in further detail.
1) Homogeneous case: In the case of homogeneous machines, we saw above that Γ = I and Y is diagonal, so
inertia disappears completely: indeed using (14) we have
kw̃k22 =

(16)

Note they are stable when λk > 0.

P
u0i
w̄(t) = Pi
f
Pi i
u
0i
= Pi
i di

•

The first-order evolution of w̄(t) implies there is no
overshoot; system frequency never deviates to a “Nadir”
further from equilibrium than its steady-state
value.
P
u
The asymptotic frequency w∞ = Pi d0ii is the ratio
i
of total disturbance to total damping, but does not
depend on the inertia m. The latter only affects the time
constant in which this asymptote is achieved.
The maximum RoCoF (rate-of-change-of frequency)
occurs at t → 0+, and is given by
P
P
u0i
d i u0i
P
= Pi
;
(19)
m
i di
i mi

n−1
X
k=1

(vkT u0 )2
.
2dλk

(21)

The cost is inversely proportional to damping, and the
direction of the disturbance u0 also matters. Recalling that
vk is the k-th Laplacian eigenvector, the worst-case for a
given magnitude |u0 | is when it is aligned to v1 , the Fiedler
eigenvector.
If the disturbance direction is chosen randomly as in
Section II-G, then (15) gives

 X
1
1 X 1
=
Tr(L† );
E kw̃k22 =
khk k22 =
2d
λk
2d
k

(22)

k

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

2) Heterogeneous, high inertia case: Assume for this
discussion that all the λk are distinct; then as m → ∞ we
have ykl → 0 for k 6= l, so Y again becomes diagonal, and
the cost has the limiting expression
m→∞
kw̃k22 −→

n−1
X
k=1

2
γkk z0k
.
2dλk

(23)

So the high inertia behavior is of a similar structure to the
homogeneous case in (21). Comparisons are not straightforward, though, since the scaling factor F affects z0k , γkk and
λk in each of the above terms.
3) Heterogeneous, low inertia case: If m → 0, then the
limiting Y matrix is not diagonal. The corresponding limiting
cost is
m→0+

kw̃k22 −→

n−1
X

k,l=1

γkl z0k z0l
.
d(λk + λl )

TURBINE DYNAMICS

We now turn to the model of Example 2, where the
representative machine is
g0 (s) =

τs + 1
.
s (mτ s2 + (m + dτ )s + d + r−1 )

(25)

The corresponding closed loop transfer functions in (7) for
k = 0, 1, . . . , n − 1 are:
hk (s) =

the system is under-damped with poles η ± jωd , and



1
1
s+γ
g0 (t) = L−1
−
d + r−1 s (s + η)2 + ωd2



(γ −η)
1
−ηt
cos(ωd t)−
1−e
sin(ωd t)
=
d + r−1
ωd
(29)
where

(24)

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 L2 price of synchrony, even in the heterogeneous
machine case.
IV. M ODEL WITH

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

2
1 1
d
d + r−1
=: ωd2 > 0
(28)
−
+
mτ
4 τ
m

τs + 1
.
mτ s3 + (m + dτ )s2 + (d + r−1 + λk τ )s + λk
(26)

It can be checked (e.g. by applying the Routh-Hurwitz
criterion) that hk (s) is stable whenever λk > 0.
A. System frequency
We can again use (11) and (25) to compute the system
frequency
P
u0i −1
L {g0 (s)} ,
(27)
w̄(t) = Pi
i fi
but now the inverse transform of g0 (s) is more involved.
Using partial fractions we first express




r −1
1
−
s
+
τ
m
1

1 −
g0 (s) =

−1
d
d + r−1 s s2 + τ1 + m
s + d+r
mτ

The first term provides the steady-state response, which is
P
P
1
i u0i
i u0i
P
;
w∞ =
=
−1
f
d
+
r
d
+
ri−1
i
i i

this is analogous to the swing equation case, except than the
droop control has been added to the damping. Again, we
observe that inertia plays no role at all in this steady-state
deviation.

η :=

1
2



1
d
+
τ
m



and γ :=



1 r−1
−
τ
m



.

(30)

The system frequency time evolution is given by
P
u0i
w̄(t) = P i −1 g0 (t),
i di +ri

(31)

with g0 (t) from (29).
A few observations are in order:
• Including the turbine model has a nontrivial effect on
the system frequency w̄(t). It is the presence of the
turbine dynamics that provides the characteristic underdamped behavior that produces a Nadir.
• We have only provided here the solution of w̄(t) for the
(practically more relevant) under-damped case.
• Interestingly, (28) 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.
1) Nadir:
In order to compute the Nadir we will use
P
u0i
||g0 ||∞ .
||w̄||∞ = Pi
i fi

Thus, one can compute the Nadir by finding the maximum
excursion of g0 (t). The following proposition summarizes
the overall calculation which can be find in Appendix II.
Proposition 3 (Nadir): Given a power system under Assumption 1 with generators containing first order turbine
dynamics (gi (s) given by (3)). Then under the under-damped
condition (28), the Nadir is given by
!
r
P
| i u0i | 1
τ r−1 −ωη (φ+ π2 )
,
1+
e d
||w̄||∞ = P
−1
m
i fi d+r
(32)

where the phase φ ∈ (− π2 , π2 ) is uniquely determined by

1
−η
m − dτ
= √
.
(33)
sin(φ) = q τ

2
1
2 mτ r−1
ωd2 + τ − η

The dependence of (32) on m is not straightforward, as
φ, η, and ωd depend on it. The next proposition shows
that the dependence is as expected by conventional power
engineering wisdom.
Proposition 4: Given a power system under Assumption 1
with generators containing first order turbine dynamics (gi (s)
given by (3)). Then under the under-damped condition (28),
the maximum frequency deviation ||w̄||∞ is a decreasing
∂
||w̄||∞ < 0.
function of m, i.e., ∂m
The proof can be found in Appendix III.
2) RoCoF:
A similar procedure as the one used to study the Nadir of the
system frequency can be used to investigate the properties
of the maximum rate of change of frequency (RoCoF).
Proposition 5 (RoCoF): Given a power system under Assumption 1 with generators containing first order turbine
dynamics (gi (s) given by (3)). Then under the under-damped
condition (28), the RoCoF is given by
P
u | 1
|
˙ ∞ = Pi 0,i
.
(34)
||w̄||
i fi m
The proof of Proposition 5 can also be found in Appendix
IV; the main difference with Proposition 3 is that, while not
trivial to establish, here the maximum is always achieved at
t = 0+, exactly as in the second order case of (19).
˙ ∞ on m is now easily addressed
The dependence of ||w̄||
and again as expected: RoCoF decreases with m.
B. Synchronization cost
The synchronization cost kw̃k2 can once again be computed through Proposition 1, which requires finding the inner
products hhk , hl i, in this case for the functions in (26).
Since the corresponding expression is in general rather
unwieldy (see Appendix I), we will present some simpler
cases, beginning with k = l; the norm, found in Appendix I
is:
m + τ (λk τ + d)
. (35)
khk k22 =
−1
2λk [m(r + d) + τ d(r−1 + λk τ + d)]
1) Homogeneous case: The above expression suffices to
analyze the case of homogeneous machines, where Γ = I
and Y is diagonal. We have from (14) that
kw̃k22 =

n−1
X
k=1

(vkT u0 )2 khk k22 ;

from (35) we see that, in contrast to the second order machine
model, the inertia m does affect the synchronization cost. A
closer look at khk k2 as a (linear fractional) function of m
shows that it is decreasing in m ∈ (0, ∞), going from
m→0+

khk k22 −→

λk τ + d
1
·
,
2λk d r−1 + λk τ + d

to
m→∞

khk k22 −→ =

d
1
·
.
2λk d r−1 + d

So higher inertia is beneficial in this case. Recalling that
the corresponding cost for the swing dynamics is 2λ1k d , we
see that this cost has been reduced. In the high inertia case,
the main change is the increased damping through the droop
coefficient r−1 .
2) Heterogeneous, high inertia case: As mentioned, the
formula for hhk , hl i for k 6= l is quite formidable, but we
can give its approximation in the limit of large m:
hhk , hl i

m→∞

∼

2(d + r−1 )
,
m(λk − λl )2

k 6= l.

This assumes λk 6= λl . So if the eigenvalues of the scaled
Laplacian LF are distinct, we see that again the matrix Y
becomes diagonal as m → ∞. The limiting cost is
m→∞

kw̃k22 −→

n−1
X
k=1

2
z0k
γkk
d
·
.
2λk d r−1 + d

This expression amounts to reducing to the cost (23) for the
second order dynamics, by the fraction r−1d+d . So the role of
the turbine in a high inertia system is again mainly a change
in the droop coefficient.
3) Heterogeneous, low inertia case: In the low inertia
m→0
limit, we find that hhk , hl i −→ N
D , where
N = 2d(d + r−1 ) + τ (2d + r−1 )(λk + λl ) + 2λk λl τ 2 ,

D = 2d(d + r−1 )2 (λk + λl ) + dτ (2d + r−1 )(λk + λl )2
+ 2dτ λk λl [2r−1 + τ (λk + λl )].
So the limiting matrix Y is not diagonal, as in the second
order case; an expression analogous to (24) can be written.
Comparisons between the two are not straightforward here,
and must be pursued by numerical experimentation.
V. C ONCLUSIONS
We have studied system-theoretic measures of synchronization performance, with the aim of covering more realistic
scenarios than the recent literature, and also closing the gap
with power engineering practice. In particular, for a family
of heterogeneous machine systems, we have focused on the
step response of the bus frequency vector, decomposed as
a system-wide weighted mean and the vector of relative
differences to it. The key mathematical assumption is a
relative proportionality in machine parameters, which may
be incorporated together with the network model to perform
an adequate diagonalization.
With this transformation, the natural system frequency
(motion of the COI) becomes independent of the network,
and its characteristics can be studied through a representative
machine. The energy of the synchronization error around this
mean depends on both network and machine models, but we
have an expression that encapsulates the latter in terms of
a matrix of inner products. From this general result special
cases can be studied.
A key question of interest to practitioners is the role of
inertia, in particular whether low inertia can compromise
performance. Our analysis shows that if a second-order,
swing equation model is used for each machine, the impact
of inertia is small. The global system frequency exhibits no

overshoot, inertia affecting mainly its time constant; while
inertia appears in the energy of oscillations, its impact is
not significant. The story changes if a dynamic model of
turbine control is adopted, where the turbine time-constant
is in the order of magnitude of the swing dynamics. In that
case, inertia does play a positive role, providing resilience of
the peak system frequency deviation and reducing the norm
of relative oscillations.
Future work will involve carrying out these calculations
for actual networks with real parameters.
A PPENDIX I
I NNER PRODUCT COMPUTATION
We show here how to evaluate the inner product in
(12) using state-space methods. We start with a state-space
realization of the representative machine,


A B
.
g0 (s) =
C 0
If follows easily that each hk in (7) has realization


Ak B
, where Ak = A − λk BC.
hk (s) =
C 0
Note that the state matrix Ak is the only one that depends
on the eigenvalue λk under consideration. By Assumption 2,
Ak is a Hurwitz matrix for any λk > 0.
Writing hk (t) = CeAk t B for the impulse response, we
compute the inner product between two such functions:
Z ∞
hk (t)hl (t)T dt
hhk , hl i =
0Z

∞
T
(36)
=C
eAk t BB T eAl t dt C T ;
0
{z
}
|
Qkl

here T denotes matrix transpose. A standard calculation
shows that Qkl satisfies the Sylvester equation
Ak Qkl +

Qkl ATl

T

+ BB = 0.

(37)

Furthermore since the eigenvalues of Ak , Al never add up to
zero it follows (see [23]) that (37) has a unique solution Qkl .
Thus, the relevant inner product can be found by solving the
above linear equation and substituting into (36).
Second order machine model
In the situation of Section III, it is easily checked that






0
1
0
B=
C= 1 0 .
Ak =
1 ;
d
−m
− λmk − m

In this case the solution to the Sylvester equation is


λk −λl
2d
1
2d
Qkl =
−λk
λk +λl ;
m(λk − λl )2 + 2(λk + λl )d2 λl 2d
2m
Substitution into (36) for the given C proves Proposition 2.

Third order machine model
Here the relevant matrices are




0
1
0
0


d
1 
1
; B = − m
; C= 1 0 0 .
Ak = − λmk − m
m
−1
0
0 − r τ − τ1

The Sylvester equations for Qkl 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


m+τ (λk +d)
0
−τ r−1
λ
k
1 

m+τ (r −1 +λk τ +d)
Qkk =

0
−r−1  ,
m
∆
−τ r−1
r−1
r−2
(38)
where ∆ = 2[m(r−1 +d)+τ d(r−1 +d+λk τ )]. By looking at
the (1, 1) element of this matrix we find the norm khk k2 =
hhk , hk i, which coincides with the expression given in (35).
Going now to the general case k 6= l, we report here
only the inner product obtained from the (1,1) element of
Qkl , itself found by solving the Sylvester equation using the
Matlab symbolic toolbox:

hhk , hl i =

N
,
D

where
N =2(2dm2 + 2m2 r−1 + 2d3 τ 2 + 4d2 mτ + 2d2 λk τ 3
+ 2d2 λl τ 3 + 2d2 r−1 τ 2 + 4dmr−1 τ + 2dλk λl τ 4
+ 2dλk mτ 2 + 2dλl mτ 2 + dλk r−1 τ 3 + dλl r−1 τ 3
+ λk mr−1 τ 2 + λl mr−1 τ 2 ),
D =4d4 λk τ 2 + 4d4 λl τ 2 + 4d3 λ2k τ 3 + 8d3 λk λl τ 3
+ 8d3 λk mτ + 8d3 λk r−1 τ 2 + 4d3 λ2l τ 3 + 8d3 λl mτ
+ 8d3 λl r−1 τ 2 + 4d2 λ2k λl τ 4 + 6d2 λ2k mτ 2
+ 2d2 λ2k r−1 τ 3 + 4d2 λk λ2l τ 4 + 4d2 λk λl mτ 2
+ 12d2 λk λl r−1 τ 3 + 4d2 λk m2 + 16d2 λk mr−1 τ
+ 4d2 λk r−2 τ 2 + 6d2 λ2l mτ 2 + 2d2 λ2l r−1 τ 3
+ 4d2 λl m2 + 16d2 λl mr−1 τ + 4d2 λl r−2 τ 2 + 2dλ3k mτ 3
− 2dλ2k λl mτ 3 + 4dλ2k m2 τ − 2dλk λ2l mτ 3 − 8dλk λl m2 τ
+ 16dλk λl mr−1 τ 2 + 8dλk m2 r−1 + 8dλk mr−2 τ

+ 2dλ3l mτ 3 + 4dλ2l m2 τ + 8dλl m2 r−1 + 8dλl mr−2 τ
+ 2λ3k λl mτ 4 + 2λ3k m2 τ 2 + λ3k mr−1 τ 3 − 4λ2k λ2l mτ 4

− 2λ2k λl m2 τ 2 − λ2k λl mr−1 τ 3 + 2λ2k m3 − 2λ2k m2 r−1 τ
+ 2λk λ3l mτ 4 − 2λk λ2l m2 τ 2 − λk λ2l mr−1 τ 3

− 4λk λl m3 + 4λk λl m2 r−1 τ + 4λk m2 r−2 + 2λ3l m2 τ 2
+ λ3l mr−1 τ 3 + 2λ2l m3 − 2λ2l m2 r−1 τ + 4λl m2 r−2 .

The limiting cases m → 0 and m → ∞, presented in Section
IV were obtained from this general formula.

A PPENDIX II
P ROOF OF P ROPOSITION 3

achieve such expression for
computing

Proof:
Since the second order system g0 (s) is stable, the maximum of the impulse response is either at t = 0 or the first
time g˙0 (t) = 0. Since ḡ0 (0) = 0, then it must be the latter.
Therefore we need to find the time instance tNadir such
that |g0 (tNadir )| = sup|g0 (t)|.
t≥0

The time derivative of (29) is given by
ġ0 (t) = L−1 {sg0 (s) − g0 (t)|t=0+ }
1p
1 + (tan(φ))2 e−ηt cos(ωd t − φ)
=
m

(39)
(40)

where φ is defined as in the (33).
Setting now ġ0 (t) = 0 in (40) gives
tk =
and since φ +
Therefore

π
2

φ+

π
2

+ kπ
,
ωd

k≥0

> 0, the first maximum is for k = 0.
tNadir =

φ + π2
ωd

which after substituting in (29) gives
||g0 ||∞ = |g0 (tNadir )| =
 



1
π  γ −η
π
−ωη (φ+ π
2)
d
=
cos
φ+
+
sin
φ+
1−e
d+r−1
2
ωd
2
!
r
−1
τr
1
− η φ+ π
1+
e ωd ( 2 )
=
−1
d+r
m
where the last step follows from (28), (30), (33) and
r
(m − dτ )2
ωd
1
−
cos(φ) = q
,
=
2
4mτ r−1
ωd2 + τ1 − η

(41)

as well as several trigonometric identities.
Therefore the nadir of the system frequency is given by
!
r
P
| i ui |
τ r−1 −ωη (φ+ π2)
1
1+
e d
||w̄||∞ = P
(42)
−1
m
i fi d + r
A PPENDIX III
P ROOF OF P ROPOSITION 4
The proof of this proposition requires the following two
lemmas.
Lemma 6: Given a power system under Assumption 1
with gi (s) given by (3), the derivative of φ with respect to
m is given by
∂φ
m + dτ
p
=
−1
∂m
2m 4r mτ − (m − dτ )2
Proof: Notice that while it is not possible to derive a
closed form condition for φ in terms of the system parameters
without the aid of a trigonometric function, it is possible to

This is achieved by first

(m−dτ )
p
−1
4r mτ −(m−dτ )2

∂
∂
tan(φ(m)) =
∂m
∂m
=

∂φ
∂m .

!

2r−1 τ (m + dτ )

3

(4r−1 mτ −(m−dτ )2 ) 2

Now using the fact that φ ∈ (− π2 , π2 ) we get
∂
∂
φ(m) =
arctan (tan(φ(m)))
∂m
∂m


∂
∂
=
arctan(x) x=tan(φ)
tan(φ(m))
∂x
∂m
∂
tan(φ(m))
= cos(φ)2
∂m


(m − dτ )2
∂
= 1−
tan(φ(m))
−1
4r mτ
∂m


(m − dτ )2
2r−1 τ (m + dτ )
= 1−
3
4r−1 mτ
(4r−1 mτ −(m−dτ )2 ) 2
m + dτ
p
=
−1
2m 4r mτ − (m − dτ )2
Lemma 7: Given a power system under Assumption 1
with gi (s) given by (3), the derivative of ωηd with respect
to m is given by
 
∂
2τ (d + r−1 )(m − dτ )
η
=
2 3
∂m ωd
(4mτ r−1 − (m − dτ ) ) 2
Proof: The proof is just by direct computation


 
∂
η
∂ 
m + dτ

q
=
∂m ωd
∂m
2
−1
4mτ r −(m−dτ )
2

=

=
=
=
=

4mτ r−1 −(m−dτ ) −(2τ r−1 −m + dτ )(m + dτ )
q
2
2
(4mτ r−1 −(m−dτ ) )(4mτ r−1 −(m−dτ ) )

4mτ r−1 −2τ r−1 (m + dτ )−(m−dτ )2 + m2 −(dτ )2
(4mτ r−1 −(m−dτ )2 ) 2
−1
−1
2mτ r −2(r τ )(dτ ) + 2mdτ −2(dτ )2
3

2

3

(4mτ r−1 −(m−dτ ) ) 2
2mτ (d + r−1 )−2(dτ )τ (d + r−1 )
(4mτ r−1 −(m−dτ )2 ) 2
2τ (d + r−1 )(m−dτ )
3

2

3

(4mτ r−1 −(m−dτ ) ) 2

We are now ready to prove Proposition 4.
Proof: Using lemmas 6 and 7, and α = π2 + φ we can
compute





∂
1 −η α
d+r−1
1
∂
√
√
||g0 ||∞ =
+ √ e ωd
∂m
∂m
m
τ r−1 
τ r−1
 



η
1
η
∂
∂φ
∂
1
η
−ω α
√
+α
−√
=e d
∂m
∂m ωd
m
m ωd ∂m

!


 !
1
η ∂φ
η
∂
where
=e d
−√
+α
d+r −1
dτ
m ωd ∂m
∂m ωd
ms+
m
.
h(s) =
r
2 + 1 + d s + d+r −1


 
−ωη α 
s
−1
τ
m
mτ
−1
η ∂φ
η
∂
τr e d
1
−2
+α
=
Thus we can compute
d+r−1
m
2
m
ωd ∂m
∂m ωd
√
η




−
α
τ r−1 e ωd
η
−1
m+dτ
g̈0 (t) = L−1 s2 g0 (s) − sg0 (t)|t=0+ − ġ0 (t)|t=0+
p
= √
−2


2 m(d+r−1 )
m
ωd 2m 4r−1 mτ −(m−dτ )2
1
1
!!
(τ − h(s)) −
= L−1
mτ
m
2τ (d+r−1 )(m−dτ )
+α
3
1
2
(4mτ r−1 −(m−dτ ) ) 2
L−1 [h(s)]
=−
mτ
√
η
−
α
η
m+dτ
τ r−1 e ωd
Therefore, it is enough to compute
√
p
=
−1−
2m m(d+r−1 )
ωd 4r−1 mτ −(m−dτ )2
!!
h(t) = L−1 [h(s)]
"
#
2τ (d+r−1 )(m−dτ )
dτ
d+r−1
+2mα
s+
3
2
m
= L−1
m
(4mτ r−1 −(m−dτ ) ) 2
2 + 1 + d s+ d+r−1
s
r
η
τ
m
− α
"
# mτ
τ r−1 e ωd
1
(m+dτ )
dτ
d+r−1
q
=
−1−
s+
d+r−1
m 2m
= L−1 m 2 m 2
4mτ r−1 −(m−dτ )2
(s+η) +ωd
!!
!
−1
d+r−1
m+dτ
2τ (d+r )(m−dτ )
−η dτ
dτ
−ηt
p
m
m
+2mα
2 3
=e
sin(ωd t)
cos(ωd t)+
4r−1 mτ −(m−dτ )2
(4mτ r−1 −(m−dτ ) ) 2
m
ωd
r




−η α
(m+dτ )2
1
τ r−1 e ωd
1
r−1
dτ −ηt
−1−
=
sin(ωd t)
ωd cos(ωd t)+
e
−η+
=
d+r−1
m 2m
4mτ r−1 −(m−dτ )2
mωd
τ
dτ




!
dτ −ηt
1
r−1
4mτ (d+r−1 )(m−dτ )
=
sin(ωd t)
ωd cos(ωd t)+
e
−η+
−α
<0
3
mωd
τ
dτ
(4mτ r−1 −(m−dτ )2 ) 2
dτ −ηt 1
e
(cos(β) cos(ωd t)+sin(β) sin(ωd t))
=
m
cos(β)
dτ −ηt cos(ωd t−β)
=
e
m
cos(β)
A PPENDIX IV
P ROOF OF P ROPOSITION 5
We are now ready to prove Proposition 5
The poof of this proposition requires the computation of
Proof: When t = 0+ , g̈0 (0) < 0. Therefore, the initial
the following lemma.
trend
the g˙0 (t) is decreasing and therefore
Lemma 8: Given g0 (s) as in (3)
1p
1
d −ηt cos(ωd t − β)
g˙0 (0) =
1 + (tan(φ))2 e0 cos(0 − φ) =
g̈0 (t) = − 2 e
m
m
m
cos(β)
is
a
local
maximum.
where the angle β ∈ (− π2 , π2 ), β > φ for φ ∈ (− π2 , π2 )
We will now show that this is indeed a global maximum.
defined by (41) and (33), and β is uniquely defined by
The function g̈0 (t) crosses zero for the first time when ωd t−
ωd
β = π2 or equivalently
cos(β) = r
2

−1
β + π2
ωd2 + τ1 − η + rdτ
(46)
t∗ =
ωd
and
1
r −1
By substituting (46) into (40) and defining
τ − η + dτ
s
.
sin(β) = r
2


2
r −1
1
1
2
r−1
2
ωd + τ − η + dτ
∆ = ωd +
−η +
τ
dτ
Proof: We first compute
we get
τ s2 + s
1
s2 g0 (s) =

(43)
−1
∗
d
1p
mτ s2 + τ1 + m
s + d+r
1 + (tan(φ))2 e−ηt cos(ωd t∗ − φ)
ġ0 (t∗ ) =
mτ
!
m
dτ
d+r −1
π
1
η
π cos(
ms+
m
1
−
2 + (β − φ))
(44)
=
τ−

ωd (β+ 2 )
e
=
d
d+r −1
1
2
mτ
s + τ + m s + mτ
m
cos(φ)
1
1 − ωη (β+ π2 ) sin(β − φ)
=
(τ − h(s))
(45)
=− e d
mτ
m
cos(φ)
1

−ωη α

−21 m−2
m

1 − ωη (β+ π2 ) sin(β) cos(φ) − cos(β) sin(φ)
e d
m
cos(φ)
1 − ωη (β+ π2 )
=− e d
(sin(β) − cos(β) tan(φ))
m



1
r −1
d
−
η
+
dτ
ωd 21 τ1 − m
1 − ωη (β+ π2 )  τ

−
=− e d
m
∆
∆
ωd

 r−1 

d
1 1
1 1
d
−
1 − ωη (β+ π2 )  2 τ − m + dτ
m 
− 2 τ
=− e d
m
∆
∆
!
−1
1 − ωη (β+ π2 ) rdτ
=− e d
m
∆
=−

Finally by simplifying
s

2
1
r−1
2
∆ = ωd +
−η +
τ
dτ
s

2  

2
d + r−1 1 1
r−1
d
1 1 d
=
+
+
−
+
−
mτ
4 τ
m
2 τ m
dτ
s




2
1 1 d r−1
r−1
r−1
+2
−
+
=
mτ
2 τ m dτ
dτ
s


2
r−1
r−1 1
=
+
dτ τ
dτ
we get
ġ0 (t∗ ) = −



1 − ωη (β+ π2 ) 
r
e d

m
r −1

1
dτ τ



r −1
dτ

+



r −1
dτ


1

2  > − m

Therefore, since for any additional instant
of time that
β+ π
2 +kπ
, the factor
ġ0 (t) achieves a local extremum t∗k =
ωd
| cos(ωd t∗k − φ)| does not change, then the maximum is
achieved at t = 0.
R EFERENCES
[1] R. D. Evans and R. C. Bergvall, “Experimental analysis of stability
and power limitations,” Journal of the A.I.E.E., vol. 43, no. 4, pp.
329–340, Apr. 1924.
[2] C. P. Steinmetz, “Power control and stability of electric generating
stations,” Journal of the A.I.E.E., vol. 39, no. 2, pp. 1215–1287, July
1920.
[3] North American Electric Reliability Corporation, “2015 Frequency
Response Annual Analysis,” Tech. Rep., Sept. 2015.
[4] Federal Energy Regulatory Commission and North American Electric
Reliability Corporation, “Arizona-Southern California Outages on
September 8, 2011,” Tech. Rep., Apr. 2012.
[5] Y.-Y. Hsu, S.-W. Shyue, and C.-C. Su, “Low frequency oscillations
in longitudinal power systems: Experience with dynamic stability of
Taiwan power system,” IEEE Trans. on Power Systems, vol. 2, no. 1,
pp. 92–98, Feb. 1987.
[6] M. Klein, G. J. Rogers, and P. Kundur, “A fundamental study of interarea oscillations in power systems,” IEEE Trans. on Power Systems,
vol. 6, no. 3, pp. 914–921, 1991.
[7] M. M. Begovic, “Inter-area oscillations in power systems: A nonlinear
and nonstationary perspective (messina, a.r.) [book reviews],” IEEE
Power and Energy Magazine, vol. 9, no. 2, pp. 76–77, 2011.
[8] N. W. Miller, M. Shao, and S. Venkataraman, “California ISO
(CAISO) Frequency Response Study ,” General Electric, Tech. Rep.,
Nov. 2011.

[9] J. Winkelman, J. Chow, B. Bowler, B. Avramovic, and P. Kokotovic,
“An analysis of interarea dynamics of multi-machine systems,” IEEE
Transactions on Power Apparatus and Systems, vol. PAS-100, no. 2,
pp. 754–763, 1981.
[10] G. Verghese, I. Pérez-Arriaga, and F. Schewppe, “Selective modal
analysis with applications to electric power systems, part ii: the
dynamic stability problem,” IEEE Transactions on Power Apparatus
and Systems, vol. PAS-101, no. 9, pp. 3126–3134, 1982.
[11] E. Tegling, B. Bamieh, and D. Gayme, “The price of synchrony:
Evaluating the resistive losses in synchronizing power networks,” IEEE
Transactions on Control of Network Systems, vol. 2, no. 3, pp. 254–
266, 2015.
[12] A. Mešanović, U. Münz, and C. Heyde, “Comparison of H∞ , H2 , and
pole optimization for power system oscillation damping with remote
renewable generation,” IFAC-PapersOnLine, vol. 49, no. 27, pp. 103–
108, 2016.
[13] E. Mallada, “iDroop: A dynamic droop controller to decouple power
grid’s steady-state and dynamic performance,” in 55th IEEE Conference on Decision and Control (CDC), 12 2016, pp. 4957–4964.
[14] B. K. Poolla, S. Bolognani, and F. Dorfler, “Optimal placement
of virtual inertia in power grids,” IEEE Transactions on Automatic
Control, 2017.
[15] M. Pirani, J. W. Simpson-Porco, and B. Fidan, “System-Theoretic
Performance Metrics for Low-Inertia Stability of Power Networks,”
ArXiv e-prints, Mar. 2017.
[16] M. Andreasson, E. Tegling, H. Sandberg, and K. H. Johansson, “Coherence in Synchronizing Power Networks with Distributed Integral
Control,” ArXiv e-prints, Mar. 2017.
[17] Y. Jiang, R. Pates, and E. Mallada, “Performance tradeoffs of dynamically controlled grid-connected inverters in low inertia power systems,”
in 56th IEEE Conference on Decision and Control (CDC), 2017.
[18] P. Kundur, Power System Stability and Control.
McGraw-Hill
Professional, 1994.
[19] A. Bergen and V. Vittal, Power System Analsysis. Prentice-Hall, 2000.
[20] G. Kou, S. W. Hadley, P. Markham, and Y. Liu, “Developing
Generic Dynamic Models for the 2030 Eastern Interconnection Grid,”
Oak Ridge National Laboratory, Tech. Rep., Dec 2013. [Online].
Available: http://www.osti.gov/scitech/
[21] S. S. Guggilam, C. Zhao, E. Dall’Anese, Y. C. Chen, and S. V. Dhople,
“Engineering inertial and primary-frequency response for distributed
energy resources,” arXiv preprint arXiv:1706.03612, 2017.
[22] D. Apostolopoulou, P. W. Sauer, and A. D. Domı́nguez-Garcı́a,
“Balancing authority area model and its application to the design of
adaptive agc systems,” IEEE Transactions on Power Systems, vol. 31,
no. 5, pp. 3756–3764, 2016.
[23] K. Zhou, J. Doyle, and K. Glover, Robust and Optimal Control.
Prentice Hall, 1996.

