JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

1

Distributed Optimal Frequency Control Considering
a Nonlinear Network-Preserving Model

arXiv:1709.01543v3 [] 13 Feb 2018

Zhaojian Wang, Feng Liu, John Z. F. Pang, Steven Low, Fellow, IEEE and Shengwei Mei, Fellow, IEEE

Abstract—This paper addresses the distributed optimal frequency control of power systems considering a networkpreserving model with nonlinear power flows and excitation
voltage dynamics. Salient features of the proposed distributed
control strategy are fourfold: i) nonlinearity is considered to
cope with large disturbances; ii) only a part of generators
are controllable; iii) no load measurement is required; iv)
communication connectivity is required only for the controllable
generators. To this end, benefiting from the concept of “virtual
load demand”, we first design the distributed controller for the
controllable generators by leveraging the primal-dual decomposition technique. We then propose a method to estimate the
virtual load demand of each controllable generator based on local
frequencies. We derive incremental passivity conditions for the
uncontrollable generators. Finally, we prove that the closed-loop
system is asymptotically stable and its equilibrium attains the
optimal solution to the associated economic dispatch problem.
Simulations, including small and large-disturbance scenarios,
are carried on the New England system, demonstrating the
effectiveness of our design.
Index Terms—Frequency control, network-preserving model,
distributed control, incremental output passivity

I. I NTRODUCTION
Frequency restoration and economic dispatch (ED) are two
important problems in power system operation. Conventionally, they are implemented hierarchically in a centralized fashion, where the former is addressed in a fast time scale while
the latter in a slow time scale [1], [2]. While this centralized
hierarchy works well for the traditional power system, it
may not be able to keep pace with the fast development of
our power system due to: 1) slow response, 2) insufficient
flexibility, 3) low privacy, 4) intense communication, and 5)
single point of failure issue. In this regard, the idea of breaking
such a hierarchy is proposed in [3], [4].
In [5], an intrinsic connection between the asymptotic
stability of the dynamical frequency control system with the
ED problem is proposed. It leads to a so-called inverseengineering methodology for designing optimal frequency
controllers, where the (partial) primal-dual gradient algorithm
This work was supported by the National Natural Science Foundation of
China (No. 51677100, U1766206, No. 51621065), the US National Science
Foundation through awards EPCN 1619352, CCF 1637598, CNS 1545096,
ARPA-E award de-ar0000699, and Skoltech through Collaboration Agreement
1075-MRA. (Corresponding author: Shengwei Mei)
Z. Wang, F. Liu and S. Mei are with the State Key Laboratory of
Power System and Department of Electrical Engineering, Tsinghua University,
Beijing, China, 100084. e-mail: (meishengwei@mail.tsinghua.edu.cn).
S. H. Low and J. Pang are with the Department of Electrical Engineering, California Institute of Technology, Pasadena, CA, USA, 91105. email:(slow@caltech.edu).

plays an essential role [1], [6]. When designing optimal frequency controllers, in choice of power flow models, including
the linear model (usually associated with DC power flow, e.g.
[5], [7]–[15]) and the nonlinear model (usually associated with
AC power flow, e.g. [16]–[21]), is crucial. While the closedloop system can be interpreted in a linear model as carrying
out a primal-dual algorithm for solving ED, this interpretation
of frequency control may not hold in a nonlinear model. In
addition to nonlinear power flow, excitation voltage dynamics
are considered in [16]–[18], making the model more realistic.
The aforementioned idea is further developed to enable the
design of distributed optimal frequency controllers. Roughly
speaking, the works of distributed optimal frequency control
can be divided into two categories in terms of different power
system models: network-reduced models e.g. [5], [12]–[16],
[18] and network-preserving models e.g. [7]–[11], [17]. In
network-reduced models, generators and/or loads are aggregated and treated as one bus or control area, which are connected to each other through tie lines. In [5], [13], aggregated
generators in each area are driven by automatic generation
control (AGC) to restore system frequency. [12], [14]–[16],
[18] further consider both the aggregated generators and load
demands in frequency control. In network-preserving models,
generator and load buses are separately handled with different
dynamic models and coupled by power flows, rendering a
set of differential algebraic equations (DAEs). In [7], an
optimal load control (OLC) problem is formulated and a
primary load-side control is derived as a partial primal-dual
gradient algorithm for solving the OLC problem. The design
approach is extended to secondary frequency control (SFC)
that restores nominal frequency in [8]. It is generalized in [9],
where passivity condition guaranteeing stability is proposed
for each local bus. Then, a unified framework combining
load and generator control is advocated in [10]. A similar
model is also utilized in [11], where only limited control
coverage is needed. Similar to [18], the Hamiltonian method
is used to analyze the network-preserving model in [17].
Compared with the network-reduced model, the networkpreserving model describes power systems more precisely
and appear more suitable for analyzing interactions between
different control areas. Therefore, we specifically consider the
network-preserving model in this work.
Almost all of aforementioned works assume that all buses
are controllable and load demands at all buses are accurately measurable, especially for those proposing secondary
frequency control. Moreover, it is usually assumed that the
communication network has the same topology of the power
grid. These assumptions are strong and arguably unrealistic

di  d̂i

Pei

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

II. P OWER S YSTEM M ODELS

Pi

g

i

 i , i , Eqi'
Pei , Qei

p j  1 q j

i ,  i ,Vi

E fi

pi  1 qi

Fig. 1. Summary of notations

buses and pure load buses. Denote controllable generator buses
as NCG = {1, 2, . . . , nCG }, uncontrollable generator buses as
NU G = {nCG + 1, nCG + 2, . . . , nCG + nU G }, and pure load
buses as NL = {nCG +nU G +1, . . . , nCG +nU G +nL }. Then
the set of generator buses is NG = NCG ∪ NU G and set of
all the buses is N = NG ∪ NL . It should be noted that load
can be connected to any bus besides pure load buses.
Let E ⊆ N × N be the set of lines, where (i, j) ∈ E if
buses i and j are connected directly. Then the whole system is
modeled as a connected graph
G = (N , E). The admittance of
√
each line is Yij := Gij + −1Bij with Gij = 0 for every line.
Denote the bus voltage by Vi ∠θi , where Vi is the amplitude
and θi is the voltage phase angle. The active and reactive
power Pij , Qij from bus i to bus j is
Pij = Vi Vj Bij sin (θi − θj )
Qij =

Bij Vi2

(1a)

− Vi Vj Bij cos (θi − θj )

(1b)

For convenience, most notations are summarized in Fig.1.
B. Synchronous Generators
For i ∈ NG , we use the standard third-order generator
model (e.g. [18], [22], [23]) (2a)-(2c). Here (2d) is the simplified governor-turbine model, and (2e) is the excitation voltage
control model:
δ̇i = ωi
ω̇i =

(Pig

(2a)
− Di ωi − Pei )/Mi

0

0

(2b)

0

Ėqi = −Eqi /Td0i + Ef i /Td0i

(2c)

Ṗig

(2d)

=

−Pig /Ti

+

ugi

Ėf i = h(Ef i , Eqi )

(2e)

In this model, Mi is the moment of inertia; Di the damping
0
constant; Td0i the d-axis transient time constant; Ti the time
constant of turbine; δi the power angle of generator i; ωi the
generator frequency deviation compared to an steady state
value; Pig the mechanical power input; pj the active load
0
demand; Pei the active power injected to network; Eqi the qaxis transient internal voltage; Eqi the q-axis internal voltage;
Ef i the excitation voltage. Eqi is given by

A. Power Network
Power systems are composed of many generators and loads,
which are integrated in different buses and connected by power
lines, forming a power network. Buses are divided into three
types, controllable generator buses, uncontrollable generator

 j ,  j ,V j

j
Pij , Qij

…

for practice. First, only a part of generators and controllable
load buses can participate in frequency control in practice.
Second, the communication network may not be identical
to the topology of the power grid. Third, it is difficult to
accurately measure the load injection on every bus. In extreme
cases, even the number of load buses is unknown in practice.
These problems highlight some reasons why theoretical work
in this area is hardly applied in practical systems.
In this context, a novel distributed frequency recovery
controller is proposed that only needs to be implemented on
controllable generator buses. To this end, a network-preserving
power system model is adopted. This work is an extension of
our former work [11]. However, in this paper, we design a
totally different controller considering a third-order nonlinear
generator model with excitation voltage dynamics and nonlinear power flow. It is also motivated partly by [17], which
adopts a similar model, although our results are significantly
different from those in [17]. Differing from [11], the controller
avoids load measurement, which greatly facilitates implementation. By using LaSalle’s invariance principle, it is proved that
the closed-loop system converges to an equilibrium point that
solves the economic dispatch problem. The salient features of
the proposed distributed optimal frequency controller are:
1) Model: The network-preserving model of power system
is used, including excitation voltage dynamics and nonlinear power flow. This, unlike work on the linear model,
returns a valid controller even under large disturbances.
2) Controllability: We allow an arbitrary subset of generator buses to be controllable.
3) Controller: The distributed controller achieves the optimal solution to economic dispatch while restoring
the nominal frequency, provided that certain sufficient
conditions on active power dynamics of uncontrollable
generators and excitation voltage dynamics of all generators are satisfied.
4) Communication: Communication is required between
neighboring controllable generators only, and the communication network can be arbitrary as long as it remains connected.
5) Measurement: No load measurement is needed, and the
controller is adaptive to unknown load changes.
The rest of this paper is organized as follows. In Section
II, we introduce the power system model. Section III formulates the optimal economic dispatch problem. The distributed
controller is proposed in Section IV, and we further prove
the optimality and stability of the corresponding equilibrium
point in Section V. The load estimation method is proposed
in Sectin VI. We confirm the performance of controllers via
simulations on a detailed power system in Section VII. Section
VIII concludes the paper.

2

0

Eqi =

xdi − xdi
xdi 0
Vi cos(δi − θi )
0 Eqi −
0
xdi
xdi

(3)
0

where xdi is the d-axis synchronous reactance, and xdi is the
d-axis transient reactance.

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

3

The active and reactive power (denoted by Qei ) injection to
the network are

III. F ORMULATION OF O PTIMAL F REQUENCY C ONTROL
A. Optimal Power-Sharing Problem in Frequency Control

0

Eqi Vi
Pei = 0 sin(δi − θi )
xdi

(4a)

0

Eqi Vi
Vi2
− 0 cos(δi − θi )
0
xdi
xdi

Qei =

(4b)

For controllable generators i ∈ NCG , the capacity limits are
g

P gi ≤ Pig ≤ P i
where

g
P gi , P i

are lower and upper limits of

(5)

The purpose of optimal frequency control is to let all the
controllable generators share power mismatch economically
when restoring frequency. Then we have the following optimization formulation, denoted by SFC.
X
fi (Pig )
(8a)
SFC:
min
g
i∈NCG
Pi ,i∈NCG
X
X
X
(8b)
Pig∗
pi −
Pig =
s.t.
i∈NCG
P gi ≤

Pig .

C. Dynamics of Voltage Phase Angles
To build a network-preserving power system model, relation
between generators and power network should be explicitly
established. In this paper, loads for bus i ∈ N are simply
modeled as constant active and reactive power injections. Then
the following equations are used to dictate the power balance
and voltage phase-angle dynamics at each bus:
θ̇i = ω̃i , i ∈ N

(6a)
X

0 = Pei − D̃i ω̃i − pi −
Pij , i ∈ NG
j∈Ni
X
0 = −D̃i ω̃i − pi −
Pij , i ∈ NL
j∈Ni
X
0 = Qei − qi −
Qij , i ∈ NG
j∈Ni
X
0 = −qi −
Qij , i ∈ NL
j∈Ni

(6b)
(6c)
(6d)
(6e)

where, pi , qi are active and reactive load demands, respectively; ω̃i the frequency deviation at bus i; Ni the set of buses
connected directly to bus i; D̃i the damping constant at bus i;
D̃i ω̃i the change of frequency-sensitive load [7].
In power system, line power flows are mainly related to
power angle difference between two buses rather than the
power angles independently. Then, we define new variables
to denote angle differences as ηii := δi − θi , i ∈ NG and
ηij := θi − θj , i, j ∈ N . The time derivative of ηii and ηij are
η̇ii = ωi − ω̃i , i ∈ NG

(7a)

η̇ij = ω̃i − ω̃j , i, j ∈ N , i 6= j

(7b)

respectively. In the following analysis, we use ηii and ηij as
new state variables instead of δi and θi .
To summarize, (1) − (4), (6b) − (6e), (7a), (7b) constitute the
nonlinear network preserving model of power systems, which
is in a form of differential-algebraic equations (DAE).
D. Communication Network
In this paper, we consider a communication graph among
the buses of controllable generators only. Denote E ⊆
NCG × NCG as the set of communication links. If generator
i and j can communicate directly to each other, we denote
(i, j) ∈ E. Obviously, edges in the communication graph E
can be different from lines in the power network E. For the
communication network, we make the following assumption:
A1: The communication graph E is undirected and connected.

i∈N
g

i∈NU G

Pig ≤ P i ,

i ∈ NCG

(8c)

where Pig∗ is the mechanical power of uncontrollable generator in the steady state. In (8a), fi (Pig ) concerns the controllable
generation Pig , satisfying the following assumption:
A2: The objective fi (Pig ) is second-order continuously dif0
ferentiable, strongly convex and fi (Pig ) is Lipschitz
continuous with Lipschitz constant li > 0. i.e. ∃ αi >
00
0, αi ≤ fi (Pig ) ≤ li .
To ensure the feasibility of the optimization problem, we
make an additional assumption.
A3: The system satisfies
X
X
X
X
g
P gi ≤
pi −
Pig∗ ≤
Pi
(9)
i∈NCG

i∈N

i∈NU G

i∈NCG

Specifically, we say A3 is strictly satisfied if all the inequalities
in (9) strictly hold.
B. Equivalent Optimization Model with Virtual Load Demands
In (8b), load demands are injected to every bus, which
sometimes cannot be measured accurately if at all. As a
consequence, the values of pi may be unknown to both the
controllable generators i, i ∈ NCG and the uncontrollable
generators i, i ∈ NU G . To circumvent such an obstacle in
design, we introduce a set of new variables, p̂i , to re-formulate
SFC as the following equivalent problem:
X
fi (Pig )
(10a)
ESFC:
min
g
i∈NCG
Pi ,i∈NCG
X
X
s.t.
Pig =
p̂i
(10b)
i∈NCG

i∈NCG

g

P gi ≤ Pig ≤ P i ,

i ∈ NCG

(10c)

where p̂i is the virtual load demand supplied byP
generator i in
the
steady
state,
which
is
a
constant,
satisfying
i∈NCG p̂i =
P
P
g∗
i∈N pi −
i∈NU G Pi . Obviously, the number of virtual
loads should be equal to that of the controllable generators.
Note that the power balance constraint (8b) only requires
that all the generators supply all the loads while it is not
necessary to figure out which loads are supplied exactly by
which generators. Hence we treat virtual load demands p̂i as
the effective demands supplied by generator i for dealing with
the issue that onlyPa part of generators
are controllable.
P
P
Simply letting i∈NCG p̂i = i∈N pi − i∈NU G Pig∗ , we
immediately have the following Lemma:
Lemma 1. The problems SFC (8) and ESFC (10) have the
same optimal solutions.

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

4

IV. C ONTROLLER D ESIGN
A. Distributed Frequency Control of Controllable Generators
1) Controller design based on primal-dual gradient algorithm: Invoking the primal-dual gradient algorithm, the
Lagrangian of the ESFC (10) is given by

X
X
X
p̂i
fi (Pig ) + µ
Pig −
L=
i∈NCG

i∈NCG
g

+ γi− (P gi − Pig ) + γi+ (Pig − P i ),

i∈NCG

i ∈ NCG (11)

where µ, γi− , γi+ are Lagrangian multipliers. Based on primaldual update, the controller for i ∈ NCG is designed as


0
ugi = Pig /Ti − kPig ωi + (fi (Pig ) + µ − γi− + γi+ ) (12a)

X
X
p̂i
(12b)
Pig −
µ̇ = kµ
=

kγi [P gi

−

+
Pig ]γ −
i
i
g +

h
γ̇i+ = kγi Pig − P i

(12c)
(12d)

γi+

where kPig , kµ , kγi are positive constants; ugi the control input;
0
fi (Pig ) the marginal cost at Pig . For any xi , ai ∈ R, the
operator is defined as: [xi ]+
ai = xi if ai > 0 or xi > 0; and
[xi ]+
ai = 0 otherwise.
2) Estimating µ by second-order consensus: In (12b), µ is
a global variable, which is a function of mechanical powers
and loads of the entire system. To circumvent the obstacle
in implementation, a second-order consensus based method is
utilized to estimate µ locally by using neighboring information
only. Specifically, for i ∈ NCG , the controller is revised to:


0
ugi = Pig /Ti − kPig ωi + fi (Pig ) + µi − γi− + γi+ (13a)


X
X
g
µ̇i = kµi Pi − p̂i −
(µi − µj ) −
zij
(13b)
j∈Nci

żij = kzi (µi − µj )
γ̇i−

=

kγi [P gi
h

−

j∈Nci

(13c)

+
Pig ]γ −
i
i
g +

γ̇i+ = kγi Pig − P i

where φωi is a positive definite function, and φωi = 0
holds only when Pig = Pig∗ .
The condition C1 on the active power dynamics of uncontrollable generators is easy to verify. As an example, the
commonly-used primary frequency controller
ugi = −ωi + ωi∗ − kωi (Pig − Pig∗ ) + Pig /Ti

(14)

satisfies C1 whenever kωi > 0. In this case, we have Sωi =
g
g∗ 2
g
g∗ 2
k1
2 (Pi − Pi ) with k1 > 0 and φ = k2 (Pi − Pi ) with
0 < k2 ≤ kωi · k1 .

i∈NCG

i∈NCG

γ̇i−

input −ωi and output Pig , i.e., there exists a continuously
differentiable, positive semidefinite function Sωi such that

Ṡωi ≤ (−ωi − (−ωi∗ )) Pig − Pig∗ − φωi (Pig − Pig∗ )

γi+

C. Excitation Voltage Dynamics of All Generators
Similar to the uncontrollable generators, the following
sufficient condition on excitation voltage dynamics of all
generators is needed to guarantee system stability, since we
do not design specific excitation voltage controllers here.
C2: The excitation voltage dynamics are strictly incrementally
output passive in terms of the input −Eqi and output
Ef i , i.e., there exists continuously differentiable, positive
semidefinite function SEi such that


∗
ṠEi ≤ −Eqi − (−Eqi
) Ef i − Ef∗i − φEi (Ef i − Ef∗i )
where φEi is a positive definite function, and φEi = 0
holds only when Ef i = Ef∗i .
C2 is also easy to satisfy. As an example, it can be verified
that the controller given in [23]
∗
h(Ef i , Eqi ) = −Ef i + Ef∗i − kEi (Eqi − Eqi
)

with kEi > 0 satisfies C2. In this case, SEi = k23 (Ef i − Ef∗i )2
with k3 > 0 and φ = k4 (Ef i − Ef∗i )2 with 0 < k4 ≤ kEi · k3 .

(13d)
(13e)

where, kµi , kzi are positive constants; Nci the set of neighbors
of bus i in the communication graph; µi the local estimation of
µ. Here, (13b) and (13c) are used to estimate µ locally, where
only neighboring information is needed. zij is an auxiliary
variable to guarantee the consistency of all µi .
For the Larangian multiplier µ, −µ is often regarded as the
marginal cost of generation. Theoretically, −µi should reach
consensus for all the generators in the steady state.
P Since µ̇i =
0 holds in the steady state, we have Pig − p̂i − j∈Nci zij = 0.
Hence, zij can be regarded as the virtual line power flow of
edge (i, j) in the communication graph.

(15)

V. O PTIMALITY AND S TABILITY
After implementing the controller on the physical power
system, the closed-loop system reads

(1) − (4), (6b) − (6e), (7a), (7b)
(16)
(13a) − (13e)
In this section, we prove the optimality and stability of the
closed-loop system (16).
A. Optimality
v(t) =
 Denote the trajectory of closed-loop system as

0
g
−
+
η(t), ω(t), ω̃(t), P (t), µ(t), z(t), γ (t), γ (t), Eq (t), V (t) .
Define the equilibrium set of (16) as
V := {v ∗ |v ∗ is an equilibrium of (16)}

(17)

B. Active Power Dynamics of Uncontrollable Generators

We first present the following Theorem.

To guarantee system stability, a sufficient condition is given
for the active power dynamics of uncontrollable generators.
C1: The active power dynamics of uncontrollable generators
are strictly incrementally output passive in terms of the

Theorem 2. Suppose assumptions A1, A2 and A3 hold. In
equilibrium of (16), following assertions are true.
g
1) The mechanical powers Pig satisfy P gi ≤ Pig∗ ≤ P i ,
∀i ∈ NCG .

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

2) System frequency recovers to the nominal value, i.e.
ωi∗ = 0, ∀i ∈ NCG ∪ NU G , and ω̃i∗ = 0, ∀i ∈ N .
0
3) The marginal generation costs satisfy fi (Pig∗ ) − γi−∗ +
0
g∗
γi+∗ = fj (Pj ) − γj−∗ + γj+∗ , i, j ∈ NCG .
4) Pig∗ is the unique optimal solution of SFC problem (8).
5) µ∗i is unique if A3 is strictly satisfied.
The detailed proof of Theorem 2 is given in Appendix.A of
this paper. It shows that the nominal frequency is recovered
at equilibrium, and marginal generation costs are identical
for all controllable generators, implying the optimality of any
equilibrium.

5

Since the voltage phase deviation between two neighboring
buses is not large in practice, A4 is usually satisfied. Detailed
explanations can be found in Appendix.B of this paper.
The following stability result can be obtained.
Theorem 3. Suppose A1–A4 and C1, C2 hold. For every v ∗ ,
there exists a neighborhood S around v ∗ where all trajectories
v(t) satisfying (16) starting in S converge to the set V. In
addition, each trajectory converges to an equilibrium point.
We can further prove that ∇2 W > 0 and Ẇ ≤ 0. Moreover,
Ẇ = 0 holds only at equilibrium point. Then the theorem
can be proved using the LaSalle’s invariance principle [24,
Theorem 4.4]. Details of the proof are given in Appendix.B.

B. Stability
In this section, the stability of the closed-loop system (16)
is proved. First we define a function as
X
X
µi (Pig − p̂i )
fi (Pig∗ ) +
L̂ :=
i∈NCG
i∈NCG
 X

X
1X
−
µi zij −
µi
(µi − µj )
i∈N
i∈NCG
j∈Ni
2
X CG g
X
g
g
g
−
+
+
γi (P i − Pi ) +
(18)
γi (Pi − P i )
i∈NCG

i∈NCG

VI. I MPLEMENTATION BASED ON F REQUENCY
M EASUREMENT
A. Estimation and Optimality
Note that virtual load demands p̂i used in the controller
(13) are difficult to directly measure or estimate
P in practice.
Lemma
1
implies
that
any
p̂
are
valid
as
long
as
i
i∈NCG p̂i =
P
P
g∗
p
−
P
.
Noticing
that
the
power
i∈N i
i∈NU G i
P imbalance
is
very
small
in
normal
operation,
we
have
i∈NCG Pei ≈
P
p̂
.
In
fact,
they
are
identical
in
steady
state. Hence,
i∈N i
we specify p̂i = Pei , which implies Pig − p̂i = Pig − Pei =
Mi ω̇i + Di ωi . That leads to an estimation algorithm of µi
X
X
µ̇i = kµi −
(µi − µj ) −
zij + Mi ω̇i + Di ωi

Denote x1 := (P ), x2 := (µ, z, γi− , γi+ ), x := (x1 , x2 ). Then
L̂(x1 , x2 ) is convex in x1 and concave in x2 .
Before giving the main result, we construct a Lyapunov
candidate function composed of four parts: the quadratic part,
the potential energy part, conditions C1 and C2 related parts,
j∈Ni
j∈Ni

0
as we now explain.
g
+ τi (−µi − fi (Pi ) + γi− − γi+ )
(23)
For i ∈ NG , the quadratic part is given by
where 0 < τi < 4/li . This way, we only need to measure
X 1
1
Wk (ω, x) =
Mi (ωi − ωi∗ )2 + (x − x∗ )T K −1 (x − x∗ ) frequencies ωi at each bus locally, other than the global load
2
2
i∈NG
demands. Since the controller only needs µi of neighboring
(19) buses in the communication graph, it is easy to implement.
Now, we reconstruct the closed-loop system by replacing
where K = diag(kPig , kµi , kzi , kγi ) is a diagonal positive
(13b)
with (23) in (16), which is
definite matrix.
0

Denoting xp = (Eqi , Vi , δi , θi ), the potential energy part is
(1) − (4), (6b) − (6e), (7a), (7b)
(24)
(13a), (13c) − (13e), (23)
∗ T
∗
∗
Wp (xp ) = W̃p (xp ) − (xp − xp ) ∇xp W̃p (xp ) − W̃p (xp ) (20)
We have the following lemma.
where,
Lemma 4. Assertions 1)-5) in Theorem 2 still hold for the
X
X
0
1
W̃p (Eqi , Ei , Vi , δi , θi ) =
Bii Vi2 +
pi θi
equilibrium of (24).
i∈N 2
i∈N
X
1X X
The proof of Lemma 4 is given in Appendix.C.
−
qi ln Vi −
Vi Vj Bij cos (θi − θj )
(21)
2
g

i∈N

i∈N j∈Ni

0
 0 2
X Eqi Vi
X
xdi
cos(δi − θi ) +
−
Eqi
0
0
0
xdi
2xdi (xdi − xdi )
i∈NG
i∈NG
P
PConditions 1C1 and C2 related parts are i∈NU G Sωi and
i∈NG T 0 (xdi −x0 ) SEi respectively.
d0i
di
The Lyapunov function is defined as
X
X
SEi
W =Wk + Wp +
Sωi +
(22)
0
0
Td0i (xdi − xdi )
i∈N
i∈N
UG

G

Then, we give the following assumption.
A4: The Hessian of Wp satisfies ∇2v Wp (v) > 0 at desired
equilbirium.

B. Discussion on Stability
Recall (2b), then (23) is derived to
X
X
µ̇i = kµi Pig − p̂i + p̂i − Pei −
(µi − µj ) −
zij
j∈Ni

+ τi (−µi −

0

fi (Pig )

+

γi−

j∈Ni

−

γi+ )



(25)

∗
Denote ρi = p̂i − Pei = Pei
− Pei , which is the difference of
electric power and its value in the steady state. We have the
following assumption
A5: The disturbance can be written as ρi = βi (t)ωi , where
|βi (t)| ≤ β̄i and β̄i is a positive constant. In addition, the
set { t < ∞ | ωi (t) = ωi∗ } has a measure zero.

DFIG
Bus12
Tw XX, XX XXXX
JOURNAL OF LATEX CLASS FILES, VOL. XX, NO.
Bus1

Bus5

Bus7

Bus8

Bus9

6

Bus10 Bus11 Bus3
T3

T1
Gen1

Gen3

Bus6
exists such βi (t).
A5 argues
Whenever ωi 6= ωi∗ , thereT2 alway
T4
Load1
Load2
Bus2
Bus4
∗
that ωi (t) = ωi only happens
at isolated Gen4
points
except
Gen2
equilibrium. Generally, this is reasonable in power system.
Denote the state variables of (24) and its equilibrium set are
ṽ and Ṽ respectively. We have following stability result
Bus1 Bus5 Bus6 Bus7
Bus10 Bus11 Bus2
Theorem 5. Suppose
A1–A5, C1,Bus8C2 Bus9hold
and
(9) is not
T1
T2
∗
binding. For every Gen1
ṽ , there exists a neighborhood Gen2
S around
ṽ ∗ where all trajectories ṽ(t) satisfying (24) starting in S
Bus15 Bus16 Bus17 Bus4
Bus14 Bus13 Bus12
converge to the set ṼBus3whenever
T4
T3
p
Gen4
Gen3
(26)
β̄i < τi Di (4 − τi li ).

G37

Controller

u gj

26

l

g

27

2
18

1

We consider the following scenario: 1) at t = 10s, there
is a step change of 60MW load demands at each of buses 3,
15, 23, 24, 25; 2) at t = 70s, there is another step change of
120MW load at bus 23. Neither the original load demands nor
their changes are known to the generators.

G35

16
15

G39

21
39

4

22

14

5
6

12

7
8

19

13

11

20
23
33

10
9

31
32

G31

34
G32

36

G33
G34

G36

Fig. 2. The New England 39-bus system
Eqi' , Vi ,  i , i , i , Pi g

Eqj' , V j ,  j ,  j ,  j , Pjg

i ,  j

Controller
g
i

Controller

u gj

E fi

u

A. Test System

B. Results under Small Disturbances

G38

24

35

uil

To test the proposed controller, the New England 39-bus
system with 10 generators as shown in Fig.2, is utilized. In the
simulation, we apply (23) to estimate the virtual load demands.
All simulations are implemented in the commercial power
system simulation software PSCAD [25], and are carried on
a notebook with 8GB memory and 2.39 GHz CPU.
We control only a subset of these generators, namely G32,
G36, G38, G39, while the remaining are equipped with the
primary frequency control given in (14). In particular, we apply
the controller (13) derived based on a simple model to a much
more realistic and complicated model in PSCAD. The detailed
electromagnetic transient model of three-phase synchronous
machines (sixth-order model) is adopted to simulate generators
with governors and exciters. All the lines and transformers
take both resistance and reactance into account. The loads are
modeled as fixed loads in PSCAD. The communication graph
is undirected and set as G32 ↔ G36 ↔ G38 ↔ G39 ↔ G32.
The objective function is set as fi = 21 ai (Pig )2 + bi Pig ,
which is the generation cost of generator i. Capacity limits of
Pig and parameters ai , bi are given in Table I.
The closed-loop system is shown in Fig.3, where each
generator only needs to measure local frequency, mechanical
power, voltage and phase angle to compute its control command. Note that there is no load measurement and only µi are
communicated between neighboring controllable generators.

38

17

3

Controller

VII.u j S IMULATION Rui ESULTS

29

28

30

Moreover, the convergence of each such trajectory is to a point.
The proof of Theorem 5 is given in Appendix.D. In
fact, the range of β(t) can be very large as long as li is
small enough. For example, set τi = 3/li , then β̄i (t) <
p
3DP
i /li . As we know, if we change the objective function
to k i∈NCG fi (Pig ), k > 0, the optimal solution will not
change. Thus, li can be very small as long as k is small
enough. In this regard, (26) is not conservative. Moreover,
g
l
g
i , Pieven
, Pi l , Pij(9)
, Pki is binding, our
Pjk Section VII that
j , Pj , Pj , Pij ,in
we will illustrate
controller still works.

37
25

G30

E fj

…
Fig. 3. Diagram of the closed-loop system
TABLE I
C APACITY LIMITS OF GENERATORS
Gi
32
36
38
39

g

[P gi , P i ] (MW)
[0, 1000]
[0, 1000]
[0, 850]
[0, 1080]

ai
0.00009
0.00014
0.00010
0.00008

bi
0.032
0.030
0.032
0.032

TABLE II
E QUILIBRIUM POINTS

Stage 1
Stage 2

g
P32
(MW)
927
968

g
P36
(MW)
610
652

g
P38
(MW)
834
850

g
P39
(MW)
1043
1080

1) Equilibrium: In steady states, the nominal frequency is
well recovered. The optimal mechanical powers are given in
Table II, which are identical to the optimal solution of (8)
computed by centralized optimization. Stage 1 is for the period
from 10s to 70s, and Stage 2 from 70s to 130s. The values
in Table II are generations at the end of each stage. In Stage
1, no generation reaches its limit, while in Stage 2 both G38
and G39 reach their upper limits. At the end of each stage, the
marginal generation cost −µi of generator i, converges identically (see Fig. 9), implying the optimality of the results. The
test results confirm the theoretical analyses and demonstrate
that our controller can automatically attain optimal operation
points even in the more realistic and sophisicated model.
2) Dynamic Performance: In this subsection, we analyze
the dynamic performance of the closed-loop system. For

G

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

Frequency (Hz)

60.1
60
59.9

59.7
59.6

10

30

50

70

Time (s)

90

G35
G36
G37
G38
G39

G30
G31
G32
G33
G34

G35
G36
G37
G38
G39

G30
G31
G32
G33
G34

59.8

110

10

130

30

50

70

Time (s)

800

1600

With AGC

Mechanical Power (MW)

Distributed control

60.2

7

110

G32
G36
G38
G39

1200
1000

600
500

G30
G31
G33
G34
G35
G37

400

800

300

600

200
400

90

700

1400

10

30

130

50

70

Time (s)

90

110

130

10

30

50

70

Time (s)

90

110

130

Fig. 6. Dynamics of mechanical powers under AGC
Fig. 4. Dynamics of frequencies
Mechanical Power (MW)

Frequency (Hz)

60.2
60.1
60
59.9

G30
G31
G32
G33
G34

59.8
59.7
59.6

G35
G36
G37
G38
G39

30

50

70

Time (s)

90

110

130

700

1000

600

900
G32
G36
G38
G39

800
700
600

500

G30
G31
G33
G34
G35
G37

400
300

500
400

10

800

1200
1100

200
10

30

50

70

Time (s)

90

110

130

10

30

50

70

Time (s)

90

110

130

Fig. 7. Dynamics of mechanical powers under the proposed control

Fig. 5. Dynamics of frequencies with the controller in [11]

Voltage (p.u.)

Generation (MW)

perfectly in coincidence with the indication given by Theorems
comparison, automatic generation control (AGC) is tested in
2, 3 and 5 and Lemma 4 in this paper.
the same scenario. In the AGC implementation, the signal of
Mechanical power dynamics under the AGC and the proarea control error (ACE) is given by ACE = Kf ω + Pij
posed controller are shown in Fig.6 and Fig.7, respectively.
[26, Chapter 11.6], where, Kf is the frequency response
The left parts show mechanical powers of G32, G36, G38,
coefficient; ω the frequency deviation; Pij the deviation of
G39, while the right parts show mechanical powers of other
tie line power. In the case studies, we can treat the whole
generators adopting conventional controller (14). With both
system as one control area, implying Pij = 0. Hence, the
controls, mechanical powers of the generators adopting (14)
control center computes ACE = Kf ω and allocates it to
remain identical in the steady state. However, there are two
AGC generators, G32, G36, G38 and G39. In this situation,
problems when adopting the AGC: 1) mechanical powers
the
P control command of each generator is ri · ACE, where are not optimal; 2) mechanical power of G39 violates the
i ri = 1. In this paper, we set ri = 0.25 for i = 1, 2, 3, 4.
capacity limit. In contrast, the proposed control can avoid these
The trajectories of frequencies are given in Fig.4, where
problems. In Stage 1 of Fig.7, no generator reaches capacity
the left one stands for the proposed controller and the right
limits. In Stage 2, both G38 and G39 reach their upper limits.
one for the AGC. It is shown in Fig.4 that the frequencies
Then, G38 and G39 stop increasing their outputs while G32
are recovered to the nominal value under both controls. The1200
800
and G36 raise their outputs to balance
the load demands. In
700
frequency drops under two controls are very similar while the1100
addition,
the
steady
states
of
in
both
stages
are optimal, which
recovery time under the proposed control is much less than1000
600
900 are the same as shown in Table II.
that under the conventional AGC.
500
G30
800
We also illustrate in G32
Fig.8 the dynamics of voltage at
G36
G31
G38
Although a number of studies have been devoted to dis- 700
G33
400
G34
buses 3, 15, 23, 24, 25.G39The voltages
converge rapidly, and
G35
tributed frequency control of power systems, most of them 600
300
G37
only
experience
small
drops
when
loads
increase.
This
result
500
assume that all the nodes are controllable except [11]. To
200
validates the effectiveness of the voltage control.
10
30
50
70
90
110
make a fair comparison, the controller proposed in [11] is 400 10 30 50 70 90 110 130
Time generation
(s)
Time
The marginal
cost of generator i, −µ
are shown
i ,(s)
adopted as a rival in our tests. As shown in Eq.(8) of [11], each
controller needs to predict the load it should supply in steady
1.1
state. However, it is hard to acquire an accurate prediction in
1.08
practice, which could lead to steady-state frequency error.
1.06
1.04
In the next case, we compare the two controls in the same
1.02
scenario as that in Section VII.B. The dynamics of frequencies
1
0.98
Bus3
with the controller given in [11] are shown in Fig.5. It is
Bus15
0.96
Bus23
observed that there is a frequency deviation in steady state, as
Bus24
0.94
Bus25
0.92
the prediction is inaccurate. Although the deviation is usually
0.9
10
30
50
70
90
110
130
quite small, it is difficult to completely eliminate. In contrast,
Time (s)
when the proposed method is adopted, there is no frequency
deviation in steady state, as is shown in Fig.4. This result is Fig. 8. Dynamics of bus voltages

130

1.05
1

10
50
70
90
130
JOURNAL
OF LAT30EX CLASS
FILES,
VOL.
XX, 110
NO. XX,
XX XXXX
Time (s)

G32
G36
G38
G39

1.05
10

30

50

70

Time (s)

90

110

130

30

50

70

Time (s)

90

110

59.5
59.4
0 10

20 30 40 50 60 70 80 90
Time (s)

1000

G32
G36
G38
G39

500

0

0 10

20 30 40 50 60 70 80 90
Time (s)

Mechanical Power (MW)

59.6

Mechanical Power (MW)

Frequency (Hz)

59.7

59.3

59.7

130

1500

G30
G31
G33
G34
G35
G36
G37
G38
G39

60
59.9

1
0.8
0.6
0.4
0.2

0

20

40

60

80

Time (s)

0

100

0

20

40

60

80

Time (s)

100

Fig. 11. Dynamics of frequencies and voltages with line trip

60.1
60
59.8

60.1

59.8

10

Fig. 9. Dynamics of −µ and z

59.9

60.2

Bus4
Bus14

1.2

G35
G36
G37
G38
G39

1200

955

1100

950

1000

945
30

50

70

90

110

G32
G36
G38
G39

G32
G36
G38
G39

1.15
1.1

900
800



1.1

G30
G31
G32
G33
G34

Voltage (p.u.)

60.3

Frequency (Hz)

1.15

0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

z39-32
z39-38
z32-36
z36-38

z



1.2

1

60.4

103

1.3
1.25

8

1.095
1.09

1.05

1.085

700

1

1.08

600

1.075

0.95

500
400

0

20

40

60

Time (s)

80

100

0.9

1.07
1.065
30

0

20

40

50

70

60

Time (s)

90

80

110

100

Fig. 10. Dynamics of frequencies and mechanical powers with generator trip

Fig. 12. Dynamics of −µ and mechanical powers with line trip

in the left part of Fig.9. They converge in both stages and the
steady-state values in Stage 2 are slightly bigger than that in
Stage 1, as the load changes lead to an increase in the marginal
generation cost. Dynamics of zij , (i, j) ∈ E are illustrated
in the right part of Fig.9, which demonstrate that the steady
state values do not change in the two stages. In addition, the
variation of z in transient is very small as the deviation of µi
is very small.

line is re-closed. At the same time, the voltages of buses 4 and
14 drop to nearly zero when the fault happens. The voltages are
then stabilized to a new steady-state value in around 10s after
the fault line is tripped. They are slightly different from their
initial values because the systems operating point has changed
due to line tripping. When the tripped line is re-closed, the
voltages recover to the initial values quickly.
When line (4, 14) is tripped, the power flow across the
power network varies accordingly. As a consequence, the line
loss also changes, causing variations of mechanical powers,
as shown in Fig.12. In the left part of Fig.12, the inset is
the dynamics of mechanical power of G39 from 30s to 110s.
Similarly, the inset in the right part is dynamics of −µ of
all generators. Mechanical powers and their marginal costs
all increase when the line is tripped. However, the proposed
controller compensates the loss change autonomously.
These simulation results demonstrate that the proposed
distributed optimal frequency controller can cope with large
disturbances such as generator tripping and short-circuit fault.

C. Performance under Large Disturbances
In this subsection, two scenarios of large disturbances are
considered. One is a generator tripping and the other is a shortcircuit fault followed by a line tripping.
1) Generator tripping: At t = 10s, G32 is tripped, followed by occurrence of certain power imbalance. Note that the
communication graph still remains connected. The output of
G32 drops to zero. Frequency and mechanical powers change
accordingly. System dynamics are illustrated in Fig.10. The
left part of Fig.10 shows the frequency dynamics, and the right
shows the mechanical power dynamics of controllable generators. It is observed that the system frequency experiences a big
drop at first, and then recovers to the nominal value quickly as
other controllable generators increase outputs to balance the
power mismatch. These results confirm that our controller can
adapt to generator tripping autonomously even if the tripped
generator is contributing to frequency control.
2) Short-circuit fault: At t = 10s, there occurs a threephase short-circuit fault on line (4, 14). At t = 10.05s, this line
is tripped by breakers. At t = 60s, the fault is removed and line
(4, 14) is re-closed. Frequency dynamics and voltage dynamics
of buses 4 and 14 are given in Fig.11, where the left part shows
frequency dynamics and the right shows voltage dynamics. It
can be seen that the frequency experiences violent oscillations
after the fault happens. And then it is stabilized quickly once
the line is tripped. Small frequency oscillation occurs when the

VIII. C ONCLUSION
In this paper, we have designed a distributed optimal frequency control using a nonlinear network-preserving model,
where only a subset of generator buses is controlled. We have
also simplified the implementation by relaxing the requirements of load measurements and communication topology.
Since nonlinearity of power flow model and dynamics of excitation voltage has been taken into account, our controllers can
cope with large disturbances. We have proved that the closedloop system asymptotically converges to the optimal solution
of the economic dispatch problem. We have also carried out
substantial simulations to verify the good performance of our
controller under both small and large disturbances.
In this work, we have not considered the constraints on
line flows, since the controllable generators are selected arbi-

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

trarily and may not suffice for congestion management. An
interesting problem is how to find out the minimal set of
controllable generators to fulfill the requirement of congestion
management, which is our future work.
R EFERENCES
[1] D. Cai, E. Mallada, and A. Wierman, “Distributed optimization decomposition for joint economic dispatch and frequency regulation,” IEEE
Trans. Power Syst., early access, 2017.
[2] D. J. Shiltz, S. Baros, M. Cvetkovi?, and A. M. Annaswamy, “Integration
of automatic generation control and demand response via a dynamic
regulation market mechanism,” IEEE Trans. Control Syst. Technology,
vol. PP, no. 99, pp. 1–16, 2017.
[3] A. Jokić, M. Lazar, and P. P. van den Bosch, “Real-time control of power
systems using nodal prices,” Int. J. Elect. Power Energy Syst., vol. 31,
no. 9, pp. 522–530, 2009.
[4] F. Dorfler, J. W. Simpson-Porco, and F. Bullo, “Breaking the hierarchy:
Distributed control and economic optimality in microgrids,” IEEE Trans.
Control Netw. Syst., vol. 3, no. 3, pp. 241–253, Sept 2016.
[5] N. Li, C. Zhao, and L. Chen, “Connecting automatic generation control
and economic dispatch from an optimization view,” IEEE Trans. Control
Netw. Syst., vol. 3, no. 3, pp. 254–264, Sept 2016.
[6] D. Feijer and F. Paganini., “Stability of primal-dual gradient dynamics
and applications to network optimization,” Automatica, vol. 46, no. 12,
pp. 1974–1981, Dec 2010.
[7] C. Zhao, U. Topcu, N. Li, and S. H.Low., “Design and stability of loadside primary frequency control in power systems,” IEEE Trans. Autom.
Control, vol. 59, no. 5, pp. 1177–1189, Jan. 2014.
[8] E. Mallada, C. Zhao, and S. H. Low, “Optimal load-side control for
frequency regulation in smart grids,” IEEE Trans. Autom. Control, early
access, 2017.
[9] A. Kasis, E. Devane, C. Spanias, and I. Lestas, “Primary frequency
regulation with load-side participation part i: stability and optimality,”
IEEE Trans. Power Syst. early access, 2016.
[10] C. Zhao, E. Mallada, S. Low, and J. Bialek, “A unified framework
for frequency control and congestion management,” in Power Systems
Computation Conference (PSCC), 2016. IEEE, 2016, pp. 1–7.
[11] J. Z. Pang, L. Guo, and S. H. Low, “Load-side frequency regulation with
limited control coverage,” ACM SIGMETRICS Performance Evaluation
Review, vol. 45, no. 2, pp. 94–96, 2017.
[12] Z. Wang, F. Liu, S. H. Low, C. Zhao, and S. Mei, “Decentralized optimal frequency control of interconnected power systems with transient
constraints,” in Proc. 55th IEEE Conf. Decis. Control (CDC), Dec 2016,
pp. 664–671.
[13] X. Zhang, N. Li, and A. Papachristodoulou, “Achieving real-time economic dispatch in power networks via a saddle point design approach,”
in 2015 IEEE PES General Meeting, July 2015, pp. 1–5.
[14] Z. Wang, F. Liu, S. H. Low, C. Zhao, and S. Mei, “Distributed frequency
control with operational constraints, part i: Per-node power balance,”
IEEE Trans. Smart Grid, early access, 2017.
[15] ——, “Distributed frequency control with operational constraints, part ii:
Network power balance,” IEEE Trans. Smart Grid, early access, 2017.
[16] T. Stegink, C. D. Persis, and A. van der Schaft, “A port-hamiltonian
approach to optimal frequency regulation in power grids,” in 54th IEEE
Conference on Decision and Control (CDC), Dec 2015, pp. 3224–3229.
[17] T. W. Stegink, C. De Persis, and A. J. van der Schaft, “Stabilization
of structure-preserving power networks with market dynamics,” arXiv
preprint arXiv:1611.04755, 2016.
[18] T. Stegink, C. D. Persis, and A. van der Schaft, “A unifying energy-based
approach to stability of power grids with market dynamics,” IEEE Trans.
Autom. Control, vol. 62, no. 6, pp. 2612–2622, June 2017.
[19] X. Zhang and A. Papachristodoulou., “A real-time control framework for
smart power networks: Design methodology and stability,” Automatica,
vol. 58, pp. 43–50, Aug 2015.
[20] C. Zhao and S. Low, “Optimal decentralized primary frequency control
in power networks,” in Proc. 53rd IEEE Conf. Decis. Control (CDC),
Dec 2014, pp. 2467–2473.
[21] S. Trip and C. De Persis, “Optimal generation in structure-preserving
power networks with second-order turbine-governor dynamics,” in Control Conference (ECC), 2016 European. IEEE, 2016, pp. 916–921.
[22] P. Kundur, Power System Stability and Control. McGraw-hill New
York, 1994, vol. 7.
[23] Y. Liu, C. Li, and Y. Wang, “Decentralized excitation control of multimachine multi-load power systems using hamiltonian function method,”
Acta Automatica Sinica, vol. 35, no. 7, pp. 919–925, 2009.

9

[24] H. K. Khalil, Nonlinear systems. Prentice hall New Jersey, 1996, vol. 3.
[25] https://hvdc.ca/pscad/.
[26] V. V. Arthur and R. Bergen, Power System Analysis. Upper Saddle
River: NJ: Pretice Hall, 2000.
[27] S. Boyd and L. Vandenberghe, Convex optimization.
Cambridge
university press, 2004.
[28] C. D. Persis, N. Monshizadeh, J. Schiffer, and F. Drfler, “A lyapunov
approach to control of microgrids with a network-preserved differentialalgebraic model,” in 2016 IEEE 55th Conference on Decision and
Control (CDC), Dec 2016, pp. 2595–2600.
[29] P. Lin, Y. Jia, and L. Li, “Distributed robust h∞ consensus control in
directed networks of agents with time-delay,” Syst. Control Lett., vol. 57,
no. 8, pp. 643–653, 2008.

A PPENDIX
A. Proof of Theorem 2
Proof. From γ̇i− = γ̇i+ = 0 in (13d), it follows that P gi ≤
g
Pig∗ ≤ P i , which is the first assertion. From żij = 0 in (13c),
we get µ∗i P
= µ∗j = µ0 . Set
P µ̇i = 0, add
P (13b) for i ∈ NCG ,
and recall i∈NCG p̂i = i∈N pi − i∈NU G Pig∗ , we have
X
X
X
Pig∗ −
pi +
Pig∗ = 0 (A.1)
i∈NCG

i∈N

i∈NU G

The right sides of (7a) and (7b) vanish in the equilibrium
points, which implies ωi∗ = ω̃i∗ = ω̃j∗ = ω0 . Set ω̇i = 0 and
add (2b), (6b) and (6c). We have
X
X
X
X
ω0
D̃i =
Pig∗ −
pi +
Pig∗ = 0 (A.2)
i∈N

i∈NCG

i∈N

i∈NU G

This implies that ω0 = 0 due to D̃i > 0, which is the second
assertion.
Combine (13a), (2d) and Ṗig = 0. We have
0

fi (Pig∗ ) − γi−∗ + γi+∗ + ωi∗ + µ∗i = 0
Since ω0 = 0, µ∗i = µ∗j = µ0 , we can obtain the third
assertion.
Next we will prove assertion 4). Since all the constraints
of SFC are linear, A3 implies that Slater’s condition holds
[27, Chapter 5.2.3]. Moreover, the objective function is strictly
convex. We only need to prove that (Pig∗ , µ0 , γi−∗ , γi+∗ ) satisfies the KKT condition of SFC in order to prove the fourth
assertion.
The KKT conditions of SFC problem (8) are
0

fi (Pig∗ ) − γi−∗ + γi+∗ + µ0 = 0
X
X
X
Pig∗ −
pi +

i∈NCG
g
P gi ≤ Pig∗ ≤ P i
γi−∗ ≥ 0, γi+∗ ≥ 0
γi−∗ (P gi − Pig∗ ) =

i∈N

(A.3a)

i∈NU G

Pig∗

=0

(A.3b)
(A.3c)
(A.3d)

0, γi+∗ (Pig∗

−

g
Pi )

=0

(A.3e)

From γ̇i− = γ̇i+ = 0, we have (A.3c), (A.3d) and (A.3e). From
the third assertion, we have (A.3a). From ω̇ = 0 and the second
assertion, we have (A.3b). Therefore, the equilibrium points of
the closed-loop system (16) satisfy the KKT conditions (A.3).
This implies the fourth assertion.
If A3 is strictly satisfied, we know ∃ i ∈ NCG that γi−∗ =
0
+∗
γi = 0. Then, µ∗i = −fi (Pig∗ ) is uniquely determined by
Pig∗ , implying the last assertion. This completes the proof.

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

10

Proof of Theorem 3. Recall (18), and dynamics of (13a) (13e) are rewritten as


(B.1a)
Ṗig = −kPig · ωi + ∂ L̂(x1 , x2 )/∂Pig
µ̇i = kµi · ∂ L̂(x1 , x2 )/∂µi

(B.1c)

g

σ + := {i ∈ NCG | γi+ = 0, Pig − P i < 0}
σ := {i ∈

= 0,

−

Pig

< 0}

(ωi −

X

− Di ωi − Pei ) −

−

(Pig

Pig∗ )(ωi

−

− ωi∗ )

Since L̂ is a convex function of x1 and a concave function of
x2 , we have
≤L̂(x∗1 , x2 ) − L̂(x1 , x2 ) + L̂(x1 , x2 ) − L̂(x1 , x∗2 )
≤0

(B.3b)

(B.7)

where the first inequality follows because L̂ is convex in x1
and concave in x2 and the second inequality follows because
(x∗1 , x∗2 ) is a saddle point. Hence, we have
X
X
Ẇk ≤ −
Di (ωi − ωi∗ )2 −
D̃i (ω̃i − ω̃i∗ )2
i∈NG

(Pig

−

Pig∗ )ωi

X

i∈N

(ωi − ω̃i )(Pei −

∗
Pei
)

i∈NG

+

X

−

(ω̃i − ω̃j )(Pij − Pij∗ )

(i,j)∈E

X
i∈NU G

(Pig

−

Pig∗ )(ωi

− ωi∗ )

(B.8)

The partial of Wp (xp ) is
1
∗
0 (Eqi − Eqi )
xdi − xdi
∇V i W p = 0

∇ E 0 Wp =

· (ωi − ωi∗ )

(B.9a)

qi

∇δi Wp = Pei −
X
X
∗
∇θ i W p =
(Pij − Pij∗ ) −
(Pei − Pei
)
∗
Pei
)

(B.9c)
(B.9d)

i∈NG

(i,j)∈E

The derivative of Wp is

(Pig − Pig∗ )(ωi − ωi∗ )

(B.4)
Ẇp =

∗
X (Eqi − Eqi
)(Ef i − Ef∗i )
0

−
−∗ T
) (P g − P g )
(γ − − γ −∗ )T [P g − P g ]+
γ − ≤ (γ − γ

X

+

0

Td0i (xdi − xdi )

i∈NG

where the inequality is due to

(ωi − ω̃i )(Pei −

∗
Pei
)

+

i∈NG

X

∗ 2
(Eqi − Eqi
)
0
0
Td0i (xdi − xdi )

(ω̃i − ω̃j )(Pij − Pij∗ )

(i,j)∈E

(B.10)

P gi −Pig

Here the inequality holds since γi− = 0 ≤ γi−∗ and
<
0 for i ∈ σ − , i.e., (γi− − γi−∗ ) · (P gi − Pig ) ≥ 0. Similarly for
i ∈ σ + , we have
g

X

−

i∈NG

= (γ − − γ −∗ )T ∇γ − L̂

g

+
+∗ T
(γ + − γ +∗ )T [P g − P ]+
) (P g − P )
γ + ≤ (γ − γ

= (γ + − γ +∗ )T ∇γ + L̂
From (6b) and (6c), it yields


X
X
∗
∗
∗
0=
(ω̃i − ω̃i ) (Pei − Pei ) −
(Pij − Pij ) (B.5)
j∈Ni

(B.9b)

∗
Pei

i∈NU G

i∈NG

(B.6)

(B.3a)

− (x1 − x∗1 )T · ∇x1 L̂ + (x2 − x∗2 )T · ∇x2 L̂
= −(x1 − x∗1 )T · ∇x1 L̂ + (x2 − x∗2 )T · ∇x2 L̂
X
X
−
Di (ωi − ωi∗ )2 −
(ωi − ωi∗ )(Pei −
i∈NG
i∈NG
X

(i,j)∈E

X

(B.2b)

i∈NCG

+

i∈NG

(ω̃i − ω̃j )(Pij − Pij∗ )

−(x1 − x∗1 )T · ∇x1 L̂(x1 , x2 ) + (x2 − x∗2 )T · ∇x2 L̂(x1 , x2 )

i∈NG

−

X

∗
(ω̃i − ωi )(Pei − Pei
)−

(B.2a)

− (x1 − x∗1 )T · ∇x1 L̂ + (x2 − x∗2 )T · ∇x2 L̂
X
∗
=
(ωi − ωi∗ )(Pig − Pig∗ − Di (ωi − ωi∗ ) − (Pei − Pei
))
Pig∗ )

X

− (x1 − x∗1 )T · ∇x1 L̂ + (x2 − x∗2 )T · ∇x2 L̂

i∈NCG

(Pig

i∈N

i∈NG

+

−

i∈NG

X

Add (B.5) to (B.4), then we have
X
X
Di (ωi − ωi∗ )2 −
D̃i (ω̃i − ω̃i∗ )2
Ẇk ≤ −

(B.1e)

i∈NG

≤

j∈Ni

+

The derivative of Wk is
X
Ẇk =
Mi (ωi − ωi∗ )ω̇i + (x − x∗ )T · K −1 ẋ
ωi∗ )(Pig

(Pij − Pij∗ )

=L̂(x∗1 , x2 ) − L̂(x∗1 , x∗2 ) + L̂(x∗1 , x∗2 ) − L̂(x1 , x∗2 )

Then (13d) and (13e) are equivalent to

g
kγi (Pig − P i ), if i ∈
/ σ+ ;
+
γ̇i =
0,
if i ∈ σ + .

kγi (P gi − Pig ), if i ∈
/ σ− ;
−
γ̇i =
0,
if i ∈ σ − .

X

i∈NL

X

(B.1d)

With regard to the closed-loop system, we first define two
sets, σ + and σ − , as follows [6].
P gi

(ω̃i − ω̃i∗ )

i∈NU G

γi

NCG | γi−

X

D̃i (ω̃i − ω̃i∗ )2 −

i∈N

(B.1b)

żij = kzi · ∂ L̂(x1 , x2 )/∂zij
h
i+
γ̇i− = kγi · ∂ L̂(x1 , x2 )/∂γi− −
γi
h
i+
+
+
γ̇i = kγi · ∂ L̂(x1 , x2 )/∂γi +

−

X

−

B. Proof of Theorem 3

The Lyapunov
is defined as W = Wk + Wp +
P
P function
1
Sωi +
S , and its derivative is
0
0
T (x −x ) Ei

i∈NU G

i∈NG

d0i

di

di

X
X
ṠEi
Ṡωi +
Ẇ = Ẇk + Ẇp +
0
0
i∈NU G
i∈NG T
d0i (xdi − xdi )
X
X
≤−
Di (ωi − ωi∗ )2 −
D̃i (ω̃i − ω̃i∗ )2
i∈NG

−

X
i∈NG

∗ 2
(Eqi − Eqi
)
0
0
Td0i (xdi − xdi )

i∈N

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

+

X

(Pig − Pig∗ )(ωi − ωi∗ ) +

Ṡωi

i∈NU G

i∈NU G

+

X

11

∗
X (Eqi − Eqi
)(Ef i − Ef∗i )
0

0

+

X

ṠEi
0
Td0i (xdi − xdi )
0

Td0i (xdi − xdi )
i∈NG
i∈NG
X
X
∗ 2
Di (ωi − ωi ) −
D̃i (ω̃i − ω̃i∗ )2
≤−
i∈NG

−

X
i∈NG

≤0

i∈N

∗ 2
X
X
(Eqi − Eqi
)
φ Ei
φωi −
−
0
0
Td0i (xdi − xdi ) i∈N
i∈N
UG

G

(B.11)

The last two inequalities are due to assumption A4.
To prove the locally asymptotic stability of the closed-loop
system, we also need to prove that W > 0, ∀v ∈ S\v ∗ .
Equivalently, we need ∇2v W > 0, ∀ v ∈ S\v ∗ , i.e. A4.
Consequently, there exists a neighborhood set { v : W (v) ≤
 } for all sufficiently small  > 0 so that ∇2v W (v) > 0.
Hence, there is a compact set S around v ∗ contained in such
neighborhood, which is forward invariant. Let Z1 := { v :
Ẇ (v) = 0 }. By LaSalle’s invariance principle, the each
of trajectories v(t) starting from S converges to the largest
invariant set Z + contained in S ∩ Z1 . From above analysis, if
Ẇ (v) = 0, v is an equilibrium point of the closed-loop system
(16). Hence, v converges to Z + ∈ V.
Finally, we will prove the convergence of each v(t) starting
from V is to a point by following the proof of Theorem 1 in
[18]. Since v(t) is bounded, its ω-limit set Ω(v) 6= ∅. By contradiction, suppose there exist two different points in Ω(v), i.e.,
v1∗ , v2∗ ∈ Ω(v), v1∗ 6= v2∗ . Since the Hessian of W1 (v), W2 (v)
is positive definite in S, there exist W1 (v), W2 (v) defined by
(22) with respect to v1∗ , v2∗ and scalars c1 > 0, c2 > 0 such
that two sets W1−1 (≤ c1 ) := {v|W1 (v) ≤ c1 }, W2−1 (≤ c2 ) :=
{v|W2 (v) ≤ c2 } are disjoint (i.e. W1−1 (≤ c1 )∩W2−1 (≤ c2 ) =
∅) and compact. In addition, W1−1 (≤ c1 ), W2−1 (≤ c2 ) are
forward invariant. By (B.11), there exists a finite time t1 > 0
such that v(t) ∈ W1−1 (≤ c1 ) for ∀t ≥ t1 . Similarly, there
exists a finite time t2 > 0 such that v(t) ∈ W2−1 (≤ c2 ) for
∀t ≥ t2 . This implies that v(t) ∈ W1−1 (≤ c1 ) ∩ W2−1 (≤ c2 )
for ∀t ≥ max(t1 , t2 ), which is a contradiction. This completes
the proof.
In this part, we discuss the reasonableness of Assumption
A4 by referring to [28, Lemma 3]. Reference [28] investigates
the control of inverter-based microgrids based on a networkpreserving model, while we extend some results to more
complicated synchronous-generator based bulk power systems.
For simplicity, we first present some notations following [28].
Comparing (1) and (4), Pei , Qei have same structures with
Pij , Qij , respectively. We can treat the reactance of generator
0
as a line with admittance as 1/xdi , which connects i ∈ NG
and inner node of the generator. We denote the inner nodes
0
of generators as NG . Then, we can get a augmented power
network, the incidence matrix of which is denoted as Ĉ. The
set of nodes in the augmented power network is denoted as
0
0
0
N = N ∪ NG . Denote V̂ = (Eq , V ), θ̂ = (δ, θ). Let |Ĉ| denote the matrix obtained from Ĉ by replacing all the elements

0

cij with |cij |. Define Γ(V̂ ) :=
Define A as

 −|Bij | cos(θi − θj ),
diag(|Bii |),
Aij =
0

−cos(δi − θj )/xdi ,

diag(|Bij |Vi Vj ), i, j ∈ N .
i 6= j, i, j ∈ N ;
i = j, i, j ∈ N ;
0
i ∈ NG , j ∈ NG .

(B.12)

For simplicity, we use the following notation. For an ndimensional vector r := {r1 , r2 , · · · , rn }, the diagonal matrix diag(r1 , r2 , · · · , rn ) is denoted in short by [r]D . And
cos(·), sin(·) are defined component-wise.
From the definition of W in (22), ∇2v W (v) > 0 if and only
if ∇2v Wp (v) > 0, i.e. the matrix


Γ(V̂ )[cos(Ĉ T θ̂)]D
[sin(Ĉ T θ̂)]D Γ(V̂ )|Ĉ|T [V̂ ]−1
D
T
[V̂ ]−1
A + H(V̂ )
D |Ĉ|Γ(V̂ )[sin(Ĉ θ̂)]D
(B.13)
is positive definite, where
" h
H(V̂ ) =

xdi
0
0
2xdi (xdi −xdi )

i
D

0

0
[qi /Vi2 ]D

#
(B.14)

In any steady state of power system (i.e.), the phase-angle
difference between two neighboring nodes is usually small.
In addition, the difference between δi and θi is also small.
This implies that the matrix in (B.13) is diagonal dominant as
well as its positive definiteness. Therefore, Assumption A4 is
usually satisfied and makes sense.
C. Proof of Lemma 4
Proof. From żij = 0, we get µ∗i = µ∗j = µ0 . Set µ̇i = 0 and
add (23) for i ∈ NCG , we have
X
X
Mi ω̇i +
(Di + τi )ωi = 0
(C.1)
i∈NCG

i∈NCG

The right sides of (7a) and (7b) vanish in the equilibrium
points, which implies ωi∗ = ω̃i∗ = ω̃j∗ = ω0 . Set ω̇ = 0 in
(C.1) and we have
X
ω0
(Di + τi ) = 0
(C.2)
i∈NCG

This implies that ω0 = 0 due to Di + τi > 0.
Other proofs are same with that in Theorem 2, which are
omitted.
D. Proof of Theorem 5
Proof of Theorem 5. We still use the Lyapunov function (22)
to analyze the stability of the closed-loop system (24). Denote
y = (ωiT , xT1 , xT2 )T , i ∈ NCG , and define the following
function


−Di ωi
 0



g
−
+

 − fi (Pi ) + µi − γi + γi




f
˜
µi
f (y) = 
 , i ∈ NCG


µ
−
µ
i
j




P gi − Pig
g
Pig − P i
(D.1)

JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. XX, XX XXXX

where fµi = Pig − p̂i −
0

P

(µi − µj ) −

j∈Nci

P

12

zij − τi µi −

j∈Nci

τi fi (Pig ) + βi ωi 1 .
The derivative of Wk is
X
Mi (ωi − ωi∗ )ω̇i + (x − x∗ )T · K −1 ẋ
Ẇk =
i∈NG

=

X

∗
(ωi − ωi∗ )(Pig − Pig∗ − Di (ωi − ωi∗ ) − (Pei − Pei
))

i∈NG

+ (x − x∗ )T · K −1 ẋ
X
∗
))
(ωi − ωi∗ )(Pig − Pig∗ − (Pei − Pei
≤
i∈NG

−

X

Di (ωi − ωi∗ )2 + (y − y ∗ )T f˜(y)

(D.2)

i∈NU G

where the inequality is due to the same reason for that in (B.4).
Divide Ẇk into two parts, and Ẇk1 = (y −y ∗ )T f˜(y), Ẇk2 =
Ẇk − Ẇk1 . Then we will analyze the sign of Ẇk1 .
Ẇk1 = (y − y ∗ )T f˜(y)
Z 1
∂
=
(y − y ∗ )T f˜(z̃(s))(y − y ∗ )ds + (y − y ∗ )T f˜(y ∗ )
∂ z̃
0

 T
Z
1 1
∂ ˜
∂ ˜
≤
f (z̃(s)) +
f (z̃(s)) (y − y ∗ )ds
(y − y ∗ )T
2 0
∂ z̃
∂ z̃
Z 1
=
(y − y ∗ )T [H(z̃(s))] (y − y ∗ )ds
(D.3)
0

where z̃(s) = y ∗ + s(y − y ∗ ).RThe second equation is from
1
the fact that f˜(y) − f˜(y ∗ ) = 0 ∂∂z̃ f˜(z̃(s))(y − y ∗ )ds. The
inequality is due to either f˜(y ∗ ) = 0 or f˜(y ∗ ) < 0, yi ≥ 0,
i.e. (y − y ∗ )T f˜(y ∗ ) ≤ 0.


D
0
0
0
0 0
 0
∇2P g f (P g )
I
0 −I I 


˜
 β τ ∇2 g f (P g ) − I τ + Lc C 0 0 
∂ f (y)
P


= −
0
−C T
0
0 0 
∂y
 0

 0
I
0
0
0 0 
0
−I
0
0
0 0
(D.4)
where D = diag(Di ), β = diag(βi ), τ = diag(τi ), I is the
identity matrix with dimension nCG , C is the incidence matrix
of the communication graph, Lc is the Laplacian matrix of the
communication graph.
Finally, H in (D.3) is


∂ ˜
1 ∂T ˜
f (y) +
f (y)
H=
2 ∂y
∂y


−D
0
− 12 β
0 0 0
 0
−∇2P g f (P g ) − τ2 ∇2P g f (P g ) 0 0 0 
 1

 − β − τ ∇2 g f (P g )
−τ − Lc
0 0 0 
2
2 P

=
 0
0
0
0 0 0 


 0
0
0
0 0 0 
0
0
0
0 0 0
(D.5)
1 Sometimes, we also use β instead of β (t) for simplification. In addition,
i
i
if (9) is not binding, we can omit γi− , γi+ in neighborhoods of equilibrium
points

H ≤ 0, if

−D
 0
− 21 β


0
− 12 β
−∇2P g f (P g ) − τ2 ∇2P g f (P g )  < 0
− τ2 ∇2P g f (P g )
−τ

(D.6)

By Schur complement [29], we only need

−1  1

 −Di

0
− 2 βi
−τi − − 21 βi − τ2i ci
<0
0
−ci
− τ2i ci
(D.7)
where ci = ∇2P g f (Pig ). Solving (D.7), we can get
i
p
p
− τi Di (4 − τi ci ) < βi < τi Di (4 − τi ci )

(D.8)

By A2, we know ci ≤ li , thus
p
p
τi Di (4 − τi li ) ≤ τi Di (4 − τi ci ),
p
p
− τi Di (4 − τi li ) ≥ − τi Di (4 − τi ci ).
Here we need τi > 0, 4 − τi li > 0, i.e. 0 < τi < 4/li . Finally,
we have
p
p
− τi Di (4 − τi li ) < βi < τi Di (4 − τi li )
(D.9)
p
i.e. β̄i < τi Di (4 − τi li ), implying (26).
Analysis of Ẇk2 and Ẇp as well as the convergence to a
point are same as those in the proof of Theorem 3, which are
omitted here.

