3518 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO.

8, AUGUST 2021

Dynamic Droop Control in Low-Inertia

Power Systems
Yan Jiang , Richard Pates , and Enrique Mallada , Senior Member, IEEE

Abstract—A widely embraced approach to mitigate the I. INTRODUCTION

dynamic degradation in low-inertia power systems is to
HE shift from conventional synchronous generation to
mimic generation response using grid-connected inverters
to restore the stiffness of the grid. In this article, we seek
to challenge this approach and advocate for a principled
T renewable converter-based sources has recently led to a
noticeable degradation of the power system frequency dynamics
design based on a systematic analysis of the performance [1]. At the center of this problem is the reduction of the system-
trade-offs of inverter-based frequency control. With this
aim, we perform a qualitative and quantitative study com- wide inertia that accentuates frequency fluctuations in response
paring the effect of conventional control strategies—droop to disturbances [2], [3]. Besides increasing the risk of frequency
control (DC) and virtual inertia (VI)—on several perfor- instabilities and blackouts [4], this dynamic degradation also
mance metrics induced by L2 and L∞ signal norms. By ex- places limits on the total amount of renewable generation that
tending a recently proposed modal decomposition method, can be sustained by the grid [5]. Ireland, for instance, is already
we capture the effect of step and stochastic power distur-
bances, and frequency measurement noise, on the overall resorting to wind curtailment whenever wind becomes larger
transient and steady-state behavior of the system. Our anal- than 50% of existing demand in order to preserve the grid
ysis unveils several limitations of these solutions, such as stability [6].
the inability of DC to improve dynamic frequency response A widely embraced approach to mitigate this problem is to
without increasing steady-state control effort, or the large mimic synchronous generation response using grid-connected
frequency variance that VI introduces in the presence of
measurement noise. We further propose a novel dynam-i-c converters [7]; that is, to introduce virtual inertia to restore the
droop controller (iDroop) that overcomes the limitations of stiffness that the system used to enjoy [8]. Notable works within
DC and VI. More precisely, we show that iDroop can be this line of research focus on leveraging computational meth-
tuned to achieve high noise rejection, fast system-wide ods [9]–[11] to efficiently allocate synthetic inertial or droop
synchronization, or frequency overshoot (Nadir) elimina- response, or analytical methods that characterize the sensitivity
tion without affecting the steady-state control effort share,
and propose a tuning recommendation that strikes a bal- of different performance metrics to global or spatial variations of
ance among these objectives. Extensive numerical experi- system parameters [12]–[14]. However, to this day, it is unclear
mentation shows that the proposed tuning is effective even whether this particular choice of control is the most suitable
when our proportionality assumption is not valid, and that for the task. On the one hand, unlike synchronous generators
the particular tuning used for Nadir elimination strikes a that leverage stored kinetic energy to modulate electric power
good trade-off among various performance metrics.
injection, converter-based controllers need to actively change
Index Terms—Frequency control, low-inertia power sys- their power injection based on noisy measurements of frequency
tems, static and dynamic performance. or power. On the other hand, converter-based control can be sig-
nificantly faster than conventional generators. Therefore, using
converters to mimic generator behavior does not take advantage
of their full potential. In this article, we seek to challenge this
Manuscript received April 19, 2020; accepted August 14, 2020. Date
of publication October 29, 2020; date of current version July 28, 2021. approach of mimicking generation response and advocate for a
This work was supported in part by ARO under Contract W911NF-17-1- principled control design perspective.
0092, in part by US DoE EERE under Award DE-EE0008006, in part by To achieve this goal, we build on recent efforts by the con-
NSF under Grant CNS 1544771, Grant EPCN 1711188, Grant AMPS
1736448, and CAREER 1752362, in part by the Swedish Foundation trol community on quantifying power network dynamic perfor-
for Strategic Research under Grant SSF RIT15-0091 SoPhy, and in part mance using L2 and L∞ norms [9], [15], and perform a system-
by the Swedish Research Council under Grant VR 2016-04764. Rec- atic study evaluating the effect of different control strategies,
ommended by Associate Editor J. Lavaei. A preliminary version of the
results in this article was presented in IEEE Conferences on Decision such as droop control (DC) [16] and virtual inertia (VI) [17], on
and Control, December 2016 and 2017. (Corresponding author: Enrique a set of static and dynamic figures of merits that are practically
Mallada.) relevant from the power engineering standpoint. More precisely,
Yan Jiang and Enrique Mallada are with the Johns Hopkins University,
Baltimore, MD 21218 USA (e-mail: yjiang@jhu.edu; mallada@jhu.edu). under a mild—yet insightful—proportionality assumption, we
Richard Pates is with the Lund University, SE-221 00 Lund, Sweden compute closed-form solutions and sensitivities of controller
(e-mail: richard.pates@control.lth.se). parameters on the steady-state control effort share, frequency
Color versions of one or more of the figures in this article are available
online at https://ieeexplore.ieee.org. Nadir, L2 -synchronization cost, and frequency variance of the
Digital Object Identifier 10.1109/TAC.2020.3034198 response of a power network to step and stochastic disturbances.

0018-9286 © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See https://www.ieee.org/publications/rights/index.html for more information.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3519

Our analysis unveils the inability of DC and VI to cope with

seemingly opposing objectives, such as synchronization cost
reduction without increasing steady-state effort share (DC),
or frequency Nadir reduction without high frequency variance
(VI). Therefore, rather than clinging to the idea of efficiently
allocating synthetic inertia or droop, we advocate the search of
a better solution.
To this end, we propose novel dynamic droop (iDroop)
control—inspired by classical lead/lag compensation—which
outperforms current control strategies (VI and DC) in an overall
sense. More precisely, the following points are provided.
r Unlike DC that sacrifices steady-state effort share to im-
prove dynamic performance, the added degrees of iDroop
allow to decouple steady-state effort from dynamic per-
formance improvement.
r Unlike VI that amplifies frequency measurement noise, the
lead/lag property of iDroop makes it less sensitive to noise
and power disturbances, as measured by the H2 norm [18] Fig. 1. Block diagram of power network.
of the input–output system defined from measurement
noise and power fluctuations to frequency deviations.
r iDroop can further be tuned to either eliminate the fre-
pin := (pin,i , i ∈ V) ∈ Rn and dp := (dp,i , i ∈ V) ∈ Rn repre-
quency Nadir (i.e., remove the frequency overshoot), by sent power injection set point changes and power fluctuations
compensating for the turbine lag, or to eliminate synchro- around the set point, respectively, and nω := (nω,i , i ∈ V) ∈
nization cost; a feature shown to be unattainable by VI. Rn represents frequency measurement noise. The weighting
All of the above properties are attained through rigorous functions Ŵp (s) and Ŵω (s) can be used to adjust the mag-
analysis on explicit expressions for performance metrics that are nitude of these disturbances in the usual way. The output signal
achieved under a mild yet insightful proportionality assumption ω := (ωi , i ∈ V) ∈ Rn represents the bus frequency deviation
that generalizes prior work [1], [2]. from its nominal value. We now discuss the dynamic elements
We further validate our analysis through extensive numerical in more detail.
simulations, performed on a low-inertia system—the Icelandic 1) Bus Dynamics: The bus dynamics that maps the net
Grid—that does not satisfy our parameter assumptions. Our nu- power bus imbalance uP := (uP,i , i ∈ V) ∈ Rn to the vector
merical results also show that iDroop with the Nadir-eliminated of frequency deviations ω can be described as a feedback loop
tuning designed based on the proportional parameters assump- that comprises a forward-path Ĝ(s) and a feedback-path Ĉ(s),
tion works well even in environments with mixed step and where Ĝ(s) := diag(ĝi (s), i ∈ V) and Ĉ(s) := diag(ĉi (s), i ∈
stochastic disturbances. V) are the transfer function matrices of generators and inverters,
The rest of this article is organized as follows. Section II respectively.
describes the power network model and defines performance a) Generator dynamics: The generator dynamics are com-
metrics. Section III introduces our assumptions and a system posed of the standard swing equations with a turbine, i.e.,
diagonalization that eases the computations and derives some
mi ω̇i = −di ωi + qr,i + qt,i + uP,i (1)
generic results that provide a foundation for further performance
analysis. Section IV analyzes both steady-state and dynamic where mi > 0 denotes the aggregate generator inertia, di > 0
performance of DC and VI, illustrates their limitations, and mo- the aggregate generator damping, qr,i the controllable input
tivates the need for a new control strategy. Section V describes power produced by the grid-connected inverter, and qt,i the
the proposed iDroop and shows how it outperforms DC and change in the mechanical power output of the turbine. The
VI from different perspectives. Section VI validates our results turbine does not react to the frequency deviation ωi until it
through detailed simulations. Section VII concludes the article. exceeds a preset threshold ω ≥ 0, i.e.,
τi q̇t,i = ϕω (ωi ) − qt,i (2)
with ⎧
II. PRELIMINARIES −1
⎪
⎨−rt,i (ωi + ω ) ωi ≤ −ω
A. Power System Model ϕω (ωi ) := 0 −ω < ωi < ω
⎪
⎩ −1
We consider a connected power network composed of n buses −rt,i (ωi − ω ) ωi ≥ ω
indexed by i ∈ V := {1, . . . , n} and transmission lines denoted where τi > 0 represents the turbine time constant and rt,i > 0
by unordered pairs {i, j} ∈ E, where E is a set of 2-element the turbine droop coefficient.
subsets of V. As illustrated by the block diagram in Fig. 1, the Two special cases of our interest are as follows.
system dynamics are modeled as a feedback interconnection Generator Dynamics 1. (Standard swing dynamics): When
of bus dynamics and network dynamics. The input signals |ωi (t)| < ω , the turbines are not triggered and the generator

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3520 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

dynamics can be described by the transfer function standard and well-justified for frequency control on transmission
1 networks [20].
ĝi (s) =
mi s + d i
(3) r Bus voltage magnitudes |Vi |’s are constant; we are not
which is exactly the standard swing dynamics. modeling the dynamics of exciters used for voltage con-
Generator Dynamics 2. (Second-order turbine dynamics): trol; these are assumed to operate at a much faster time-
When ω = 0, the turbines are constantly triggered and the scale.
generator dynamics can be described by the transfer function
r Lines {i, j} are lossless.
r Reactive power flows do not affect bus voltage phase
τi s + 1
ĝi (s) = 2 −1 . (4) angles and frequencies.
mi τi s + (mi + di τi ) s + di + rt,i r Without loss of generality, the equilibrium angle differ-
b) Inverter dynamics: We consider the most commonly
ence (θ0,i − θ0,j ) across each line is less than π/2.
used type of inverters, i.e., the grid-following inverters, which For a first principle derivation of the model, we refer to [21,
adjust their power injection to the network in response to fre- Section VII]. For applications of similar models for frequency
quency deviations. Since power electronics are significantly control within the control literature (see, e.g., [22]–[24]).
faster than the electro-mechanical dynamics of generators, we Remark 2. (Internal stability of (8)): Throughout this arti-
assume that each inverter measures the local grid frequency cle, we consider feedback interconnections of positive real and
deviation ωi and instantaneously updates the output power qr,i . strictly positive real subsystems. Internal stability follows from
Different control laws can be used to map ωi to qr,i . We represent classical results [25]. Since the focus of this article is on perfor-
such laws using a transfer function ĉi (s). The two most common mance, we do not discuss internal stability here in detail. We refer
ones are as follows. to the reader to [26], for a thorough treatment of similar feedback
Inverter Dynamics 1. (Droop Control): This control law can interconnections. From now on, a standing assumption—that
provide additional droop capabilities and is given by can be verified—is that feedback interconnection described in
−1
ĉi (s) = −rr,i (5) Fig. 1 is internally stable.
where rr,i > 0 is the droop coefficient. Remark 3. (Model extensions and limitations): Our model
Inverter Dynamics 2. (Virtual Inertia): Besides providing accepts generalizations accounting for constant power loads via
additional droop capabilities, this control law can compensate Kron reduction [20] and lossy lines with uniform R/X ratios [27],
the loss of inertia and is given by which are not presented due to space constraints. Also, we would
 −1
 like to point out that constant power loads are good surrogates
ĉi (s) = − mv,i s + rr,i (6)
to resistive loads, since the demand variations of resistive loads
where mv,i > 0 is the virtual inertia constant.
are mostly driven by voltage fluctuations which are ignorable in
2) Network Dynamics: The network power fluctuations
the context of frequency control on transmission networks. Ad-
pe := (pe,i , i ∈ V) ∈ Rn are given by a linearized model of the
mittedly, our model cannot capture frequency-dependent loads
power flow equations [19]:
or lossy lines with heterogeneous R/X ratios, which are possible
LB directions of future research.
p̂e (s) = ω̂(s) (7)
s
where p̂e (s) and ω̂(s) denote the Laplace transforms of pe
and ω, respectively.1 The matrix LB is an undirected weighted B. Performance Metrics
Laplacian matrix of the network with elements Having considered the model of the power network, we are
n
now ready to introduce performance metrics used in this article
LB,ij = ∂θj |Vi ||Vj |bij sin(θi − θj ) .
θ=θ0 to compare different inverter control laws.
j=1
1) Steady-State Effort Share: This metric measures the
Here, θ := (θi , i ∈ V) ∈ Rn denotes the angle deviation from its fraction of the power imbalance addressed by inverters, which is
nominal, θ0 := (θ0,i , i ∈ V) ∈ Rn are the equilibrium angles, calculated as the absolute value of the ratio between the inverter
|Vi | is the (constant) voltage magnitude at bus i, and bij is the steady-state output power and the total power imbalance, i.e.,
line {i, j} susceptance.
3) Closed-Loop Dynamics: We will investigate the closed- n
i=1 ĉi (0)ωss,i
loop responses of the system in Fig. 1 from the power injection ES := n +
(9)
i=1 pin,i (0 )
set point changes pin , the power fluctuations around the set
point dp , and frequency measurement noise nω to frequency when the system T̂ωp undergoes a step change in power ex-
deviations ω, which can be described compactly by the transfer citation. Here, ĉi (0) is the dc gain of the inverter and ωss,i
function matrix is the steady-state frequency deviation. A higher steady-state
T̂ (s) := T̂ωp (s) T̂ωdn (s) := T̂ωd (s) T̂ωn (s) . (8) effort share indicates that a larger amount of steady-state power
Remark 1. (Model assumptions): The linearized network output is required from inverters in the process of handling
model (8) implicitly makes the following assumptions which are certain power imbalance. Since this necessary headroom is
usually achieved by either curtailment or additional storage [28],
1 We use hat to distinguish the Laplace transform from its time domain [29], a lower steady-state effort can be associated with lower
counterpart. operational costs and it is therefore desired.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3521

2) Power Fluctuations and Measurement Noise: This

metric measures how the relative intensity of power fluctua-
tions and measurement noise affect the frequency deviations, as
quantified by the H2 norm of the transfer function T̂ωdn
T̂ωdn 2H2
1
∞ ∗
:= 2π −∞ tr(T̂ωdn (jω) T̂ωdn (jω))dω if T̂ωdn is stable
∞ otherwise2 .
(10)
Fig. 2. Equivalent block diagrams of power network under proportion-
The quantity T̂ωdn H2 has several standard interpretations in ality assumption.
terms of the input–output behavior of the system T̂ωdn [18]. In
particular, in the stochastic setting, when the disturbance signals
dp,i and nω,i are independent, zero mean, unit variance, white A. Diagonalization
noise, then limt→∞ E[ω(t)T ω(t)] = T̂ωdn 2H2 .
This means that the sum of the steady-state variances in In order to make the analysis tractable, we require the closed-
the output of T̂ωdn in response to these disturbance equals the loop transfer functions to be diagonalizable. This is ensured by
squared H2 norm of T̂ωdn . Thus, the H2 norm gives a precise the following assumption, which is a generalization of [13], [15].
measure of how the intensity of power fluctuations and measure- Assumption 1. (Proportionality): There exists a proportion-
ment noise affects the frequency deviations of the system. ality matrix F := diag(fi , i ∈ V) ∈ Rn×n≥0 such that
3) Synchronization Cost: This metric measures the size Ĝ(s) = ĝo (s)F −1 and Ĉ(s) = ĉo (s)F
of individual bus deviations from the synchronous response where ĝo (s) and ĉo (s) are called the representative generator
when the system T̂ωp is subject to a step change in power and the representative inverter, respectively.
excitation given by pin = u0 1t≥0 ∈ Rn , where u0 ∈ Rn is a Remark 4. (Proportionality parameters): The parameters fi ’s
given vector direction that allows for power disturbances of represent the individual machine rating. This definition is rather
different magnitudes at individual nodes and 1t≥0 is the unit-step arbitrary for our analysis, provided that Assumption 1 is satis-
function [15]. This is quantified by the squared L2 norm of the fied. Other alternatives could include fi = mi or fi = mi /m
vector of deviations ω̃ := ω − ω̄1n ∈ Rn , i.e., where m is, for example, either the average or maximum gener-
n  ∞
ator inertia. The practical relevance of Assumption 1 is justified,
ω̃22 := ω̃i (t)2 dt . (11)
0 for example, by the empirical values reported in [30], which
i=1
n n show that, at least in regards of order of magnitude, Assumption
Here, ω̄ := ( mi ωi )/(
i=1 mi ) is the system frequency
i=1 1 is a reasonable first-cut approximation to heterogeneity.
that corresponds to the inertia-weighted average of bus fre-
Under Assumption 1, the representative generators of (3) and
quency deviations and 1n ∈ Rn is the vector of all ones.
(4) are given by
4) Nadir: This metric measures the minimum postcontin-
1
gency frequency of a power system, which can be quantified by ĝo (s) = (13)
the L∞ norm of the system frequency ω̄, i.e., ms + d
and
ω̄∞ := max |ω̄(t)| (12) τs + 1
t≥0
ĝo (s) = (14)
when the system T̂ωp has as input a step change in power mτ s2 + (m + dτ ) s + d + rt−1
excitation [15], i.e., pin = u0 1t≥0 ∈ Rn . This quantity matters respectively, with mi = fi m, di = fi d, rt,i = rt /fi , and τi =
in that deeper Nadir increases the risk of underfrequency load τ .3
shedding and cascading outrages. Similarly, the representative inverters of DC (5) and VI (6)
are given by
III. RESULTS ĉo (s) = −rr−1 (15)
In this section, we show that under a simplifying assumption, and
 
it is possible to compute all performance metrics introduced in ĉo (s) = − mv s + rr−1 (16)
Section II-B analytically as functions of the system parameters,
with mv,i = fi mv and rr,i = rr /fi .
which pave us a way to formally compare the conventional
Using Assumption 1, we can derive a diagonalized version
control laws DC and VI in Section IV as well as suggest an
of (8). First, we rewrite Ĝ(s) = F − 2 [ĝo (s)In ]F − 2 and Ĉ(s) =
1 1

improved control law iDroop in Section V. We remark that 1 1

the assumptions are only used in the analysis, but as shown in F 2 [ĉo (s)In ]F 2 , and after a loop transformation obtain Fig. 2(a).
Section VI, the insights and advantages of the proposed solution Then, we define the scaled Laplacian matrix
LF := F − 2 LB F − 2
1 1
are still there when these assumptions do not hold. (17)

2 j represents the imaginary unit which satisfies j 2 = −1 and ω represents 3 We use variables without subscript i to denote parameters of the representa-
the frequency variable. tive generator and inverter.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3522 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

transfer function without making a further declaration in the

rest of this article. For example, we may add “T” if the turbine
is triggered and “DC” if the inverter operates in DC mode as in
ĥp,k,T,DC (s).

B. Generic Results for Performance Metrics

We now derive some important building blocks required for
the performance analysis of the system T̂ described in (21). As
described in Section II-B, the sensitivity to power fluctuations
and measurement noise can be evaluated through the H2 norm
of the system T̂ωdn , while the steady-state effort share, syn-
Fig. 3. Diagonalized block diagram of power network.
chronization cost, and Nadir can all be characterized by a step
response of the system T̂ωp . There are two scenarios that are of
our interest.
by grouping the terms in the upper block of Fig. 2(a). Moreover,
Assumption 2. (Proportional weighting scenario):
since LF ∈ Rn×n is symmetric positive semidefinite, it is real r The noise weighting functions are given by
orthogonally diagonalizable with non-negative eigenvalues [31].
Ŵω (s) = κω F − 2
1 1
Thus, there exists an orthogonal matrix V ∈ Rn×n with V T V = Ŵp (s) = κp F 2 and
V V T = In , such that where κp > 0 and κω > 0 are weighting constants.
LF = V ΛV T (18)
r |ωi (t)| < ω , ∀i ∈ V and t ≥ 0 such that turbines will not
n×n
where Λ := diag(λk , k ∈ {1, . . . , n}) ∈ R≥0 with λk being be triggered.
the kth eigenvalue of LF ordered non-decreasingly (0 = λ1 < Assumption 3. (Step input scenario):
r There is a step change as defined in Section II-B on the
λ2 ≤ · · · ≤ λn )4 and V := ( ni=1 fi )− 2 F 2 1n V⊥ with
1 1

power injection set point, i.e., pin = u0 1t≥0 , dp = 0n , and

V⊥ := v2 ··· vn composed by the eigenvector vk asso- nω = 0n with 0n ∈ Rn being the vector of all zeros.
r ω = 0 such that turbines are constantly triggered.
ciated with λk . Now, applying (17) and (18) to Fig. 2(a) and
5
Remark 5. (Weighting assumption): As a natural counterpart
rearranging blocks of V and V T results in Fig. 2(b). Finally,
of Assumption 1, we look at the case when the power fluctuations
moving the block of ĉo (s)In ahead of the summing junction
and measurement noise are weighted directly and inversely
and combining the two parallel paths produces Fig. 3, where the
proportional to the square root of the bus ratings, respectively.
boxed part is fully diagonalized.
In the case of Ŵp (s), this is equivalent to assuming that demand
Now, by defining the closed-loop with a forward-path ĝo (s)In
fluctuation variances are proportional to the bus ratings, which
and a feedback-path (Λ/s − ĉo (s)In ) as
  is in agreement with the central limit theorem. For Ŵω (s), this
Ĥp (s) = diag ĥp,k (s), k ∈ {1, . . . , n} is equivalent to assuming that the frequency measurement noise
variances are inversely proportional to the bus ratings, which is
where
in line with the inverse relationship between jitter variance and
ĝo (s)
ĥp,k (s) = (19) power consumption for an oscillator in phase-locked-loop [32].
1 + ĝo (s) (λk /s − ĉo (s))
1) Steady-State Effort Share: As indicated by (9), the key
and Ĥω (s) = ĉo (s)Ĥp (s), i.e., of computing the steady-state effort share lies in computing the
 
Ĥω (s) = diag ĥω,k (s), k ∈ {1, . . . , n} steady-state frequency deviation ωss of the system T̂ωp . When
the system synchronizes, the steady-state frequency deviation
where is given by ωss = ωsyn 1n and ωsyn is called the synchronous
ĥω,k (s) = ĉo (s)ĥp,k (s) (20) frequency. In the absence of a secondary control layer, e.g.,
the closed-loop transfer functions from pin , dp , and nω to ω automatic generation control [33], the system can synchronize
become with a nontrivial frequency deviation, i.e., ωsyn = 0.
The following lemma provides a general expression for ωsyn
T̂ωp (s) = F − 2 V Ĥp (s)V T F − 2
1 1
(21a)
in our setting.
T̂ωd (s) = F − 2 V Ĥp (s)V T F − 2 Ŵp (s) Lemma 1. (Synchronous frequency): Let Assumption 3 hold.
1 1
(21b)
If qr,i is determined by a control law ĉi (s), then the output ω
T̂ωn (s) = F V Ĥω (s)V T F Ŵω (s)
− 12 1
2 (21c) of the system T̂ωp synchronizes to the steady-state frequency
respectively. deviation ωss = ωsyn 1n with
n
Note that depending on the specific generator and inverter i=1u0,i
dynamics involved, we may add subscripts in the name of a ωsyn = n  −1
. (22)
i=1 di + rt,i − ĉi (0)
4 Recall Proof: See Appendix A. 
that we assume the power network is connected, which means that
LF has a single eigenvalue at the origin. Now, the theorem below provides an explicit expression for
5 We use k and l to index dynamic modes but i and j to index bus numbers. the steady-state effort share.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3523

Theorem 1. (Steady-state effort share): Let Assumption 3 where

hold. If qr,i is determined by a control law ĉi (s), then the ζ0 := a2 a3 − a1 , ζ1 := a0 a3 , ζ2 := a0 a1 ,
steady-state effort share of the system T̂ωp is given by
n ζ3 := a0 a1 a2 − a20 a3 , ζ4 := −2a0 (a1 b1 b3 + a3 b0 b2 ) .
 i=1 ĉi (0)  (25)
ES = n −1 . (23)
i=1 di + rt,i − ĉi (0)
Otherwise, ĥ2H2 = ∞.
Proof: It follows directly from Lemma 1 that ωss,i = Proof: See Appendix B. 
ωsyn and ni=1 u0,i = ωsyn ni=1 (di + rt,i −1
− ĉi (0)). Plugging Remark 6. (H2 norm of a transfer function lower than fourth-
these two equations to the definition of ES in (9) yields (23).  order): Although Lemma 2 is stated for a fourth-order transfer
2) Power Fluctuations and Measurement Noise: We seek function, it can also be used to find the H2 norm of third-,
to characterize the effect of power fluctuations and frequency second-, and first-order transfer functions by considering appro-
measurement noise on the frequency variance, i.e., the H2 norm priate limits. For example, setting a0 = b0 = and considering
of the system T̂ωdn . the limit → 0, (24) gives the H2 norm of a generic third-order
We first show that the squared H2 norm of T̂ωdn is a weighted transfer function. This process shows that given a stable transfer
sum of the squared H2 norm of each ĥp,k and ĥω,k in the function ĥ(s), if b4 = 0 and
diagonalized system (21). r (third-order) a0 = b0 = 0, then
Theorem 2. (Frequency variance): Define Γ := V T F −1 V . If
a3 b21 + a1 b22 + a1 a2 b23 − 2a1 b1 b3
Assumptions 1 and 2 hold, then ĥ2H2 =
n   2a1 (a2 a3 − a1 )

T̂ωdn 2H2 = Γkk κ2p ĥp,k 2H2 + κ2ω ĥω,k 2H2 . r (second-order) a0 = b0 = a1 = b1 = 0, then
k=1 b2 + a2 b23
Proof: It follows from (8) and (10) that ĥ2H2 = 2
 ∞  2a2 a3
1  r (first-order) a0 = b0 = a1 = b1 = a2 = b2 = 0, then
T̂ωdn 2H2 = tr T̂ωd (jω)∗ T̂ωd (jω) dω
2π −∞ b2
 ∞   ĥ2H2 = 3
1 2a3
+ tr T̂ωn (jω)∗ T̂ωn (jω) dω
2π −∞ otherwise ĥ2H2 = ∞.
=: T̂ωd 2H2 + T̂ωn 2H2 . Remark 7. (Well-definedness by the stability): Note that the
stability of ĥ(s) guarantees that the denominators in all the
We now compute T̂ωd 2H2 . Using (21b) and the fact above H2 norm expressions are nonzero by the Routh–Hurwitz
1
that Ŵp (s) = κp F 2 by Assumption 2, we get T̂ωd (s) = stability criterion.
κp F − 2 V Ĥp (s)V T . Therefore,
1
3) Synchronization Cost: The computation of the synchro-
nization cost defined in (11) for the system T̂ωp in the absence of
T̂ωd (jω)∗ T̂ωd (jω) = κ2p V Ĥp (jω)∗ V T F −1 V Ĥp (jω)V T .
inverter control can be found in [13]. Taking this into account,
Using the cyclic property of the trace, this implies that we can get corresponding results for the system with any control
   
tr T̂ωd (jω)∗ T̂ωd (jω) = κ2p tr Ĥp (jω)∗ ΓĤp (jω) law readily.
Lemma 3. (Synchronization cost): Let Assumptions 1 and 3
where Γ := V T F −1 V . Therefore, it follows that hold. Define ũ0 := V⊥T F − 2 u0 and Γ̃ := V⊥T F −1 V⊥ . Then, the
1
 ∞  
1
T̂ωd 2H2 = κ2p tr Ĥp (jω)∗ ΓĤp (jω) dω synchronization cost of the system T̂ωp is given by
2π −∞  
n  n
ω̃22 = ũT0 Γ̃ ◦ H̃ ũ0
 κ2p Γkk ∞ 2 
= ĥp,k (jω) dω = κp 2
Γkk ĥp,k 2H2 . where ◦ denotes the Hadamard product and H̃ ∈ R(n−1)×(n−1)
2π −∞
k=1 k=1
is the matrix with entries
The result follows from a similar argument on T̂ωn 2H2 .  ∞
H̃kl := hu,k (t)hu,l (t) dt , ∀k, l ∈ {1, . . . , n − 1}
Theorem 2 allows us to compute the H2 norm of T̂ωdn by 0
means of computing the norms of a set of simple scalar transfer with ĥu,k (s) := ĥp,k+1,T (s)/s and ĥp,k,T (s) being a specified
functions. However, for different controllers, the transfer func-
case of the transfer function ĥp,k (s) defined in (19), i.e., when
tions ĥp,k and ĥω,k will change. Since in all the cases these
the turbine is triggered.
transfer functions are of fourth-order or lower, the following
Proof: This is a direct extension of [13, Prop. 2]. 
lemma will suffice for the purpose of our comparison.
Lemma 3 shows that the computation of the synchronization
Lemma 2. (H2 norm of a fourth-order transfer function): Let
cost requires knowing the inner products H̃kl . Yet, the general
b3 s3 + b2 s2 + b1 s + b0 expressions of these inner products for an arbitrary combination
ĥ(s) = 4 + b4
s + a3 s 3 + a2 s 2 + a1 s + a0 of k and l are too tedious to be useful in analysis. Thus, we will
be a stable transfer function. If b4 = 0, then investigate instead bounds on the synchronization cost in terms
ζ0 b20 + ζ1 b21 + ζ2 b22 + ζ3 b23 + ζ4 of the inner products H̃kl when k = l, which are exactly the H2
ĥ2H2 = (24)
2a0 (a1 a2 a3 − a21 − a0 a23 ) norms of transfer functions ĥu,k (s).

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3524 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

Lemma 4. (Bounds for Hadamard product): Let P ∈ Rn×n where the conditions in braces jointly imply ξ > 1.
be a symmetric matrix with minimum and maximum eigenvalues Proof: Basically, Nadir must occur at some non-negative
given by λmin (P ) and λmax (P ), respectively. Then, ∀x, y ∈ Rn finite time instant tnadir , such that ṗu (tnadir ) = 0 and pu (tnadir )
n
 n is a maximum, where pu (t) denotes the unit-step response of
T
  T 
λmin (P ) xk yk ≤ x P ◦ yy
2 2
x ≤ λmax (P ) x2k yk2 . ĥ(s), i.e., p̂u (s) := ĥ(s)/s. We consider three cases based on
k=1 k=1 the value of damping ratio ξ separately.
Proof: See Appendix C.  1) Under damped case (0 ≤ ξ < 1): The output is
Lemma 4 implies the bounds on the synchronization cost. 
Kz 1 s + ξωn
Theorem 3. (Bounds on synchronization cost): Let Assump- p̂u (s) = 2 −
ωn s (s + ξωn )2 + ωd2
tions 1 and 3 hold. Then, the synchronization cost of the system 
T̂ωp is bounded by ω̃22 ≤ ω̃22 ≤ ω̃22 , where ξωn − ωn2 z −1
−
n−1 (s + ξωn )2 + ωd2
k=1 ũ20,k ĥu,k 2H2 
ω̃22 := and with ωd := ωn 1 − ξ 2 , which gives the time domain
maxi∈V (fi )
response
n−1
ũ20,k ĥu,k 2H2 Kz  
ω̃22 :=
k=1
. pu (t) = 2 1 − e−ξωn t η0 sin (ωd t + φ)
mini∈V (fi ) ωn
Proof: By Lemma 3, where
 ∞   2
  ξωn − ωn2 z −1 ωd
ω̃2 =
2
ũT0 Γ̃ ◦ hu (t)hu (t)T ũ0 dt η0 = 1 + and tan φ = .
0 ωd2 ξωn − ωn2 z −1
 ∞ n−1
 Clearly, the above response must have oscillations. There-
≥ λmin (Γ̃) ũ20,k hu,k (t)2 dt fore, for the case 0 ≤ ξ < 1, Nadir always exists.
0 k=1 2) Critically damped case (ξ = 1): The output is
 
n−1
 Kz 1 1 ωn − ωn2 z −1
= λmin (Γ̃) ũ20,k ĥu,k 2H2 p̂u (s) = 2 − −
ωn s s + ωn (s + ωn )2
k=1
which gives the time domain response
n−1
 n−1
ũ20,k ĥu,k 2H2 Kz     
≥ λmin (F −1 ) ũ20,k ĥu,k 2H2 =
k=1 pu (t) = 2 1 − e−ωn t 1 + ωn − ωn2 z −1 t .
maxi∈V (fi ) ωn
k=1
Thus,
which concludes the proof of the lower bound. The first in-   
equality follows from Lemma 4 by setting P = Γ̃, x = ũ0 , ṗu (t) = Kze−ωn t 1 − ωn z −1 t + z −1 .
and y = hu (t) := (hu,k (t), k ∈ {1, . . . , n − 1}) ∈ Rn−1 . The Letting ṗu (t) = 0 yields
     
second inequality follows from the interlacing theorem [31, Th. ωn e−ωn t 1 + ωn − ωn2 z −1 t = e−ωn t ωn − ωn2 z −1
4.3.17]. The proof of the upper bound is similar.  which has a non-negative finite solution
Remark 8. (Synchronization cost in homogeneous case): In
z −1
the system with homogeneous parameters, i.e., F = f In for tnadir =
some f > 0, the identical lower and upper bounds on the syn- ωn z −1 − 1
−1
chronization cost imply that whenever ωn z > 1. For any > 0, it holds that
 
n−1
 ṗu (tnadir − ) = Kze−ωn (tnadir −) ωn z −1 − 1 > 0
ω̃22 = f −1 ũ20,k ĥu,k 2H2 .  
k=1
ṗu (tnadir + ) = Kze−ωn (tnadir +) 1 − ωn z −1 < 0 .
4) Nadir: A deep Nadir poses a threat to the reliable opera- Clearly, Nadir occurs at tnadir . Therefore, for the case
tion of a power system. Hence, one of the goals of inverter control ξ = 1, Nadir is eliminated if and only if ωn z −1 ≤ 1. To
laws is the reduction of Nadir. We seek to evaluate the ability of put it more succinctly, we combine the two conditions
different control laws to eliminate Nadir. To this end, we provide into
a necessary and sufficient condition for Nadir elimination in a 1 = ξ ≤ z/ωn . (27)
second-order system with a zero. 3) Over damped case (ξ > 1): The output is
Theorem 4. (Nadir elimination for a second-order system):  
Kz 1 η1 η2
Assume K > 0, z > 0, ξ ≥ 0, ωn > 0. The step response of a p̂u (s) = 2 − −
ωn s s + σ 1 s + σ2
second-order system with transfer function given by
with
K (s + z)    1 ξ − ωn z −1
ĥ(s) = 2
s + 2ξωn s + ωn2 σ1,2 = ωn ξ ± ξ 2 − 1 and η1,2 = ∓ 
2 2 ξ2 − 1
has no Nadir if and only if
 which gives the time domain response
1 ≤ ξ ≤ z/ωn or
ξ > z/ωn Kz  
ξ ≥ (z/ωn + ωn /z) /2
(26) pu (t) = 2 1 − η1 e−σ1 t − η2 e−σ2 t .
ωn

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3525

Thus, i.e., ωss = ωsyn 1n with

Kz   n
u0,i
ṗu (t) = 2 σ1 η1 e−σ1 t + σ2 η2 e−σ2 t . ωsyn = n  i=1
−1 −1 .
 (29)
ωn i=1 di + rt,i + rr,i
Letting ṗu (t) = 0 yields σ1 η1 e−σ1 t = −σ2 η2 e−σ2 t , Proof: The result follows directly from Lemma 1. 
which has a non-negative finite solution Now, the corollary below gives the expression for the steady-
   state effort share when inverters are under the control law DC
1 1 − ωn z −1 ξ + ξ 2 − 1 or VI.
tnadir =  ln   
2ωn ξ 2 − 1 1 − ωn z −1 ξ − ξ 2 − 1 Corollary 2. (Steady-state effort share of DC and VI): Let
 Assumption 3 hold. If qr,i is under the control law (5) or (6),
whenever 1 − ωn z −1 (ξ − ξ 2 − 1) < 0. For any > 0, then the steady-state effort share of the system T̂ωp,DC or T̂ωp,VI
it holds that is given by
Kz  n −1
ṗu (tnadir − ) > 2 eσ1  σ1 η1 e−σ1 tnadir i=1 rr,i
ωn ES = n  −1

−1 . (30)
 i=1 di + rt,i + rr,i
+σ2 η2 e−σ2 tnadir Proof: The result follows directly from Theorem 1 applied to
(5) and (6). 
= eσ1  ṗu (tnadir ) = 0
Corollary 2 indicates that DC and VI have the same steady-
Kz −σ1   −1
state effort share, which increases as rr,i −1
increase. However, rr,i
ṗu (tnadir + ) < e σ1 η1 e−σ1 tnadir
ωn2 are parameters that also directly affect the dynamic performance
 of the power system, which can be seen clearly from the dynamic
+σ2 η2 e−σ2 tnadir
performance analysis.
= e−σ1  ṗu (tnadir ) = 0
since σ1 > σ2 > 0 and one can show that σ2 η2 < 0. B. Power Fluctuations and Measurement Noise
Clearly, Nadir occurs at tnadir . Therefore, for the case Using Theorem 2 and Lemma 2, it is possible to get closed-
ξ> 1, Nadir is eliminated
 if and only if 1 − ωn z −1 (ξ − form expressions of H2 norms for systems T̂ωdn,DC and T̂ωdn,VI .
ξ 2 − 1) ≥ 0, i.e., ξ 2 − 1 ≥ ξ − z/ωn , which holds Corollary 3. (Frequency variance under DC and VI): Let
if and only if Assumptions 1 and 2 hold. The squared H2 norm of T̂ωdn,DC
 and T̂ωdn,VI is given by
ξ > z/ωn
ξ ≤ z/ωn or n
κ2p + rr−2 κ2ω
ξ ≥ (z/ωn + ωn /z) /2. T̂ωdn,DC H2 =
2
Γkk (31a)
k=1
2mdˇ
Thus, we get the conditions
⎧
⎪ T̂ωdn,VI 2H2=∞ (31b)
⎨ξ > 1
ˇ −1
1 < ξ ≤ z/ωn or ξ > z/ωn (28) respectively, where d := d + rr .
⎪
⎩
ξ ≥ (z/ωn + ωn /z) /2. Proof: We study the two cases separately.
√ We begin with T̂ωdn,DC 2H2 . Applying (13) and (15)
Finally, since ∀a, b ≥ 0, (a + b)/2 ≥ ab with equality only to (19) and (20) shows ĥp,k,DC (s) is a transfer function
when a = b, it follows that the second condition in (28) can only with b4 = a0 = b0 = a1 = b1 = 0, a2 = λk /m, b2 = 0, a3 =
hold when ξ > 1. Thus, we can combine (27) and (28) to yield ˇ
d/m, b3 = 1/m, while ĥω,k,DC (s) is a transfer function
(26). 
with b4 = a0 = b0 = a1 = b1 = 0, a2 = λk /m, b2 = 0, a3 =
ˇ
d/m, b3 = −rr−1 /m. Thus, by Lemma 2
IV. THE NEED FOR A BETTER SOLUTION
1 r−2
We now apply the results in Section III to illustrate the ĥp,k,DC 2H2 = and ĥω,k,DC 2H2 = r .
2mdˇ 2mdˇ
performance limitations of the traditional control laws DC and Then, (31a) follows from Theorem 2.
VI. With this aim, we seek to quantify the frequency variance We now turn to show that T̂ωdn,VI 2H2 is infinite. Define
(10) under DC and VI through the H2 norm of T̂ωdn,DC and m̌ := m + mv . Applying (13) and (16) to (20) yields
T̂ωdn,VI , as well as the steady-state effort share (9), synchro- mv s2 + rr−1 s
nization cost (11), and Nadir (12) through the step response ĥω,k,VI (s) = −
ˇ + λk
m̌s2 + ds
characterizations of T̂ωp,DC and T̂ωp,VI .
which by Lemma 2 has b4 = −mv /m̌ = 0 and thus
ĥω,k,DC 2H2 = ∞. Then, (31b) follows from Theorem 2. 
A. Steady-State Effort Share
Corollary 4. (Optimal rr−1 for T̂ωdn,DC 2H2 ): Let Assump-
Corollary 1. (Synchronous frequency under DC and VI): Let tions 1 and 2 hold. Then
Assumption 3 hold. When qr,i is defined by the control law DC 
(5) or VI (6), the steady-state frequency deviation of the system rr−1 := argminrr−1 >0 T̂ωdn,DC 2H2 = −d + d2 + (κp /κω )2 .
T̂ωp,DC or T̂ωp,VI synchronizes to the synchronous frequency, (32)

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3526 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

Proof: The partial derivative of T̂ωdn,DC 2H2 with respect to Since ĥu,k,VI 2H2 is a function of rr−1 and mv , in what follows
rr−1 is we denote it by ρ(rr−1 , mv ). In order to have an insight on
n
 κ2ω rr−2 + 2dκ2ω rr−1 − κ2p how ĥu,k,VI 2H2 changes with rr−1 and mv , we take partial
∂rr−1 T̂ωdn,DC 2H2 = Γkk . (33) derivatives of ρ(rr−1 , mv ) with respect to rr−1 and mv , i.e.,
k=1
2mdˇ2
∂rr−1 ρ(rr−1 , mv )
By equating (33) to  0, we can solve the corresponding rr−1   2
as rr−1 ± = −d ± d2 + (κp /κω )2 . The only positive root m̌ + τ λk+1 τ + dˇ + λk+1 τ 3 rt−1
 = −    2
is therefore rr−1 := −d + d2 + (κp /κω )2 . We now show 2λk+1 τ dˇ λk+1 τ + dˇ + r−1 + m̌(dˇ + r−1 )
t t
that Γkk > 0, ∀k ∈ {1, . . . , n}. Recall that Γ := V T F −1 V . We
know Γkk = nj=1 (vk,j 2
/fj ). Since vk is an eigenvector, ∀k ∈ ∂mv ρ(rr−1 , mv )
{1, . . . , n}, there must exist at least one j ∈ V such that vk,j = τ 2 rt−1
0. Since fi > 0, ∀i, we have that Γkk > 0, ∀k ∈ {1, . . . , n}. = −    2 .
In addition, since the denominator of (33) is always positive 2 τ dˇ λk+1 τ + dˇ + rt−1 + m̌(dˇ + rt−1 )
and the highest order coefficient of the numerator is positive, Clearly, for all rr−1 ≥ 0, ∂rr−1 ρ(rr−1 , mv ) < 0, which means
whenever 0 < rr−1 < rr−1 , then ∂rr−1 T̂ωdn,DC 2H2 < 0, and if that ρ(rr−1 , mv ) is a monotonically decreasing function of rr−1 .
rr−1 > rr−1 , then ∂rr−1 T̂ωdn,DC 2H2 > 0. Therefore, rr−1 is the Similarly, for all mv ≥ 0, ∂mv ρ(rr−1 , mv ) < 0, which means
minimizer of T̂ωdn,DC 2H2 .  that ρ(rr−1 , mv ) is a monotonically decreasing function of mv .
Two main observations can be made from Corollary 3. First, Therefore, given rr−1 > 0, ∀mv > 0, it holds that
the control parameter rr−1 of DC has an direct effect on the lim ρ(rr−1 , mv ) < ρ(rr−1 , mv ) < ρ(rr−1 , 0) < ρ(0, 0).
mv →∞
size of the frequency variance in the system, which makes it
impossible to require DC to bear an assigned amount of steady- Recall that ĥu,k,VI 2H2 = ρ(rr−1 , mv ), ĥu,k,DC 2H2 =
state effort share and reduce the frequency variance at the same ρ(rr−1 , 0), and ĥu,k,SW 2H2 = ρ(0, 0). The result follows. 
time. The other important point is that VI will induce unbounded Corollary 5. (Comparison of synchronization cost in homoge-
frequency variance, which poses a threat to the operation of the neous case): Denote the synchronization cost of the open-loop
power system. Therefore, neither DC nor VI is a good solution to system as ω̃SW 22 . Then, under Assumptions 1 and 3, given
improve the frequency variance without sacrificing the steady- rr−1 > 0, ∀mv > 0, we can order the synchronization cost when
state effort share. F = f In as
 
n−1 2
k=1 ũ 0,k /λ k+1
C. Synchronization Cost   < ω̃VI 22 < ω̃DC 22 < ω̃SW 22 .
2f d + rt−1
ˇ
Theorem 3 implies that the synchronization cost of T̂ωp,DC
and T̂ωp,VI is bounded by a weighted sum of ĥu,k,DC 2H2 and Proof: The result follows by combining Remark 8 and The-
orem 5. 
ĥu,k,VI 2H2 , respectively. Hence, in order to see the limited Corollary 6. (Lower bound of synchronization cost under DC
ability of DC and VI to reduce the synchronization cost, we need and VI): Under Assumptions 1 and 3, the ordering of the size
to gain a deeper understanding of ĥu,k,DC 2H2 and ĥu,k,VI 2H2 of the bounds on the synchronization cost of open-loop, DC,
first. and VI depends on the parameter values. Thus, we cannot order
Theorem 5. (Bounds of ĥu,k,DC 2H2 and ĥu,k,VI 2H2 ): Let ω̃VI 22 , ω̃DC 22 , and ω̃SW 22 strictly. Instead, we highlight
Assumptions 1 and 3 hold. Then, given rr−1 > 0, ∀mv > 0, that, given rr−1 > 0, the synchronization cost under DC and VI
1 is bounded below by
  < ĥu,k,VI 2H2  
2λk+1 dˇ + rt−1 n−1
ũ 2
/λ
k=1 0,k k+1
 .
< ĥu,k,DC 2H2 < ĥu,k,SW 2H2 2 maxi∈V (fi ) dˇ + rt−1
where ĥu,k,SW 2H2 represents the inner products of the open- Proof: The result follows from Theorems 3 and 5. 
loop system with no additional control from inverters. Corollary 5 provides both upper and lower bounds for the
Proof: Considering that DC can be viewed as VI with synchronization cost under DC and VI in homogeneous case.
mv = 0 and the open-loop system can be viewed as VI The upper bound verifies that DC and VI do reduce the synchro-
with mv = rr−1 = 0, we only compute ĥu,k,VI 2H2 , which nization cost by adding damping and inertia while the lower
straightforwardly implies the other two. Applying (14) and bound indicates that the reduction of the synchronization cost
(16) to (19) shows ĥu,k,VI (s) = ĥp,k+1,T,VI (s)/s is a transfer through DC and VI is limited by certain value that is dependent
function with b4 = a0 = b0 = 0, a1 = λk+1 /(m̌τ ), b1 = on rr−1 . Corollary 6 implies that in the proportional case, the
1/(m̌τ ), a2 = (dˇ + rt−1 + λk+1 τ )/(m̌τ ), b2 = 1/m̌, a3 = synchronization cost under DC and VI is also bounded below
ˇ )/(m̌τ ), b3 = 0. Then, it follows from Lemma 2 that
(m̌ + dτ by a value that is dependent on rr−1 . The fact that the lower
  bound of the synchronization cost under DC and VI is reduced
m̌ + τ λk+1 τ + dˇ as rr−1 increases is unsatisfactory. This is because increasing rr−1
ĥu,k,VI H2 =
2
     .
2λk+1 τ dˇ λk+1 τ + dˇ + rt−1 + m̌ dˇ + rt−1 also increases the steady-state effort share, which is usually not

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3527

desired. However, given a small rr−1 , even if the inertia is very addition of mv makes the tuning region in (35) more accessible,
high, i.e., mv → ∞, the synchronization cost ω̃VI 22 can never which indicates that VI can help a low-inertia power system
reach zero, not to mention ω̃DC 22 . improve Nadir.
We end this section by summarizing the pros and cons of each
D. Nadir controller.
r Droop control: With only one parameter rr−1 , DC can
Finally, with the help of Theorem 4, we can determine the
neither reduce frequency variance or synchronization cost
conditions that the parameters of DC and VI must satisfy to
without affecting steady-state effort share. Moreover, for
eliminate Nadir of the system frequency.
low-inertia systems, DC cannot eliminate Nadir.
Theorem 6. (Nadir elimination under DC and VI): Under r Virtual inertia: VI can use its additional dynamic param-
Assumptions 1 and 3:
r for T̂ωp,DC , the tuning region that eliminates Nadir eter mv to eliminate system Nadir and relatively improve
synchronization cost. However, this comes at the price of
through DC is rr−1 such that
   introducing large frequency variance in response to noise,
rr−1 ≤ m τ −1 − 2 τ −1 rt−1 /m − d (34) and cannot be decoupled from increases in the steady-state
effort share.
r for T̂ωp,VI , the tuning region that eliminates Nadir through
VI is (rr−1 , mv ) such that V. DYNAM-I-C DROOP CONTROL (IDROOP)
  
rr−1 ≤ (m + mv ) τ −1 − 2 τ −1 rt−1 / (m + mv ) − d. We now show how, by moving away from the broadly pro-
posed approach of mimicking generators response, one can
(35) overcome the weaknesses presented in the previous section.
Proof: We start by deriving the Nadir elimination condition With this aim, we introduce an alternative iDroop controller that
for VI. The system frequency of T̂ωp,VI is given by [15] uses dynamic feedback to make a trade-off among the several
n different objectives described in Section II-B. The proposed
u0,i
ω̄VI (t) = i=1
n pu,VI (t) solution is described below.
i=1 fi
Inverter Dynamics 3. (Dynamic Droop Control): The dynam-
where pu,VI (t) is the unit-step response of ĥp,1,T,VI (s). Clearly, ics of an inverter with iDroop is given by the transfer function
as long as pu,VI (t) has no Nadir, neither does ω̄VI (t). Thus, as −1
νi s + δi rr,i
shown later, the core is to apply Theorem 4 to ĥp,1,T,VI (s). ĉi (s) = − (37)
Substituting (14) and (16) to (19) yields s + δi
1 s + τ −1 where δi > 0 and νi > 0 are tunable parameters.
ĥp,1,T,VI (s) = Similarly to (13) and (14), one can define a representative
m̌ s2 + 2ξωn s + ωn2
 iDroop inverter controller as
dˇ + rt−1 τ −1 + d/ˇ m̌
where ωn := , ξ :=  . Now νs + δrr−1
m̌τ ĉo (s) = − (38)
2 (dˇ + rt−1 )/(m̌τ ) s+δ
we are ready to search the Nadir elimination tuning region with νi = fi ν, rr,i = rr /fi , and δi = δ.
by means of Theorem 4. An easy computation shows the In the rest of this section, we expose iDroop to the same
following inequality: 2ξωn − τ −1 = d/ ˇ m̌ < (dˇ + r−1 )/m̌ = performance analysis done for DC and VI in Section IV.
t
2
ωn τ . Equivalently, it holds that ξ < [1/(ωn τ ) + ωn τ ]/2, which
indicates that the second set of conditions in (26) cannot be A. Steady-State Effort Share
satisfied. Hence, we turn to the first set of conditions in (26),
which holds if and only if ξ ≥ 1 and ξωn ≤ τ −1 . Via simple We can show that iDroop is able to preserve the steady-state
algebraic computations, this is equivalent to behavior given by DC and VI.
Corollary 7. (Synchronous frequency under iDroop): Let As-
τ dˇ2 /m̌ − 2dˇ + τ −1 m̌ − 4rt−1 ≥ 0 and d/ ˇ m̌ ≤ τ −1 . (36)
sumption 3 hold. If qr,i is under the control law (37), then
The first condition in (36) can be viewed as a quadratic inequality the steady-state frequency deviation of the system T̂ωp,iDroop
with respect to d,ˇ which holds if and only if
⎛  ⎞ ⎛  ⎞ synchronizes to the synchronous frequency given by (29).
r −1
r −1 Proof: The result follows directly from Lemma 1. 
dˇ ≤ m̌ ⎝τ −1 − 2 ⎠ or dˇ ≥ m̌ ⎝τ −1 + 2 ⎠.
t t
Corollary 8. (Steady-state effort share of iDroop): Let As-
m̌τ m̌τ
sumption 3 hold. If qr,i is under the control law (37), then the
However, only the former region satisfies the second condition steady-state effort share of the system T̂ωp,iDroop is given by
in (36). This concludes the proof of the second statement. The (30).
first statement follows trivially by setting mv = 0.  Proof: The result follows from Theorem 1 applied to (37).
Important inferences can be made from Theorem 6. The fact Corollaries 7 and 8 suggest that iDroop achieves the same
that a small m tends to make the term on the right-hand side synchronous frequency and steady-state effort share as DC and
−1 −1
of (34) negative implies that in a low-inertia power system, it is VI do, which depend on rr,i . Note that besides rr,i , iDroop
impossible to eliminate Nadir using only DC. Undoubtedly, the provides us with two more degrees of freedom by δi and νi .

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3528 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

B. Power Fluctuations and Measurement Noise By Lemma 6, for a given ν, if α1 (ν) < 0, then
The next theorem quantifies the frequency variance under T̂ωdn,iDroop 2H2 always decreases as δ increases. However,
iDroop through the squared H2 norm of the system T̂ωdn,iDroop . according to Lemma 5, even if δ → ∞, we can only ob-
Corollary 9. (Frequency variance under iDroop): Let As- tain T̂ωdn,iDroop 2H2 = T̂ωdn,DC 2H2 . Similarly, if α1 (ν) = 0,
sumptions 1 and 2 hold. The squared H2 norm of T̂ωdn,iDroop then T̂ωdn,iDroop 2H2 keeps constant as δ increases, which
is given by means whatever δ is, we will always obtain T̂ωdn,iDroop 2H2 =
T̂ωdn,iDroop 2H2 T̂ωdn,DC 2H2 . Therefore, iDroop cannot outperform DC when
  α1 (ν) ≤ 0. To put it another way, Lemmas 5 and 6 imply that
n
 (κ2p + rr−2 κ2ω )mδ 2 + (κ2p + ν 2 κ2ω ) dδ
ˇ + λk in order to improve the frequency variance through iDroop, one
= Γkk    .
2m dmδ ˇ 2 + (d + ν) dδ ˇ + λk needs to set ν such that α1 (ν) > 0 and δ as small as practically
k=1
possible. The following lemma characterizes the minimizer ν 
(39)
of T̂ωdn,iDroop 2H2 when δ = 0.
Proof: The proof is based on the Theorem 2 and Lemma 2. Lemma 7. (Minimizer ν  of T̂ωdn,iDroop 2H2 when δ = 0):
Applying (13) and (38) to (19) and (20) shows ĥp,k,iDroop (s) Let Assumptions 1 and 2 hold. Then
is a transfer function with b4 = a0 = b0 = 0, a1 =
(λk δ)/m, b1 = 0, a2 = (dδˇ + λk )/m, b2 = δ/m, a3 = (mδ + ν  := argminδ=0,ν>0 T̂ωdn,iDroop 2H2

d + ν)/m, b3 = 1/m, while ĥω,k,iDroop (s) is a transfer
= − d + d2 + (κp /κω )2 . (40)
function with b4 = a0 = b0 = 0, a1 = (λk δ)/m, b1 = 0, a2 =
ˇ + λk )/m, b2 = −(r−1 δ)/m, a3 = (mδ + d + ν)/m, b3 =
(dδ r
Proof: See Appendix F. 
−ν/m. Thus, by Lemma 2, We are now ready to prove the next theorem.
Theorem 7. (T̂ωdn,iDroop 2H2 optimal tuning): Let Assump-
ˇ + λk
mδ 2 + dδ
ĥp,k,iDroop 2H2 =    tions 1 and 2 hold. Define ν  as in (40). Then,
ˇ 2 + (d + ν) dδ
2m dmδ ˇ + λk r whenever (κp /κω )2 = 2rr−1 d + rr−2 , for any δ > 0 and ν
  such that
rr−2 mδ 2 + ν 2 dδ
ˇ + λk
ĥω,k,iDroop 2H2 =    . ν ∈ [ν  , rr−1 ) or ν ∈ (rr−1 , ν  ] (41)
ˇ 2 + (d + ν) dδ
2m dmδ ˇ + λk
iDroop outperforms DC in terms of frequency variance,
Then, (39) follows from Theorem 2.  i.e.,
The explicit expression of T̂ωdn,iDroop 2H2 given in Corol-
lary 9 is useful to show that iDroop can reduce the frequency T̂ωdn,iDroop 2H2 < T̂ωdn,DC 2H2 .
variance relative to DC and VI. Given the fact that T̂ωdn,VI 2H2 Moreover, the global minimum of T̂ωdn,iDroop 2H2 is
is infinite, the question indeed lies in whether we can find a set obtained by setting δ → 0 and ν → ν  .
of values for parameters δ and ν that ensure T̂ωdn,iDroop 2H2 ≤ r if (κp /κω )2 = 2rr−1 d + rr−2 , then for any δ > 0, by set-
T̂ωdn,DC 2H2 . Fortunately, we can not only find such a set but ting ν → ν  = rr−1 , iDroop matches DC in terms of fre-
also the optimal setting for (39). The following three lemmas set quency variance, i.e.,
the foundation of this important result which is given as Theorem T̂ωdn,iDroop 2H2 = T̂ωdn,DC 2H2 .
7.
Lemma 5. (Limit of T̂ωdn,iDroop 2H2 ): Let Assumptions 1 Proof: As discussed before, to guarantee T̂ωdn,iDroop 2H2 <
and 2 hold. If δ → ∞, then T̂ωdn,iDroop 2H2 = T̂ωdn,DC 2H2 . T̂ωdn,DC 2H2 , one requires to set ν such that α1 (ν) > 0. In
Proof: See Appendix D.  this case, T̂ωdn,iDroop 2H2 always increases as δ increases, so
Lemma 5 shows that T̂ωdn,iDroop 2H2 asymptotically con- choosing δ arbitrarily small is optimal for any fixed ν.
verges to T̂ωdn,DC 2H2 as δ → ∞. The next lemma shows that We now look for the values of ν that satisfy the requirement
this convergence is monotonically from either above or below α1 (ν) > 0. Since the denominator of α1 (ν) is always positive,
depending on the value of the parameter ν. the sign of α1 (ν) only depends on its numerator. Denote the
Lemma 6. (ν-dependent monotonicity of T̂ωdn,iDroop 2H2 numerator of α1 (ν) as Nα1 (ν). Clearly, Nα1 (ν) is a univariate
with respect to δ ): Let Assumptions 1 and 2 hold. Define quadratic function in ν, whose roots are as follows: ν1 = rr−1 and
  ν2 = [(κp /κω )2 − rr−1 d]/d.ˇ Provided that the highest order co-
−dκˇ 2 ν 2 + κ2 + r−2 κ2 ν + dr−2 κ2 − r−1 κ2 efficient of Nα1 (ν) is negative, the graph of Nα1 (ν) is a parabola
ω p r ω r ω r p
α1 (ν) := .
d+ν that opens downwards. Therefore, if ν1 < ν2 , then ν ∈ (ν1 , ν2 )
guarantees α1 (ν) > 0; if ν1 > ν2 , then ν ∈ (ν2 , ν1 ) ∩ (0, ∞)
Then guarantees α1 (ν) > 0. Notably, if ν1 = ν2 , there is no feasible
r T̂ωdn,iDroop 2 is a monotonically increasing or decreas- point of ν to make α1 (ν) > 0.
H2
ing function of δ > 0 if and only if α1 (ν) is positive or The condition ν1 = ν2 happens only if (κp /κω )2 = 2rr−1 d +
negative, respectively. rr , from which it follows that ν  = rr−1 = ν1 = ν2 . Then
−2
r T̂ωdn,iDroop 2 is independent of δ > 0 if and only if α1 (ν  ) = α1 (rr−1 ) = 0. Therefore, by setting ν → ν  = rr−1 ,
H2
α1 (ν) is zero. we get T̂ωdn,iDroop 2H2 = T̂ωdn,DC 2H2 . This concludes the
Proof: See Appendix E.  proof of the second part.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3529

We now focus on the case where the set S = (ν1 , ν2 ) ∪ δ(dˇ + rt−1 + λk+1 τ ) + λk+1 δτ + 1
{(ν2 , ν1 ) ∩ (0, ∞)} is nonempty. Recall from the proof of a1 = , b1 =
mτ mτ
Lemma 6 that T̂ωdn,iDroop H2 = Π(δ, ν). For any fixed ν ∈ S, ˇ −1
δ(m + dτ ) + d + rt + λk+1 τ + ν 1
it holds that α1 (ν) > 0 and thus Π(δ, ν) > Π(0, ν) for any a2 = , b2 =
δ > 0. Recall from the proof of Lemma 7 that ν  is the minimizer mτ m
of Π(0, ν). Hence, (0, ν  ) globally minimizes Π(δ, ν) as long mδτ + m + dτ + ντ
a3 = , b3 = 0 , b4 = 0 .
as ν  ∈ S. In fact, we will show next that ν  is always within S mτ
whenever S = ∅. Considering that a0 → 0 and b0 → 0 as δ → 0 and ν → ∞, we
First, we consider the case when ν1 < ν2 , which implies can employ the H2 norm computation formula for the third-order
−1
2
 (κp /κω ) > 2rr d +
that rr−2 . Then, we have ν  > −d + transfer function in Remark 6. Then
d + 2rr d + rr = rr = ν1 . We also want to show ν  <
2 −1 −2 −1
lim ĥu,k,iDroop 2H2
δ→0,ν→∞
ν2 which holds if and only if
 
  1 2
(κp /κω )2 − rr−1 d (κp /κω )2 + d2 ν
m
1 2
mτ + λmτ
k+1
m
d2 + (κp /κω )2 < +d= = lim = 0.
dˇ dˇ δ→0,ν→∞ 2 λk+1 ( ν ν − λk+1 )
mτ mτ m mτ
which is equivalent to 1 < d2 + (κp /κω )2 /d. ˇ This always
Thus, by Theorem 3, ω̃iDroop 22 = ω̃iDroop 22 = 0, which
holds since (κp /κω )2 > 2rr−1 d + rr−2 . Thus, ν1 < ν  < ν2 .
Similarly, we can prove that in the case when ν1 > ν2 , ν2 < forces ω̃iDroop 22 = 0.
ν  < ν1 holds and thus ν  ∈ (ν2 , ν1 ) ∩ (0, ∞). It follows that Theorem 8 shows that unlike DC and VI that require changes
(0, ν  ) is the global minimizer of Π(δ, ν). on rr−1 to arbitrarily reduce the synchronization cost, iDroop
Finally, by Lemma 5, T̂ωdn,DC 2H2 = Π(∞, ν). The condi- can achieve zero synchronization cost without affecting the
tion (41) actually guarantees ν ∈ S and thus α1 (ν) > 0. Then, steady-state effort share. Naturally, δ ≈ 0 may lead to slow
by Lemma 6, we have T̂ωdn,DC 2H2 = Π(∞, ν) > Π(δ, ν). response and ν → ∞ may hinder robustness. Thus, this result
This concludes the proof of the first part.  should be appreciated from the viewpoint of the additional
Theorem 7 shows that, to optimally improve the frequency tuning flexibility that iDroop provides.
variance, iDroop needs to first set δ arbitrarily close to zero.
Interestingly, this implies that the transfer function ĉo (s) ≈ −ν D. Nadir
except for ĉo (0) = −rr−1 . In other words, iDroop uses its first- Finally, we show that with δ and ν tuned appropriately, iDroop
order lead/lag property to effectively decouple the dc gain ĉo (0) enables the system frequency of T̂ωp,iDroop to evolve as a first-
from the gain at all the other frequencies such that ĉo (jω) ≈ −ν. order response to step power disturbances, which effectively
This decouple is particularly easy to understand in two special makes Nadir disappear. The following theorem summarizes this
regimes: i) If κp  κω , the system is dominated by measure- idea.
ment noise and, therefore, ν  ≈ 0 < rr−1 , which makes iDroop Theorem 9. (Nadir elimination with iDroop): Let Assump-
a lag compensator. Thus, by using lag compensation (setting tions 1 and 3 hold. By setting δ = τ −1 and ν = rr−1 + rt−1 , Nadir
ν < rr−1 ), iDroop can attenuate frequency noise; ii) If κp  κω , (12) of T̂ωp,iDroop disappears.
the system is dominated by power fluctuations and therefore Proof: The system frequency of T̂ωp,iDroop is given by [15]
ν  ≈ κp /κω > rr−1 which makes iDroop a lead compensator. n
u0,i
Thus, by using lead compensation (setting ν > rr−1 ), iDroop ω̄iDroop (t) = i=1 n pu,iDroop (t) (42)
can mitigate power fluctuations. i=1 fi
where pu,iDroop (t) is the unit-step response of ĥp,1,T,iDroop (s).
C. Synchronization Cost If we set δ = τ −1 and ν = rr−1 + rt−1 , then (38) becomes
Theorem 3 implies that the bounds on the synchronization r−1  
ĉo (s) = t − rr−1 + rt−1 . (43)
cost of T̂ωp,iDroop are closely related to ĥu,k,iDroop 2H2 . If we τs + 1
can find a tuning that forces ĥu,k,iDroop 2H2 to be zero, then both Applying (14) and (43) to (19) yields
lower and upper bounds on the synchronization cost converge to 1
ĥp,1,T,iDroop (s) =
zero. Then, the zero synchronization cost is achieved naturally. ms + dˇ + rt−1
The next theorem addresses this problem. whose unit-step response pu,iDroop (t) is a first-order evolution.
Theorem 8. (Zero synchronization cost tuning of iDroop): Let Thus, (42) indicates that Nadir of the system frequency disap-
Assumptions 1 and 3 hold. Then, a zero synchronization cost of pears. 
the system T̂ωp,iDroop , i.e., ω̃iDroop 22 = 0, can be achieved by
setting δ → 0 and ν → ∞. VI. NUMERICAL ILLUSTRATIONS
Proof: Since the key is to show that ĥu,k,iDroop 2H2 → 0 as
δ → 0 and ν → ∞, we can use Lemma 2. Applying (14) and In this section, we present simulation results that compare
(38) to (19) shows ĥu,k,iDroop (s) = ĥp,k+1,T,iDroop (s)/s is a iDroop with DC and VI. The simulations are performed on
transfer function with the Icelandic Power Network taken from the Power Systems
λk+1 δ δ Test Case Archive [34]. Instead of the linear network model
a0 = , b0 = , used in the analysis, the simulations are built upon a nonlinear
mτ mτ

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3530 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

TABLE I
PARAMETERS OF REPRESENTATIVE GENERATOR AND INVERTER

setup including voltage dynamics and nonlinear power flows.

The dynamic model contains 35 generator buses and 83 load
buses. Each of generator buses for simplicity is distinctly in-
dexed by some i ∈ {1, . . . , 35} here. Even though our previous
analysis requires the proportionality assumption (Assumption
1), in the simulations, for every generator bus i, the generator
inertia coefficient, the turbine time constant, and the turbine
droop coefficient are directly obtained from the dataset, i.e.,
mi = md,i , τi = τd,i , and rt,i = rt,d,i .6,7 In addition, turbine
governor deadbands are taken into account such that turbines are
only responsive to frequency deviations exceeding ±0.036Hz.
Given that the values of generator damping coefficients are
not provided by the dataset, we set di = fi d with d being the
representative generator damping coefficient and fi := mi /m
being the proportionality parameters, where m is the repre-
sentative generator inertia defined as the mean of mi ’s, i.e.,
m := ( ni=1 mi )/n. For every load bus, the damping coefficient
is chosen as the mean of all generator damping coefficients.
We refer to this system without inverter control as “SW” in the
simulations.
We then add an inverter to each generator bus i, whose control
law is either one of DC, (filtered) VI, and iDroop.8 The design
of controller parameters will be based on the representative
generator parameters. Hence, besides m and d, we define τ :=
( ni=1 τd,i )/n and rt := ( ni=1 fi )( ni=1 rt,d,i
−1
). Note that to
keep the synchronous frequency unchanged, once inverters are
−1
added, we halve the inverse turbine droop rt,i and assign the
representative inverter droop coefficient rr a value such that the
−1
inverse inverter droop rr,i := fi rr−1 should exactly compensate
this decreased rt−1 in the absence of turbine governor deadbands.
The values of all the representative parameters mentioned above Fig. 4. Comparison between controllers when a −0.3 p.u. step change
are given in Table I. in power injection is introduced to bus number 2 and power fluctuations
Although our current theoretical analysis does not contem- and measurement noise are introduced with κp = 10−3 and κω = 10−4 .
(a) Frequency deviations. (b) Control effort. (c) System frequency and
plate jointly the effects of step and stochastic disturbances, we synchronization cost. (d) Empirical probability distribution function (PDF)
illustrate here that the Nadir-eliminated tuning of Theorem 9 of frequency deviations.
for iDroop can perform quite well in more realistic scenarios
with combined step and stochastic disturbances. In Fig. 4, we
show how different controllers perform when the system is number 2 at time t = 1s as well as power fluctuations and
subject to a step drop of −0.3 p.u. in power injection at bus measurement noise. Since in reality power fluctuations are larger
than measurement noise, we focus on the case dominated by
6 Throughout this section, we use the subscript d, i to denote the original power fluctuations, where κp = 10−3 and κω = 10−4 . As for
parameters of the ith generator bus from the dataset. the representative inverter, we turn δ = τ −1 = 0.218s−1 and
7 For illustrative purpose only, we reassign a part of the droop r
t,d,i ’s on ν = rr−1 + rt−1 = 0.004 srad−1 in iDroop such that Nadir of the
turbines in the dataset to let there be a deeper Nadir in the system frequency.
8 We add a low-pass filter with time constant 0.0001s to VI in simulation since system frequency disappears as suggested by Theorem 9 and we
it is required in reality for a controller with a derivative term. tune mv = 0.0445s2 rad−1 in VI such that the system frequency

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3531

is critically damped. The inverter parameters on each bus i are M ω̇ = − Dω − LB θ + qr + qt + pin (44b)
defined as follows: δi := δ, νi := fi ν, and mv,i = fi mv . where M := diag(mi , i ∈ V) ∈ Rn×n D := diag(di , i ∈
≥0 ,
Some observations are in order. First, all three controllers lead V) ∈ Rn×n , q := (q , i ∈ V) ∈ R , n
and q := (qt,i , i ∈ V) ∈
≥0 r r,i t
to the same synchronous frequency as predicted by Corollaries 1 Rn . In steady state, (44) yields
and 7. Second, although both of VI and iDroop succeed in
eliminating Nadir of the system frequency—which is better LB ωss t = −Dωss − LB θss0 + qr,ss + qt,ss + u0 (45)
than what DC does—the system synchronizes with faster rate where (θss0 + ωss t, ωss , qr,ss , qt,ss ) denotes the steady-state so-
and lower synchronization cost under iDroop than VI. Inter- lution of (44). Equation (45) indicates that LB ωss t is constant
estingly, although our simulation violates the homogeneous and thus LB ωss = 0n . It follows that ωss = ωsyn 1n . Therefore,
assumption, the synchronization costs under VI, DC, and SW (45) becomes
are in ascending order, which aligns with Corollary 5. Third, 0n = − Dωsyn 1n − LB θss0 + qr,ss + qt,ss + u0 (46)
there is no surprise that the frequency variance under VI is n
where qr,ss = (ĉi (0)ωsyn , i ∈ V) ∈ R and qt,ss =
no longer unbounded due to the added filter, yet, for a filter −1
(−rt,i ωsyn , i ∈ V) ∈ Rn when ω = 0 by (2). Premultiplying
with a bandwidth that does not affect the system dynamics, the
system under VI is still more sensitive to noise compared to the (46) by 1Tn and using the property that 1Tn LB = 0Tn , we get the
one under DC or iDroop as illustrated by the empirical PDF. desired result in (22).
Moreover, we highlight the huge control effort required by VI
when compared with DC and iDroop. Last but not least, a bonus B. Proof of Lemma 2
of the Nadir-eliminated tuning is that iDroop outperforms DC First recall that given any state-space realization of ĥ(s), the
in frequency variance as well. This can be understood through H2 norm can be calculated by solving a particular Lyapunov
Theorem 7. Provided that κp  κω , we know from the definition equation. More specifically, suppose
in Lemma 7 that ν  ≈ κp /κω . Thus, for realistic values of sys-  
AB
tem parameters, ν   rr−1 always holds. It follows directly that Σĥ(s) =
CD
ν = rr−1 + rt−1 ∈ (rr−1 , ν  ]. By Theorem 7, iDroop performs
better than DC in terms of frequency variance. Further, the and let X denote the solution to the Lyapunov equation
preceding simulation results suggest that the Nadir-eliminated AX + XAT = −BB T . (47)
tuning of iDroop designed based on the proportional parameters
If ĥ(s) is stable, then
assumption works relatively well even when parameters are 
non-proportional. ∞ if D = 0
ĥ2H2 = T
(48)
CXC otherwise.
VII. CONCLUSION Consider the observable canonical form [35, Ch. 4.1.6] of
This article studies the effect of grid-connected inverter-based ĥ(s) given by
⎡ ⎤
control on the power system performance. When it comes to 0 0 0 −a0 b0
⎢ 1 0 0 −a1 b1 ⎥
the existing two common control strategies, we show that DC ⎢ ⎥
cannot decouple the dynamic performance improvement from Σĥ(s) = ⎢
⎢ 0 1 0 −a2 b2 ⎥ .
⎥ (49)
the steady-state effort share and VI can introduce unbounded fre- ⎣ 0 0 1 −a3 b3 ⎦
quency variance. Thus, we propose a new control strategy named 0 0 0 1 b4
iDroop, which is able to enhance the dynamic performance and Since D = b4 , it is trivial to see from (48) that if b4 = 0,
preserve the steady-state effort share at the same time. We show then ĥ2H2 = ∞. Hence, in the rest of the proof, we assume
that iDroop can be tuned to achieve strong noise rejection, zero b4 = 0. We will now solve the Lyapunov equation analytically
synchronization cost, and frequency Nadir elimination when the for the realization (49). X must be symmetric and thus can be
system parameters satisfy the proportionality assumption. We parameterized as
illustrate numerically that the Nadir-eliminated tuning designed
X = [xij ] ∈ R4×4 , with xij = xji . (50)
based on the proportional parameters assumption strikes a good
T
trade-off among various performance metrics even if parameters Since it is easy to see that CXC = x44 , the problem becomes
are nonproportional. This article illustrates the superiority of solving for x44 . Substituting (49) and (50) into (47) yields the
principled control design when compared with naive approaches following equations:
such as VI. 2a0 x14 = b20 (51a)
x12 − a2 x14 − a0 x34 = − b0 b2 (51b)
APPENDIX
2(x12 − a1 x24 ) = − b21 (51c)
A. Proof of Lemma 1
x23 − a3 x24 + x14 − a1 x44 = − b1 b3 (51d)
Combining (1) and (7) through uP = pin − pe , we get the
(partial) state-space representation of the system T̂ωp as 2(x23 − a2 x34 ) = − b22 (51e)

θ̇ = ω , (44a) 2(x34 − a3 x44 ) = − b23 . (51f)

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
3532 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 66, NO. 8, AUGUST 2021

Through standard algebra, we can solve for x44 , which gives F. Proof of Lemma 7
the right hand side of (24) the denominator is guaranteed to
Recall from the proof of Lemma 6 that T̂ωdn,iDroop 2H2 =
be nonzero by the Routh–Hurwitz criterion. This concludes the
Π(δ, ν). Then, we have
proof. n
κ2p + ν 2 κ2ω 
Π(0, ν) = Γkk
2m(d + ν)
C. Proof of Lemma 4 k=1
whose derivative with respect to ν is given by
First note that n
     κ2ω ν 2 + 2dκ2ω ν − κ2p 
xT P ◦ yy T x = tr P T (x ◦ y) (x ◦ y)T Π (0, ν) = Γkk . (52)
2m(d + ν)2
k=1

= (x ◦ y)T P T (x ◦ y) . Note that (52) and (33) are in the same form. Thus, ν  is
determined in the same way as in the proof of Corollary 4.
Let w := x ◦ y. Since P is symmetric, by Rayleigh [31]
   ACKNOWLEDGMENT
λmin (P )wT w ≤ xT P ◦ yy T x ≤ λmax (P )wT w.
n The authors would like to acknowledge and thank Fernando
Observing that wT w = k=1 x2k yk2 completes the proof. Paganini, Petr Vorobev, and Janusz Bialek for their insightful
comments that helped improve earlier versions of this article.
D. Proof of Lemma 5
REFERENCES
The limit of (39) as δ → ∞ can be computed as
n
[1] B. Kroposki et al., “Achieving a 100% renewable grid: Operating electric
 κ2p + rr−2 κ2ω power systems with extremely high levels of variable renewable energy,”
lim T̂ωdn,iDroop 2H2 = Γkk = T̂ωdn,DC 2H2 IEEE Power Energy Mag., vol. 15, no. 2, pp. 61–73, Mar. 2017.
δ→∞
k=1
2mdˇ [2] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, “Foundations and
challenges of low-inertia systems (invited paper),” in Proc. Power Syst.
where the second equality follows from (31a). Comput. Conf., Jun. 2018, pp. 1–25.
[3] T. Ackermann, T. Prevost, V. Vittal, A. J. Roscoe, J. Matevosyan, and N.
Miller, “Paving the way: A future without inertia is closer than you think,”
E. Proof of Lemma 6 IEEE Power Energy Mag., vol. 15, no. 6, pp. 61–69, Nov. 2017.
[4] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of low rotational
Provided that T̂ωdn,iDroop 2H2 is a function of δ and ν, in inertia on power system stability and operation,” in Proc. IFAC World
what follows, we denote it by Π(δ, ν). To make it clear how Congr., Aug. 2014, pp. 7290–7297.
[5] A. S. Ahmadyar, S. Riaz, G. Verbič, A. Chapman, and D. J. Hill, “A frame-
Π(δ, ν) changes with δ, we first put it into the equivalent form work for assessing renewable integration limits with respect to frequency
of performance,” IEEE Trans. Power Syst., vol. 33, no. 4, pp. 4444–4453,

n   Jul. 2018.
α1 (ν)δ 2 [6] J. O’Sullivan, A. Rogers, D. Flynn, P. Smith, A. Mullane, and M. O’Malley,
Π(δ, ν) = Γkk + α 5 (ν) “Studying the maximum instantaneous non-synchronous generation in an
α2 δ 2 + α3 (ν)δ + α4 (ν, λk ) island system—Frequency stability challenges in Ireland,” IEEE Trans.
k=1
Power Syst., vol. 29, no. 6, pp. 2943–2951, Nov. 2014.
with [7] B. K. Bose, “Global energy scenario and impact of power electronics in
  21st century,” IEEE Trans. Ind. Electron., vol. 60, no. 7, pp. 2638–2651,
ˇ 2 ν 2 + κ2 + r−2 κ2 ν + dr−2 κ2 − r−1 κ2
−dκ Jul. 2013.
ω p r ω r ω r p
α1 (ν) := [8] R. Ofir, U. Markovic, P. Aristidou, and G. Hug, “Droop vs. virtual inertia:
d+ν Comparison from the perspective of converter operation mode,” in Proc.
α2 := 2mdˇ , α3 (ν) := 2(d + ν)dˇ IEEE Int. Energy Conf., Jun. 2018, pp. 1–6.
[9] B. K. Poolla, D. Groß, and F. Dörfler, “Placement and implementation
of grid-forming and grid-following virtual inertia and fast frequency re-
κ2p + ν 2 κ2ω
α4 (ν, λk ) := 2(d + ν)λk , α5 (ν) := . sponse,” IEEE Trans. Power Syst., vol. 34, no. 4, pp. 3035–3046, Jul. 2019.
2m(d + ν) [10] S. S. Guggilam, C. Zhao, E. Dall’Anese, Y. C. Chen, and S. V. Dhople, “Op-
timizing DER participation in inertial and primary-frequency response,”
We then take the partial derivative of Π(δ, ν) with respect to δ IEEE Trans. Power Syst., vol. 33, no. 5, pp. 5194–5205, Sep. 2018.
[11] U. Markovic, Z. Chu, P. Aristidou, and G. Hug, “LQR-based adaptive
as virtual synchronous machine for power systems with high inverter pen-
n
   etration,” IEEE Trans. Sustain. Energy, vol. 10, no. 3, pp. 1501–1512,
α3 (ν)δ 2 + 2α4 (ν, λk )δ Jul. 2019.
∂δ Π(δ, ν) = α1 (ν) Γkk .
(α2 δ 2 + α3 (ν)δ + α4 (ν, λk ))2 [12] L. Guo, C. Zhao, and S. H. Low, “Graph Laplacian spectrum and primary
k=1 frequency regulation,” in Proc. IEEE Conf. Decis. Control, Dec. 2018,
pp. 158–165.
Since m > 0, d > 0, ν > 0, and rr−1 > 0, α2 and α3 (ν) [13] F. Paganini and E. Mallada, “Global analysis of synchronization perfor-
are positive. Also, given that all the eigenvalues of the scaled mance for power systems: Bridging the theory–practice gap,” IEEE Trans.
Laplacian matrix LF are non-negative, α4 (ν, λk ) must be non- Autom. Control, vol. 65, no. 7, pp. 3007–3022, Jul. 2020.
[14] L. Pagnier and P. Jacquod, “Optimal placement of inertia and primary
negative. Thus, ∀δ > 0, (α3 (ν)δ 2 + 2α4 (ν, λk )δ)/(α2 δ 2 + control: A matrix perturbation theory approach,” IEEE Access, vol. 7,
α3 (ν)δ + α4 (ν, λk ))2 > 0. pp. 145889–145900, 2019.
Recall from the proof of Corollary 4 that Γkk > [15] F. Paganini and E. Mallada, “Global performance metrics for synchro-
nization of heterogeneously rated power systems: The role of machine
0, ∀k ∈ {1, . . . , n}. Therefore, ∀δ > 0, sign(∂δ Π(δ, ν)) = models and inertia,” in Proc. Allerton Conf. Commun. Control, Comput.,
sign(α1 (ν)). Oct. 2017, pp. 324–331.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.
JIANG et al.: DYNAMIC DROOP CONTROL IN LOW-INERTIA POWER SYSTEMS 3533

[16] K. D. Brabandere, B. Bolsens, J. V. den Keybus, A. Woyte, J. Driesen, and Yan Jiang received the B.Eng. degree in electri-
R. Belmans, “A voltage and frequency droop control method for parallel cal engineering and automation from the Harbin
inverters,” IEEE Trans. Power Electron., vol. 22, no. 4, pp. 1107–1115, Institute of Technology, China, in 2013, and the
Jul. 2007. M.S. degree in electrical engineering from the
[17] H. Beck and R. Hesse, “Virtual synchronous machine,” in Proc. Int. Conf. Huazhong University of Science and Technol-
Elect. Power Qual. Utilisation, Oct. 2007, pp. 1–6. ogy, in 2016. She is currently working toward
[18] E. Tegling, B. Bamieh, and D. F. Gayme, “The price of synchrony: the Ph.D. degree at the Department of Electrical
Evaluating the resistive losses in synchronizing power networks,” IEEE and Computer Engineering and the M.S.E. de-
Trans. Control Netw. Syst., vol. 2, no. 3, pp. 254–266, Sep. 2015. gree at the Department of Applied Mathematics
[19] K. Purchala, L. Meeus, D. Van Dommelen, and R. Belmans, “Usefulness and Statistics, Johns Hopkins University, Balti-
of DC power flow for active power flow analysis,” in Proc. IEEE Power more, MD, USA.
Eng. Soc. Gen. Meeting, Jun. 2005, pp. 454–459. Her research interests include the area of control of power systems.
[20] P. Kundur, Power System Stability and Control. New York, NY, USA:
McGraw-Hill, 1994.
[21] C. Zhao, U. Topcu, N. Li, and S. H. Low, “Power system dynam-
ics as primal-dual algorithm for optimal load control,” May 2013,
arXiv:1305.0585v1.
[22] C. Zhao, U. Topcu, N. Li, and S. Low, “Design and stability of load- Richard Pates received the bachelors in me-
side primary frequency control in power systems,” IEEE Trans. Automat. chanical engineering, the M.Eng. degree in
Control, vol. 59, no. 5, pp. 1177–1189, May 2014. computer engineering, and the Ph.D. degree in
[23] N. Li, C. Zhao, and L. Chen, “Connecting automatic generation control control engineering from the University of Cam-
and economic dispatch from an optimization view,” IEEE Trans. Control bridge, U.K, 2009 and 2014, respectively.
Netw. Syst., vol. 3, no. 3, pp. 254–264, Sep. 2016. He is currently an Associate Professor with
[24] E. Mallada, C. Zhao, and S. Low, “Optimal load-side control for frequency Lund University, Lund, Sweden. His research
regulation in smart grids,” IEEE Trans. Autom. Control, vol. 62, no. 12, interests include modular methods for control
pp. 6294–6309, Dec. 2017. system design, stability and control of electrical
[25] H. K. Khalil, Nonlinear Systems, 3rd ed. Englewood Cliffs, NJ, USA: power systems, and fundamental performance
Prentice Hall, 2002. limitations in large-scale systems.
[26] R. Pates and E. Mallada, “Robust scale-free synthesis for frequency
control in power systems,” IEEE Trans. Control Netw. Syst., vol. 6, no. 3,
pp. 1174–1184, Sep. 2019.
[27] Young-Hyun Moon, Byoung-Hoon Cho, Tae-Hoon Rho, and Byoung-Kon
Choi, “The development of equivalent system technique for deriving an
energy function reflecting transfer conductances,” IEEE Trans. Power
Syst., vol. 14, no. 4, pp. 1335–1341, Nov. 1999. Enrique Mallada (Senior Member, IEEE) re-
[28] L. Bird et al., “Wind and solar energy curtailment: A review of interna- ceived the Ingeniero en Telecomunicaciones
tional experience,” Renewable Sustain. Energy Rev., vol. 65, pp. 577–586, degree from Universidad ORT, Uruguay, in
Nov. 2016. 2005, and the Ph.D. degree in electrical and
[29] V. Knap, S. K. Chaudhary, D. Stroe, M. Swierczynski, B. Craciun, and R. computer engineering with a minor in applied
Teodorescu, “Sizing of an energy storage system for grid inertial response mathematics from Cornell University, USA, in
and primary frequency reserve,” IEEE Trans. Power Syst., vol. 31, no. 5, 2014.
pp. 3447–3456, Sep. 2016. Since 2016, he has been an Assistant Pro-
[30] G. Kou, S. W. Hadley, P. Markham, and Y. Liu, “Developing generic fessor of Electrical and Computer Engineering
dynamic models for the 2030 eastern interconnection grid,” Oak Ridge with the Johns Hopkins University, Baltimore,
National Laboratory, Tech. Rep., Dec. 2013. MD, USA. He was a Postdoctoral Fellow with the
[31] R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. Cambridge, U.K.: Center for the Mathematics of Information at Caltech from 2014 to 2016.
Cambridge University Press, 2012. Dr. Mallada received the NSF CAREER Award in 2018, the ECE
[32] T. C. Weigandt, B. Kim, and P. R. Gray, “Analysis of timing jitter in Director’s Ph.D. Thesis Research Award for his dissertation in 2014, the
CMOS ring oscillators,” in Proc. IEEE Int. Symp. Circuits Syst., May 1994, Center for the Mathematics of Information (CMI) Fellowship from Caltech
pp. 27–30. in 2014, and the Cornell University Jacobs Fellowship in 2011. His
[33] F. P. deMello, R. J. Mills, and W. F. B’Rells, “Automatic generation research interests include control, dynamical systems, and optimization,
control part II—Digital control techniques,” IEEE Trans. Power App. Syst., with applications to engineering networks such as power systems and
vol. PAS-92, no. 2, pp. 716–724, Mar. 1973. the Internet.
[34] University of Edinburgh. Power systems test case archive. Jan. 2011.
[Online]. Available: https://www.maths.ed.ac.uk/optenergy/Network
Data/icelandDyn/
[35] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Anal-
ysis and Design, 2nd ed. Hoboken, NJ, USA: Wiley, 2005.

Authorized licensed use limited to: Johns Hopkins University. Downloaded on July 31,2021 at 23:20:16 UTC from IEEE Xplore. Restrictions apply.