1

Distributed Frequency Control with Operational
Constraints, Part II: Network Power Balance

arXiv:1703.00083v1 [] 28 Feb 2017

Zhaojian Wang, Feng Liu, Steven H. Low, Fellow, IEEE, Changhong Zhao, and Shengwei Mei Fellow, IEEE

Abstract—In Part I of this paper we propose a decentralized
optimal frequency control of multi-area power system with operational constraints, where the tie-line powers remain unchanged
in the steady state and the power mismatch is balanced within
individual control areas. In Part II of the paper, we propose a
distributed controller for optimal frequency control in the network power balance case, where the power mismatch is balanced
over the whole system. With the proposed controller, the tieline powers remain within the acceptable range at equilibrium,
while the regulation capacity constraints are satisfied both at
equilibrium and during transient. It is revealed that the closedloop system with the proposed controller carries out primal-dual
updates with saturation for solving an associated optimization
problem. To cope with discontinuous dynamics of the closed-loop
system, we deploy the invariance principle for nonpathological
Lyapunov function to prove its asymptotic stability. Simulation
results are provided to show the effectiveness of our controller.
Index Terms—Power system dynamics, frequency control;
network power balance; distributed control.

I. I NTRODUCTION
In Part I of the paper we have investigated the optimal
frequency control of multi-area power system with operational
constraints [1]. In that case, the tie-line powers are required to
be unchanged in the steady state after load disturbances, which
implies the power mismatch in each area has to be balanced
individually. It is referred to as the per-node power balance
case. In Part II of the paper, we consider the transmission
congestion in the distributed optimal frequency design.
The per-node balance case in Part I mainly considers the
situation where the power delivered from one area to another is
fixed, e.g. contract power, which should not be violated in normal operation. However, in some circumstances, control areas
may cooperate for better frequency recovery or regulation cost
reduction. In this case, power mismatch may be balanced by all
generations and controllable loads among all control areas in
cooperation. Similar situations also appear in one control area
with multiple generators and controllable loads that cooperate
to eliminate the power mismatch in the area. It is referred
to as the network power balance case. Compared with the
This work was supported by the National Natural Science Foundation of
China ( No. 51377092, No. 51677100, No. 51621065), Foundation of Chinese
Scholarship Council (CSC No. 201506215034), 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.
Z. Wang, F. Liu and S. Mei are with the Department of Electrical Engineering, Tsinghua University, Beijing, China, 100084 e-mail:
(lfeng@tsinghua.edu.cn).
S. H. Low and C. Zhao are with the Department of Electrical Engineering, California Institute of Technology, Pasadena, CA, USA, 91105 email:(slow@caltech.edu)

per-node balance case, the most challenging problem in this
case is that the tie-line powers may change and congestions
may occur. In addition, local information is not sufficient
and neighboring information turns to be helpful. As for the
constraints, the tie-line power constraints are not hard limits,
which only need to be satisfied at equilibrium. The capacity
limits on the generations and controllable loads also required
to be satisfied both in steady state and during transient.
In the recent literature of frequency control, tie-line power
constraints are considered in [2]–[7]. In [2], tie-line power
constraints are included in the load-side secondary frequency
control. A virtual variable is used to estimate the tie-line
power, whose value is identical to the tie-line power at
equilibrium. In [3], [4], an optimal economic dispatch problem
including tie-line power constraints is formulated, then the
solution dynamics derived from a primal-dual algorithm is
shaped as a port-Hamiltonian form. The power system dynamics also have a port-Hamiltonian form, which are interconnected with the solution dynamics to constitute a closedloop Hamiltonian system. Then, the optimality and stability are
proved. In [6], a unified method is proposed for primary and
secondary frequency control, where the congestion management is implemented in the secondary control. In [7], a realtime control framework is proposed for tree power networks,
where transmission capacities are considered.
Similar to the per-node balance case, hard limits, such as
capacity constraints of power injections on buses, are enforced
only in the steady-state in the literature, which may fail if such
constraints are violated in transient. Here we construct a fully
distributed control to recover nominal frequency while eliminating congestion. Differing from the literature, it enforces
regulation capacity constraints not only at equilibrium, but
also during transient. We show that the controllers together
with the physical dynamics serve as primal-dual updates with
saturation for solving the optimization problem. The optimal
solution of the optimization problem and the equilibrium point
of closed-loop system are identical.
The enforcement of capacity constraints during transient and
tie-line power limits simultaneously makes the stability proof
difficult. Specifically, the Lyapunov function is not continuous
anymore, as in the per-node case in Part I of the paper. In
this situation, the conventional LaSalle’s invariance principle
does not apply. To overcome the difficulty, we construct a
nonpathological Lyapunov function to mitigate the impacts of
nonsmooth dynamics. The salient features of the controller
are:
1) Control goals: the controller restores the nominal frequency and balance the power mismatch in the whole

2

system after unknown load disturbance while minimizing the regulation costs;
2) Constraints: the regulation capacity constraints are always enforced even during transient and the congestions
can be eliminated automatically;
3) Communication: only neighborhood communication is
needed in the network balance case;
4) Measurement: the controller is adaptive to unknown
load disturbances automatically with no need of load
measurement.
The rest of this paper is organized as follows. In Section II, we describe our model. Section III formulates the
optimal frequency control problem in the network balance
case, presents the distributed frequency controller and proves
the optimality, uniqueness and stability of the closed-loop
equilibrium. Simulation results are given in Section IV. Section
V concludes the paper.

nominal or scheduled values so that, e.g., ω j (t) = 0 means the
frequency is at its nominal value.
In this paper, all nodes cooperate to rebalance power over
the entire network after a disturbance. The power flows Pi j on
the tie lines may deviate from their scheduled values and we
require that they satisfy line limits, i.e.,
∀(i, j) ∈ E

Pi j ≤ Pi j ≤ Pi j

for some upper and lower bounds Pi j , Pi j .
In DC approximation the power flow on line (i, j) is given
by Pi j = Bi j (θi − θ j ). Hence line flow constraints in the pernode balance case is θi = θ j for all (i, j) ∈ E and in the network
balance case is:
∀(i, j) ∈ E

θ i j ≤ θi − θ j ≤ θ i j

(3)

II. N ETWORK MODEL
We summarize the notation used in Part I [1]. The power
network is model by a directed graph G := (N, E) where N =
{0, 1, 2, ...n} is the set of nodes (control areas) and E ⊆ N ×
N is the set of edges (tie lines). If a pair of nodes i and j
are connected by a tie line directly, we denote the tie line
by (i, j) ∈ E. Let m := |E| denote the number of tie lines.
Use (i, j) ∈ E or i → j interchangeably to denote a directed
edge from i to j. Assume the graph is connected and node 0
is a reference node. For each node j ∈ N, θ j (t) denotes the
rotor angle at node j at time t and ω j (t) is the frequency.
Pjg (t) denotes the (aggregate) generation at node j at time t
and ugj (t) is its generation control command. Pjl (t) denotes
the (aggregate) controllable load and ulj (t) is its load control
command. p j is the (aggregate) uncontrollable load.
The power system dynamics for each node j ∈ N is
θ̇ j = ω j (t)

where θ i j = Pi j /Bi j , θ i j = Pi j /Bi j are lower and upper bounds
on angle differences.
As the generation Pjg and load Pjl in each area can increase
or decrease, and a line flow Pi j can in either direction, we
make the following assumption.
A1:

g

l

1) Pgj < 0 < P j and Plj < 0 < P j for ∀ j ∈ N.
2) θ i j ≤ 0 ≤ θ i j for (i, j) ∈ E.
3) θ0 (t) := 0 and φ0 (t) := 0 for all t ≥ 0.

Here φ is a vector variable that represents virtual phase angles
in the network balance case in Section III. The assumption
θ0 ≡ φ0 ≡ 0 amounts to using (θ0 (t), φ0 (t)) as reference angles.
It is made merely for notational convenience: as we will see,
the equilibrium point will be unique with this assumption (or
unique up to reference angles without this assumption).

(1a)

M j ω̇ j = Pjg (t) − Pjl (t) − p j − D j ω j (t)
Bi j (θi (t) − θ j (t)) −
B jk (θ j (t) − θk (t))
i:i→ j
k: j→k
= −Pjg (t) + ugj (t) − ω j (t)/R j
= −Pjl (t) + ulj (t)

+
T jg Ṗjg
T jl Ṗjl

∑

∑

(1b)
(1c)
(1d)

where D j > 0 are damping constants, R j > 0 are droop
parameters, and B jk > 0 are line parameters that depend on
the reactance of the line ( j, k). Let x := (θ , ω, Pg , Pl ) denote
the state of the network and u := (ug , ul ) denote the control.
The capacity constraints are:
Pgj

Pjg (t)

≤

Plj

≤ Pjl (t)

g
≤ Pj,
l
≤ P j,

III. C ONTROLLER D ESIGN FOR T HE N ETWORK P OWER
BALANCE C ASE

j∈N

(2a)

j∈N

(2b)

Here (2a) and (2b) are hard limits on the regulation capacities
of generation and controllable load at each node, which should
not be violated at any time even during transient.
The system operates in a steady state initially, i.e., the
generation and the load are balanced and the frequency is at
its nominal value. All variables represent deviations from their

In the per-node balance case, individual control areas rebalance power within their own areas after disturbances. However,
in many circumstances, it may be more efficient for all control
areas to eliminate power imbalance of the overall system in a
coordinated manner. This can be modeled as the condition:

∑ j Pjg

=

∑j



Pjl + p j

(4)

In this case the tie-line flows may not be restored to their predisturbance values. To ensure that they are within operational
limits, the constraints (3) are imposed.
Even though the philosophy of the controller design as well
as the proofs are similar to the per-node case, the details are
much more complicated. Our presentation will however be
brief where there is no confusion.

3

A. Control goals
In the network power balance case, the control goals are
formalized as the following optimization problem.
NBO: min

z2j
α j  2
β j  2
Dj
∑ 2 Pjg + ∑ 2 Pjl + ∑ 2 ω 2j + ∑ 2
j
j
j
j
(5a)
g

l

g

l

over

x := (θ , φ , ω, P , P ) and u := (u , u )

s. t.

(2),
Pjg = Pjl + p j +U j (θ , ω)
Pjg

= Pjl

+ p j + Û j (φ )

θ i j ≤ φi − φ j ≤ θ i j ,
Pjg
Pjl

=
=

ugj ,
ulj ,

j∈N

Even though neither the network balance condition (4) nor
the line limits (3) are explicitly enforced in (5), they are
satisfied at optimality (Theorem 2 below). Indeed, the virtual
phase angles φ and the conditions (5b)–(5d) are carefully
designed to enforce these conditions as well as to restore
the nominal frequency ω ∗ = 0 at optimality. This technique
is previously used in [2].

(5b)

j∈N

(5c)

(i, j) ∈ E

(5d)

j∈N

(5e)

j∈N

(5f)

where α j > 0, β j > 0 are constant weights; z j is a shorthand
defined for convenience as
zj

2) Even though any feasible point of (5) has z j = 0, the
objective function is augmented with z2j to improve
convergence (see [8]).

:= Pjg − Pjl − p j − Û j (φ )

U(θ , ω) := Dω +CBCT θ , and Û(φ ) := CBCT φ .
While θ represents the phase angles in the physical power
network, φ is a cyber quantity that can be interpreted as virtual
phase angles (see remarks below). The matrices D, C and B
are defined as in the previous section.
As in the per-node case, we define the variables θ̃i j := θi −
θ j , or in vector form, θ̃ := CT θ . As we fix θ0 := 0 to be
a reference angle under assumption A1, θ̃ = CT θ defines a
bijection between θ and θ̃ . Similarly we define φ̃i j := φi − φ j
or φ̃ := CT φ , and φ0 := 0 so there is a bijection between φ
and φ̃ . Note that both θ̃ and φ̃ are restricted to the column
space of CT . We will use (θ , φ ) and (θ̃ , φ̃ ) interchangeably.
For instance we will abuse notation and write Û(φ ) := CBCT φ
or Û(φ̃ ) := CBφ̃ .
We now summarize some of the interesting properties
of NBO (5) that will be proved formally in the next two
subsections. We first compare NBO (5) with PBO in [1] for
the per-node balance case.
Remark 1 (Comparison of NBO and PBO).
1)
Intuitively the network balance condition (4) is a
relaxation of the per-node balance condition (3) in [1],
and hence we expect that the optimal cost of NBO
lower bounds that of PBO. This is indeed the case,
as we now argue. Constraint (5c) implies that any
feasible point of (5) has z j = 0 and hence these two
optimization problems have the same objective function.
Their variables and constraints are different in that PBO
directly enforces the per-node balance condition while
NBO (5) has the additional variable φ and constraints
∗
(5c)(5d). Any optimal point (θ ∗ , ω ∗ , Pg∗ , Pl ) for PBO
however defines a feasible point (θ ∗ , φ , ω ∗ , Pg∗ , Pl∗ )
for NBO (5) with the same cost where φ = θ ∗ . The
point (θ ∗ , φ , ω ∗ , Pg∗ , Pl∗ ) satisfies (5c)(5d) because
∗
(θ ∗ , ω ∗ , Pg∗ , Pl ) satisfies (5b), ω ∗ = 0 and θi∗ = θ j∗ by
Theorem 2 in [1], and θ i j ≤ 0 ≤ θ i j by assumption A1.

Remark 2 (Virtual phase angles φ ).
1) Under mild conditions, ω ∗ = 0 at optimality for both PBO and NBO.
For NBO, this is a consequence of the constraint (5c)
on φ ; see Lemma A.1. Summing (5c) over all j ∈ N
also implies the network balance condition (4) since
1T Û(φ ) = 1T CBCT φ = 0.
2) In PBO, θi∗ = θ j∗ at optimality (i.e., tie-line flows are
restored Pi∗j = 0) and U(θ ∗ , ω ∗ ) = 0. This does not
necessarily hold in NBO. However φ is regarded as
virtual phase angles because, at optimality, φ ∗ differs
from the real phase angles θ ∗ only by a constant,
φ ∗ − θ ∗ = 1(φ0 − θ0 ) (Lemma A.1 in the appendix).
Hence CT φ ∗ − CT θ ∗ = CT · 1(φ0 − θ0 ) = 0, implying
φ̃i∗j = θ̃i∗j . Then the constraints (5d) are exactly the flow
constraints (3). In other words, we impose the flow
constraints on θ̃ indirectly by enforcing such constraints
on the virtual angle φ̃ . .

B. Distributed controller
Our control laws are:


λ̇ j = γ λj Pjg (t) − Pjl (t) − p j − Û j (φ̃ (t)) , j ∈ N

(6a)

η̇i+j

=

γiηj [φ̃i j (t) − θ i j ]+
η+

∀(i, j) ∈ E

(6b)

η̇i−j

= γiηj [θ i j − φ̃i j (t)]+
η−

∀(i, j) ∈ E

(6c)

φ̃˙i j

ij
ij

φ̃

= γi j (Bi j [λi (t) − λ j (t) + zi (t) − z j (t)]

+ ηi−j (t) − ηi+j (t)
∀(i, j) ∈ E
(6d)
h

iPgj
ugj (t) = Pjg (t) − γ gj α j Pjg (t) + ω j (t) + z j (t) + λ j (t)
g
Pj

+ ω j (t)/R j ,
j∈N
(6e)

iPlj
h
ulj (t) = Pjl (t) − γ lj β j Pjl (t) − ω j (t) − z j (t) − λ j (t)
l
Pj

j∈N

(6f)

where γ λj , γiηj , γi j , γ gj , γ lj are positive constants. For any xi , ai ∈
R, the operator [xi ]+
ai is defined by
φ̃

[xi ]+
ai


:=

xi if ai > 0 or xi > 0;
0, otherwise.

For a vector case, [x]+
a is defined accordingly componentwise [9].

4

Here we assume that each node i updates a set of internal
1
states (λi (t), ηi+j (t), η −
ji (t), φ̃i j (t)) according to (6a)–(6d). In
contrast to the completely decentralized control derived in the
per-node balance case, here the control is distributed where
each node i updates (λi (t), ηi+j (t), ηi−j (t)) using only local
measurements or computation but requires the information
(λ j (t), z j (t)) from its neighbors j to update φ̃i j (t). Note that
zi (t) is not a variable but a shorthand for (function) Pjg (t) −
Pjl (t) − p j − Û j (φ̃ (t)). The control inputs (ugi (t), uli (t)) in
(6e)(6f) are functions of the network state (Pig (t), Pil (t), ωi (t))
and the internal state (λi (t), ηi+j (t), ηi−j (t), φ̃i j (t)). We write ugj
and ulj as functions of (Pjg , Pjl , ω j , λ j ): for j ∈ N


(7a)
ugj (t) := ugj Pjg (t), ω j (t), λ j (t), z j (t)


ulj (t) := ulj Pjl (t), ω j (t), λ j (t), z j (t)
(7b)
where the functions are defined by the right-hand sides of
(6e)(6f).
Now we comment on the implementation of the control (6).
Remark 3 (Implementation).
1) As discussed above,
communication is needed only between neighboring
nodes (areas) to update the variables φ̃i j (t).
2) Similar to the per-node power balance case, we can
avoid measuring the load change p j by using (1b) and
the definition of z j (t) to replace (6a) with
z j (t) = M j ω̇ j + D j ω j (t) −

∑ Pi j (t) + ∑

i:i→ j

Pjk (t)

k: j→k

− Û j (φ̃ (t))
= γ λj z j (t)

λ̇ j

C. Design rationale
The controller design (6) is also motivated by a (partial)
primal-dual algorithm for (5), as for the per-node power
balance case.
Primal-dual algorithms. The optimization problem in the
network balance case differs from that in the per-node balance case in the inequalities (5d) on φ̃ . Consider a general
constrained convex optimization with inequality constraints:
min f (x) s.t. g(x) = 0, h(x) ≤ 0

where Γx , Γρ1 , Γρ2 are strictly positive diagonal gain matrices.
Here, if a is a scalar then (a)+ := max{a, 0} and if a is a
vector then (a)+ is defined accordingly componentwise. For a
dual algorithm, (8a) is replaced by
x(t) := min L(x; ρ(t))
x∈X

(8d)

As for the per-node balance case, all variables in x(t) are
updated according to (8a) except ω(t) which is updated
according to (8d), as we see below.
The set X in (8a) is defined by the constraints (2):
n
o
g l
(9)
X := (Pg , Pl ) : (Pg , Pl ) ≤ (Pg , Pl ) ≤ (P , P )
Controller (6) design. Let ρ1 := (λ , µ) be the Lagrange multipliers associated with constraints (5c) and (5b) respectively,
ρ2 := (η + , η − ) the multipliers associated with constraints
(5d), and ρ := (ρ1 , ρ2 ). Define the Lagrangian of (5) by (10).
Note that it is only a function of (x, ρ) and independent of
u := (ug , ul ) as we treat u as a function of (x, ρ) defined by
the right-hand sides of (6e)(6f).
The closed-loop dynamics (1)(6) carry out an approximate
primal-dual algorithm (8) for solving (5) in real time over the
coupled physical power network and cyber computation. Since
the reasoning is similar to the per-node balance case, we only
provide a summary. Rewrite the Lagrangian L2 in vector form

1 g T g g
L2 (x; ρ) =
(P ) A P + (Pl )T Al Pl + ω T Dω + zT z
2


+ λ T Pg − Pl − p −CBφ̃


+ µ T Pg − Pl − p − Dω −CBθ̃


+ (η + )T φ̃ − θ + (η − )T θ − φ̃
(11)
where Al := diag(β j , j ∈ N), B := diag(Bi j , (i, j) ∈ E).
First, the control (6b)(6c) can be interpreted as a continuoustime version of the dual update (8c) on the dual variable ρ2 :=
(η + (t), η − (t)):

+
η̇ + = Γη ∇η + L2 (x(t); ρ(t)) η + (t)
+
− 
η̇ − = Γη ∇η − L2 (x(t); ρ(t)) η −

x∈X

Rn

where f :
→ R, g : Rn → Rk1 , h : Rn → Rk2 , and X ⊆ Rn
is closed and convex. Here an inequality constraint h(x) ≤ 0
is imposed explicitly. Let ρ1 ∈ Rk1 be the Lagrange multiplier
associated with the equality constraint g(x) = 0, ρ2 ∈ Rk2 that
associated with the inequality constraint h(x) ≤ 0, and ρ :=
(ρ1 , ρ2 ). Define the Lagrangian L(x; ρ) := f (x) + ρ1T g(x) +
ρ2T h(x). A standard primal-dual algorithm takes the form:
x(t + 1) := ProjX (x(t) − Γx ∇x L(x(t); ρ(t))) (8a)
ρ1 (t + 1) := ρ1 (t) + Γρ1 ∇ρ1 L(x(t); ρ(t))
ρ2 (t + 1) :=
1 For

+
ρ2 (t) + Γ ∇ρ2 L(x(t); ρ(t))
ρ2

(8b)

where Γη := diag(γiηj , (i, j) ∈ E).
Second, the control (6a) carries out the dual update (8b) on
λ (t):
λ̇

= Γλ ∇λ L2 (x(t), ρ(t))

(12a)

where Γλ := diag(γ λj , j ∈ N). The swing dynamic (1b) carries
out the dual update (8b) on µ(t) because, as in the per-node
balance case, we can identify µ(t) ≡ ω(t) so that
µ̇

= ω̇

= M −1 ∇µ L2 (x(t); ρ(t))

(12b)

(8c)

each (directed) link (i, j) ∈ E we assume that only node i maintains
the variables (ηi+j (t), η −
ji (t), φ̃i j (t)). In practice, node j will probably maintain
symmetric variables to reduce communication burden or for other reasons
outside our mathematical model here.

where M := diag(M j , j ∈ N).
Finally we show that (1a), (1c), (1d), and (6d) implement
a mix of the primal updates (8a) and (8d) on the primal variables x := (θ̃ (t); φ̃ (t); ω(t); Pg (t); Pl (t)). Setting ω(t) ≡ µ(t)

5

 2
∑ α j Pjg +

1
L2 (x; ρ) =
2

j∈N

 2
∑ β j Pjl +

j∈N

!

∑ D j ω 2j + ∑ z2j

j∈N

!
+

j

∑ µj

Pjg − Pjl

− p j − D jω j +

j∈N

∑

Bi j θi j −

i:i→ j

∑

B jk θ jk

k: j→k

!
+

∑ λj

Pjg − Pjl

− pj +

j∈N

∑

Bi j φi j −

i:i→ j

∑

B jk φ jk

k: j→k

∑


ηi−j θ i j − φi j (t) +

g

(12c)

ηi+j φi j (t) − θ i j

∑



(10)

(i, j)∈E

(i, j)∈E

is equivalent to the primal update (8d) on ω(t), as in the pernode balance case. Moreover the control laws (6e)(6f) are then
equivalent to
T g Ṗg = [Pg (t) − Γg ∇Pg L2 (x(t), ρ(t))]PPg − Pg (t)
h
iPl
T l Ṗl = Pl (t) − Γl ∇Pl L2 (x(t), ρ(t)) l − Pl (t)

+

2) The nominal frequency is restored, i.e., ω ∗j = 0 for all
j ∈ N; moreover φ̃i∗j = θ̃i∗j for all (i, j) ∈ E.
3) The network balance condition (4) is satisfied by x∗ .
4) The line limits (3) are satisfied by x∗ , implying Pi j ≤
Pi j ≤ Pi j on every tie line (i, j) ∈ E.

where Γφ̃ := diag(γi j , (i, j) ∈ E) .

Theorem 2 shows that the equilibrium point has a simple
yet intuitive structure. Moreover, Theorem 2 implies that the
closed-loop system can autonomously eliminate congestions
on tie lines. This feature has important implications. It means
our distributed frequency control is capable of serving as a
corrective re-dispatch without the coordination of dispatch
centers if a congestion arises. This can enlarge the feasible
region for economic dispatch, since corrective re-dispatch has
been naturally taken into account.
The proofs of Theorem 1 and 2 are given in Appendix A.

D. Optimality of equilibrium point

E. Asymptotic stability

P

(12d)

i.e., the generator and controllable load at each node j carry
out the primal update (8a). For (θ̃ , φ̃ ), (1a) and (6d) are
equivalent to the primal update (8a):
θ̃˙
φ̃˙

= −B−1 ∇θ̃ L2 (x(t), ρ(t))
φ̃

= −Γ ∇φ̃ L2 (x(t), ρ(t))

(12e)
(12f)

φ̃

In this subsection, we address the optimality of the equilibrium point of the closed-loop system (1)(6). Given an
(x, ρ) := (θ̃ , φ̃ , ω, Pg , Pl ), (λ , µ), (η − , η + ) , recall that the
control input u(x, ρ1 , ρ2 ) is given by (7).
Definition 1. A point (x∗ , ρ ∗ ) := (θ̃ ∗ , φ̃ ∗ , ω ∗ , Pg∗ , Pl∗ , λ ∗ , η +∗ ,
η −∗ , µ ∗ ) is an equilibrium point or an equilibrium of the
closed-loop system (1)(6) if
1) The right-hand side of (1) vanishes at x∗ and u(x∗ , ρ ∗ ).
2) The right-hand side of (6a)–(6d) vanishes at (x∗ , ρ ∗ ).
Definition 2. A point (x∗ , ρ ∗ ) is primal-dual optimal if
(x∗ , u(x∗ , ρ ∗ )) is optimal for (5) and ρ ∗ is optimal for its dual
problem.
We make the following assumption:
A2: The problem (5) is feasible.
The following theorem characterizes the correspondence
between the equilibrium of the closed-loop system (1)(6) and
the primal-dual optimal solution of (5).
Theorem 1. Suppose A2 holds. A point (x∗ , ρ ∗ ) is primal-dual
optimal if and only if (x∗ , ρ ∗ ) is an equilibrium of closed-loop
system (1)(6) satisfying (2) and µ ∗ = 0.
Next result says that, at equilibrium, the network balance
condition (4) and line limits (3) are satisfied and the nominal
frequency is restored. Moreover the equilibrium is unique.
(x∗ , ρ ∗ )

Theorem 2. Suppose A1 and A2 hold. Let
be primaldual optimal. Then
1) The equilibrium (x∗ , µ ∗ ) is unique, with (θ ∗ , φ ∗ ) being
unique up to (equilibrium) reference angles (θ0 , φ0 ).

In this subsection, we address the asymptotic stability of the
closed-loop system (1)(6), under an additional assumption:
A3: The initial state of the closed-loop system (1)(6) is finite,
and pgj (0), plj (0) satisfy constraint (2).
As in the per-node balance case the closed-loop system
(1)(6) satisfies constraint (2) even during transient.
Lemma 3. Suppose A1 and A3 hold. Then constraint (2) is
satisfied for all t > 0, i.e. (Pg (t), Pl (t)) ∈ X for all t ≥ 0 where
X is defined in (9).
The proof is exactly the same as that for Lemma 3 in [1]
and omitted.
Similar to the per-node balance case, we first rewrite the
closed-loop system using states θ̃ , φ̃ instead of θ , φ (they are
equivalent under assumption A1). Setting µ ≡ ω, the closedloop system (1)(6) is equivalent to (in vector form):
θ̃˙ (t) = CT ω(t)
(13a)


−1
g
l
ω̇(t) = M
P (t) − P (t) − p − Dω(t) −CBθ̃ (t)
(13b)
g

g −1

l

l −1

Ṗ (t) = (T )
Ṗ (t) = (T )

g

g

(−P (t) + û (t))


−Pl (t) + ûl (t)

η̇ + (t) = Γη [φ̃ (t) − θ ]+
η+
−

η

− φ̃ (t)]+
η−

η̇ (t) = Γ [θ


λ̇ (t) = Γλ Pg (t) − Pl (t) − p −CBφ̃ (t)
φ̃˙ (t) = Γφ̃

(13c)
(13d)
(13e)
(13f)

(13g)

BCT λ (t) + BCT z(t) + η − (t) − η + (t) (13h)

6

where z(t) := Pg (t) − Pl (t) − p −CBφ̃ (t) and
iPg
h

ûg (t) = Pg (t) − Γg Ag Pjg (t) + ω(t) + z(t) + λ (t)
g
P

h

iPl
ûlj (t) = Pl (t) − Γl Al Pjl (t) − ω(t) − z(t) − λ (t)
l
P

:= (θ̃ , ω, Pg , Pl , λ , η + , η − , φ̃ ).

Denote w
Note that the right-hand sides of (13e)(13f) are discontinuous due to projection to the nonnegative quadrant for
(η + (t), η − (t)). The system (13) is called a projected dynamical system and we adopt the concept of Caratheodory
solutions for such a system where a trajectory (w(t),t ≥ 0)
is called a Caratheodory solution, or just a solution, to (13)
if it is absolutely continuous in t and satisfies (13) almost
everywhere. The result in [10, Theorems 2 and 3] implies that,
given any initial state, there exists a unique solution trajectory
to the closed-loop system (13) as the unprojected system is
Lipschitz and the nonnegative quadrant is closed and convex.
See [11, Theorem 3.1] for extension of this result to the Hilbert
space.
With regard to system (13), we first define two sets, σ + and
−
σ , as follows [8].
σ+

:= {(i, j) ∈ E | ηi+j = 0, φ̃i j − θ i j < 0}

σ−

:= {(i, j) ∈ E | ηi−j = 0, θ i j − φ̃i j < 0}

Then (6b) and (6c) are equivalent to
 η
γi j (φ̃i j − θ i j ),
η̇i+j =
0,
 η
γi j (θ i j − φ̃i j ),
η̇i−j =
0,

if(i, j) ∈
/ σ +;
if(i, j) ∈ σ + .

(14a)

if(i, j) ∈
/ σ −;
if(i, j) ∈ σ − .

(14b)

In a fixed σ + , σ − , define F(w).

−B1/2CT ω

 −M −1/2 Pg − Pl − p − Dω −CBθ̃

 (T g )−1 (Ag Pg + ω + z + λ )


 (T l )−1 Al Pl − ω − z − λ

F(w) =  −(Γη )1/2 [φ̃ − θ ]+

η+

+
η
1/2
 −(Γ ) [θ − φ̃ ]η −


 −(Γλ )1/2 Pg − Pl − p −CBφ̃

−(Γφ̃ )1/2 BCT λ + BCT z + η − − η +








(15)






If σ + and σ − do not change, F(w) is continuously differentiable in w.
Similarly, we define S := Rm+n+1 × X × R2m+n+1+m , where
the closed convex set X is defined in (9). Then for any w we
define the projection of w − F(w) onto S as
H(w) := ProjS (w − F(w)) := arg min ky − (w − F(w))k2
y∈S

Then the closed-loop system (13) is equivalent to
ẇ(t) = Γ2 (H(w(t)) − w(t))
where the positive definite gain matrix is :

Γ2 = diag B−1/2 , M −1/2 , (T g )−1 , (T l )−1 ,

λ 1/2
η 1/2
η 1/2
φ̃ 1/2
(Γ ) , (Γ ) , (Γ ) , (Γ )

(16)

Note that the projection operation H has an effect only on
(Pg ; Pl ) and Lemma 3 indicates that w(t) ∈ S for all t > 0,
justifying the equivalence of (13) and (16).
A point w∗ ∈ S is an equilibrium of the closed-loop system
(16) if and only if it is a fixed point of the projection H(w∗ ) =
w∗ . Let E2 := { w ∈ S | H(w(t)) − w(t) = 0 } be the set of
equilibrium points. Then we have the following theorem.
Theorem 4. Suppose A1, A2 and A3 hold. Starting from any
initial point w(0), w(t) remains in a bounded set for all t and
w(t) → w∗ as t → ∞ for some equilibrium w∗ ∈ E2 that is
optimal for problem (5).
For any equilibrium point w∗ , we define the following
function taking the same form as the per-node case.
1
Ṽ2 (w) = −(H(w) − w)T F(w) − ||H(w) − w||22
2
1
∗ T −2
+ k(w − w ) Γ2 (w − w∗ )
2

(17)

where k is small enough such that Γ2 − kΓ−1
2 > 0 is strictly
positive definite.
For any fixed σ + and σ − , Ṽ2 is continuously differentiable
as F(w) is continuously differentiable in this situation. Similar
to V1 (w) used in Part I of the paper, we know Ṽ2 (w) ≥ 0 on
S and Ṽ2 (w) = 0 holds only at any equilibrium w∗ = H(w∗ )
[12]. Moreover, Ṽ2 is nonincreasing for fixed σ + and σ − , as
we prove in Appendix B.
It is worth to note that the index sets σ + and σ − may change
sometimes, resulting in discontinuity of Ṽ2 (w). To circumvent
such an issue, we slightly modify the definition of V2 (w) at
the discontinuous points as:
1) V2 (w) := Ṽ2 (w), if Ṽ2 (w) is continuous at w;
2) V2 (w) := lim sup Ṽ2 (v), if Ṽ2 (w) is discontinuous at w.
v→w

Then V2 (w) is upper semi-continuous in w, and V2 (w) ≥ 0 on
S and V2 (w) = 0 holds only at any equilibrium w∗ = H(w∗ ).
As V2 (w) is not differentiable for w at discontinuous points,
we use the Clarke gradient as the gradient at these points [13].
Note that Ṽ2 is continuous almost everywhere except the
switching points. Hence both V2 (w) is nonpathological [14],
[15]. With these definitions and notations above, we can prove
Theorem 4. The detail of proof is provided in Appendix B.
IV. C ASE STUDIES
A. System configuration
A four-area system based on Kundur’s four-machine, twoarea system [16] [17] is used to test our optimal frequency
controller. There are one (aggregate) generator (Gen1∼Gen4),
one controllable (aggregate) load (L1c∼L4c) and one uncontrollable (aggregate) load (L1∼L4) in each area, which is
shown in Fig.1. The parameters of generators and controllable
loads are given in Table I. The total uncontrollable load in
each area are identically 480MW. At time t = 10s, we add
step changes on the uncontrollable loads in four areas to test
the performance of our controllers.
All the simulations are implemented in PSCAD [18] with
8GB memory and 2.39 GHz CPU. We use the detailed
electromagnetic transient model of three-phase synchronous

Bus7

Bus5

Bus8

Bus9

Bus10 Bus11 Bus3

T3
machines to simulate generators with both governors
and
exciters.
The uncontrollable load L1-L4 are modelled
Gen1
Gen3 by
the fixed load
in Bus6
PSCAD, while controllable load
L1c-L4c
T2
T4
Load1
Load2
are formulated
by the
self-defined controlled
current
Bus2
Bus4 source.
Gen4 We need
The closed-loop Gen2
system diagram is shown in Fig.2.
measure loacal frequency, generation, controllable load and
tie-line power flows to compute control demands. Only φ̃i j
are exchanged between neighbors. All variables are added by
their initial steady state values to explicitly show the actual
values.
T1

Bus1

Bus5

Bus7

Bus6

Bus8

Bus9

Bus10 Bus11 Bus2

T1

T2

Gen1

1) Stability and optimality: The dynamics of local frequencies and tie-line power flows are illustrated in Fig.3.
The frequencies are well restored in all four control areas
while the tie line powers are remained within their acceptable
ranges. The generations and controllable loads are different
from that before disturbance, indicating that the system is
stabilized at a new steady state. The resulting equilibrium
point is given in Table IV, which is identical to the optimal
solution of (5) computed by centralized optimization using
CVX. These simulation results confirm that our controller
can autonomously guarantee the frequency stability while
achieving optimal operating point in the overall system.

Gen2
L1c

L2

L1

Area2

L2c

Bus15 Bus16 Bus17 Bus4

Bus3 Bus14 Bus13 Bus12

T4

T3

Gen4

Gen3
Area3

L3c

L3

Area4

L4c

L4

60

50.2

Frequency (Hz)

Area1

80

50.3
Area 1
Area 2
Area 3
Area 4

50.1
50
49.9
49.8
49.7
0

Fig. 1: Four-area power system

Tie line power (MW)

Bus1

7

40
20
0

Line 21
Line 31
Line 32
Line 42

-20
-40

-60
20

40

60

80

100

120 140

160

-80
0

20

40

Time (s)

ij

Controller

Controller

u lj

u gj

uig

TABLE I:
S YSTEMPARAMETERS
, P g , Pl , P , P

, P , P , Pij , Pjk

u

l
j

i

Area j
1
2
3
4

Controller
g
j

u lj

Dj
0.04
0.045
0.05
0.055

Rj
0.04
0.06
0.05
0.045

i

i

ij

100

120 140

160

Pjg∗ (MW)
Pjl∗ (MW)

Area 1
620
23.6

Area 2
596
59.8

Area 3
660
23.6

Area 4
580
39.7

Tie line
Pi∗j (MW)

(2,1)
-39.94

(3,1)
13.35

(3,2)
53.27

(4,2)
-59.6

2) Dynamic performance: In this subsection,
 j , Pjgwe
, Pjlanalyze
, Pij , Pjk
the impacts of operational (capacity and line power) constraints on the dynamic property. Similarly, we compare the
responses of frequency controllers with and without
considerController
ing input saturations. The trajectories of mechanical power of
g
l
turbines and controllable loads are shown in uFig.4
and Fig.5, u j
j
respectively. In this case, the system frequency is restored, and
the same optimal equilibrium point is achieved.

ki

T jg

αj
βj
T jl
2
2.5
4
4
2.5
4
6
5
Controller
1.5
2.5
5
4
3
5.5 l 5
g 3

ui

80

Time (s)

TABLE III:
E QUILIBRIUM P OINTS

uil

Fig. 2: Closed-loop system diagram

g
j

60

Fig. 3: Dynamics of frequency (left) and tie-line flows (right)
in network balance case

i , Pi g , Pi l , Pij , Pki

 j , Pjg , Pjl , Pij , Pjk

ui

B. Network Power Balance Case
In this case, the generations in each area are initiated as
(560.9, 548.7, 581.2, 540.6) MW and the controllable loads
(70.8, 89.6, 71.3, 79.4) MW. The load changes are identical to
those in Table II, which are also unknown to the controllers.
We use method in Remark 4 to estimate the load changes.
Operational constraints on generations, controllable loads and
tie lines are shown in Table II.

900

Mechanical pwoer (MW)

j

DFIG
Bus12
Tw

850
Area 1
Area 2
Area 3
Area 4

800
750

Area 1
Area 2
Area 3
Area 4

700
650
600
550

500

TABLE II:

450
0

C APACITY LIMITS IN NETWORK CASE
Area 2
[530, 680]

Area 3
[550, 700]

Area 4
[530, 670]

[Plj ,P j ] (MW)

[20, 80]

[60, 100]

[20, 80]

[35, 80]

Tie line
[Pi j ,Pi j ] (MW)

(2,1)
[-65, 65]

(3,1)
[-65, 65]

(3,2)
[-65, 65]

(4,2)
[-65, 65]

l

40

60

80

Time (s)

Area 1
[550, 710]

g

[Pgj ,P j ] (MW)

20

100

120 140 160 0

20

40

60

80

100

120 140 160

Time (s)

Fig. 4: Mechanical outputs with(left)/without(right) capacity
constraints
3) Congestion analysis: In this scenario, we reduce tieline power constraints to Pi j = −Pi j = 50MW, which causes
congestions in tie-line (2,3) and (2,4). The steady states



8

120
Area 1
Area 2
Area 3
Area 4

Pl (MW)

100

Area 1
Area 2
Area 3
Area 4

80
60
40
20

0

0

20

40

60

80

100

120

140

160 0

20

40

Time (s)

60

80

100

120

140

160

Time (s)

Fig. 5: Controllable loads with(left)/without(right) capacity
constraints
under the distributed control are listed in Table IV. Note that
∗ +
Pjg∗ − p j + Pjl∗ − ∑ Pjk
∑ Pi∗j = 0 hold in each area.
i:i→ j

k: j→k

The dynamics of tie-line powers in two different scenarios
shown in Fig.6 indicate that (2,4) reaches the limit in steady
state. However, by adopting the proposed fully distributed
optimal frequency control, the congestion is eliminated and all
the tie line powers remain within the limits. Thus, congestion
control is achieved optimally in a distributed manner.
TABLE IV:
S IMULATION R ESULTS WITH C ONGESTION
Pjg∗ (MW)
Pjl∗ (MW)

Area 1
618
25.1

Area 2
595
60.7

Area 3
658
25.1

Area 4
585
34.9

Tie line
Pi∗j (MW)

(2,1)
-36.4

(3,1)
13.1

(3,2)
49.5

(4,2)
-49.9

80

Tie line power (MW)

60
40

20
Line 21
Line 31
Line 32
Line 42

0
-20

Line 21
Line 31
Line 32
Line 42

-40
-60

-80
0

20

40

60

80

Time (s)

100

120 140

160 0

20

40

60

80

100

120 140

160

Time (s)

Fig. 6: Tie line power with(left)/without(right) capacity constraints

V. C ONCLUDING R EMARKS
In this paper, we have devised a distributed optimal frequency control in the network balance case, which can autonomously restore the nominal frequencies after unknown
load disturbances while minimizing the regulation costs. The
capacity constraints on the generations and controllable loads
can also be satisfied even during transient. In addition, congestions can be eliminated automatically, implying tie-line powers
can be remained within given ranges. Only neighborhood
communication is required in this case. Like the per-node
case, here the closed-loop system again carries out a primaldual algorithm with saturation to solve the associated optimal
problem. To cope with the discontinuity introduced due to

enforcing different types of constraints, we have constructed
a nonpathological Lyapunov function to prove the asymptotically stability of the closed-loop systems. Simulations on a
modified Kundur’s power system validate the effectiveness of
our controller.
This approach is also applicable to other problem involving
frequency regulation, e.g. standalone microgrid or demand
side management. We highlight two crucial implications of
our work: First, our distributed frequency control is capable
of serving as an automatically corrective re-dispatch without
the coordination of dispatch center when certain congestion
happens; Second, the feasible region of economic decisions
can be enlarged benefiting from the corrective re-dispatch.
In this sense, our work may provide a systematic way to
bridge the gap between the (secondary) frequency control in a
fast timescale and the economic dispatch in a slow timescale,
hence breaking the traditional hierarchy of the power system
frequency control and economic dispatch.
R EFERENCES
[1] Z. Wang, F. Liu, S. H. Low, C. Zhao, and S. Mei, “Distributed frequency
control with operational constraints, part i: Per-node power balance,”
https://arxiv.org/abs/1702.07965, 2017.
[2] E. Mallada, C. Zhao, and S. H. Low, “Optimal load-side control for
frequency regulation in smart grids,” IEEE Trans. Autom. Control, to
appear, 2017.
[3] T. Stegink, C. D. Persis, and A. van der Schaft, “A unifying energybased approach to stability of power grids with market dynamics,” arXiv
preprint arXiv:1604.05200, 2016.
[4] T. Stegink, C. D. Persis, and A. van der Schaft, “A port-hamiltonian
approach to optimal frequency regulation in power grids,” in Proc. 54th
IEEE Conference on Decision and Control (CDC), Dec 2015, pp. 3224–
3229.
[5] P. Yi, Y. Hong, and F. Liu, “Distributed gradient algorithm for constrained optimization with application to load sharing in power systems,”
Syst. Control Lett., vol. 83, pp. 45–52, Sept 2015.
[6] 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.
[7] X. Zhang and A. Papachristodoulou, “A real-time control framework for
smart power networks: Design methodology and stability,” Automatica,
vol. 58, pp. 43–50, 2015.
[8] 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.
[9] A. Cherukuri, E. Mallada, and J. Corts, “Asymptotic convergence of
constrained primaldual dynamics,” Syst. Control Lett., vol. 87, pp. 10 –
15, 2016.
[10] P. Dupuis and A. Nagurney, “Dynamical systems and variational inequalities,” Annals of Operations Research, vol. 44, pp. 9–42, 1993.
[11] M.-G. Cojocaru and L. B. Jonker, “Existence of solutions to projected
differential equations in hilbert spaces,” Proc. Amer. Math. Soc., vol.
132, pp. 183–193, 2004.
[12] M. Fukushima, “Equivalent differentiable optimization problems and
descent methods for asymmetric variational inequality problems,” Math.
programming, vol. 53, no. 1-3, pp. 99–110, Jan 1992.
[13] F. H. Clarke, Optimization and nonsmooth analysis. Siam, 1990, vol. 5.
[14] A. Bacciotti and F. Ceragioli, “Nonpathological lyapunov functions and
discontinuous carathodory systems,” Automatica, vol. 42, no. 3, pp. 453
– 458, 2006.
[15] A. Bacciotti and F. Ceragioli, “Stability and stabilization of discontinuous systems and nonsmooth lyapunov functions,” ESAIM: Control,
Optimisation and Calculus of Variations, vol. 4, pp. 361–376, 1999.
[16] J. Fang, W. Yao, Z. Chen, J. Wen, and S. Cheng, “Design of anti-windup
compensator for energy storage-based damping controller to enhance
power system stability,” IEEE Trans. Power Syst., vol. 29, no. 3, pp.
1175–1185, May 2014.
[17] P. Kundur, Power System Stability and Control. McGraw-hill New
York, 1994, vol. 7.
[18] https://hvdc.ca/pscad/.

9

Since (x∗ , ρ ∗ ) is a saddle point, we must have
∇φ̃ L2 (x∗ , ρ ∗ ) = 0, implying φ̃˙ = 0.

A PPENDIX A
P ROOFS OF T HEOREM 1 AND T HEOREM 2
We start with a lemma.
Lemma A.1. Suppose (x∗ , u∗ ) is optimal for (5). Then
1) ω ∗ = 0, i.e., the nominal frequency is restored;
2) the network balance condition (4) is satisfied by x∗ ;
3) φ ∗ − θ ∗ = (φ0∗ − θ0∗ ) 1;
4) θ i j ≤ θi∗j ≤ θ i j , i.e., the line limits (3) are satisfied.
Proof. Suppose (x∗ , u∗ ) is optimal but ω ∗ 6= 0. Then (5c)
implies
Pg∗ − Pl∗ − p

= CBCT φ ∗

(A.1)

Consider x̂ := (θ̂ , φ ∗ , ω̂, Pg∗ , Pl∗ ) with θ̂ := φ ∗ and ω̂ := 0.
Then x̂ satisfies (5b)(5c) due to (A.1). Hence (x̂, u∗ ) is feasible
for (5) but has a strictly lower cost, contradicting the optimality
of (x∗ , u∗ ). Hence ω ∗ = 0.
Multiplier both sides of (A.1) by 1T yields the network
balance condition (4), proving 2).
To prove 3), setting ω ∗ = 0 in (5b) and combining with (5c)
yield
CBCT θ ∗

= Pg∗ − Pl∗ − p = CBCT φ ∗

Since CBCT is an (n + 1) × (n + 1) matrix with rank n, its null
space has dimension 1 and is spanned by 1 because CT 1 = 0.
Hence CBCT (φ ∗ − θ ∗ ) = 0 implies that φ ∗ − θ ∗ = (φ0∗ − θ0∗ ) 1.
To prove 4), note that φ̃ = CT φ and θ̃ = CT θ and hence
φ̃ ∗ − θ̃ ∗

We have the following result.
Lemma A.2. Suppose (x∗ , ρ ∗ ) is primal-dual optimal. Then
h

iPgj
g∗
g∗
g
g∗
∗
∗
∗
ug∗
=
P
=
P
−
γ
α
P
+
ω
+
z
+
λ
j j
j
j
j
j
j
j
j
g
Pj

=

⇐: Suppose now (x∗ , ρ ∗ ) is an equilibrium of the closed-loop
system (1)(6) that satisfies (2) and µ ∗ = 0. Since (5) is convex
with linear constraints, (x∗ , ρ ∗ ) is a primal-dual optimal if and
only if (x∗ , u(x∗ , ρ ∗ )) is primal feasible and satisfies
x∗ = arg min{L2 (x; ρ ∗ )|(x, u(x, ρ ∗ )) satisfies(2), (5e), (5f)}
x

(A.2)
This is because ∇ρ1 L2 (x∗ , ρ ∗ ) = 0 since µ̇ = λ̇ = 0, η +∗ ≥
0, η −∗ ≥ 0, and the complementary slackness condition
−∗
∗
+
∗
ηi+∗
j (φ̃i j − θ i j ) = 0, ηi j (θ i j − φ̃i j ) = 0 is satisfied since η̇ =
−
η̇ = 0.
To show that (x∗ , u(x∗ , ρ ∗ )) is primal feasible, note that
since (x∗ , u(x∗ , ρ ∗ )) is an equilibrium of (1), it satisfies ω ∗ = 0
and hence (5e)(5f), in addition to (2). Moreover ω̇ = 0 means
(x∗ , u(x∗ , ρ ∗ )) satisfies (5b), λ̇ = 0 implies (5c), η̇ + = η̇ − = 0
implies (5d).
To show that (x∗ , ρ ∗ ) satisfies (A.2), note that (5e)(5f) and
(6e)(6f) imply that

Pjl∗

h

iPlj
= Pjl∗ − γ lj β j Pjl∗ − ω ∗j − z∗j − λ j∗
l
Pj

for any γ gj > 0 and γ lj > 0.
Lemma A.2 shows that the saturation of control input does
not impact the optimal solution of optimization problem (5).
With Lemma 3, Lemma A.1 and Lemma A.2, we now can
prove Theorem 1 and Theorem 2.
Proof of Theorem 1. ⇒: Suppose (x∗ , ρ ∗ ) is primal-dual optimal. Then x∗ satisfies the operational constraints (2). Moreover
the right-hand side of (1) vanishes because:
• θ̇ = 0 since ω ∗ = 0 from Lemma A.1.
• ω̇ = 0 since constraint (5b) holds for x∗ .
• Ṗg = Ṗl = 0 since ω ∗ = 0 and x∗ satisfies (5e) and (5f).
• λ̇ = 0 since (5c) holds for x∗ .
• η̇ + = η̇ − = 0 since (5d) holds for x∗ .
• From (12f) we have
φ̃˙

= −Γφ̃ ∇φ̃ L2 (x∗ , ρ ∗ )

h

iPgj
Pjg∗ − γ gj α j Pjg∗ + ω ∗j + z∗j + λ j∗
g

Pjg∗

=

Pjl∗

h

iPlj
= Pjl∗ − γ lj β j Pjl∗ − ω ∗j − z∗j − λ j∗
l

= CT (φ ∗ − θ ∗ ) = (φ0∗ − θ0∗ )CT 1 = 0

i.e. φ̃ ∗ = θ̃ ∗ . We conclude from (5d) that θ i j ≤ θi∗j ≤ θ i j . This
completes the proof.

ul∗
j

Hence (x∗ , ρ ∗ ) is an equilibrium of the closed-loop system
(1)(6) that satisfies the operational constraints (2). Moreover
µ ∗ = ω ∗ = 0 since ∂∂ ωL2j (x∗ , ρ ∗ ) = D j (ω ∗j − µ ∗j ) = 0 and D j > 0
for all j ∈ N.

Pj

Pj

The rest of the proof follows the same line of argument as that
in Theorem 1 in [1]. This proves that (x∗ , ρ ∗ ) is primal-dual
optimal and completes the proof of Theorem 1.
Next we prove Theorem 2.
Proof of Theorem 2. Let (x∗ , ρ ∗ ) = (θ̃ ∗ , φ̃ ∗ , ω ∗ , Pg∗ , Pl∗ ,
λ ∗ ,η +∗ , η −∗ , µ ∗ ) be primal-dual optimal. For the uniqueness
of x∗ , (ω ∗ , Pg∗ , Pl∗ ) are unique because the objective function
in (5) is strictly convex in (ω, Pg , Pl ). Hence µ ∗ = ω ∗ is unique
as well. Assumption A1 that φ0∗ := 0 and (A.1) imply that φ ∗
is uniquely determined by the equilibrium (ω ∗ , Pg∗ , Pl∗ ). Since
θ ∗ −φ ∗ = (θ0∗ −φ0∗ )1, assumption A1 that θ0∗ := 0 then implies
that θ ∗ is unique. This proves the uniqueness of (x∗ , µ ∗ ).
The remaining three parts of the theorem follow from
Lemma A.1.
A PPENDIX B
P ROOF T HEOREM 4
We start with a lemma.
Lemma B.1. Suppose A1, A4 and A5 hold. Then
1) V̇2 (w(t)) ≤ 0 in a fixed σ + and σ − .
2) The trajectory w(t) is bounded, i.e., there exists w such
that kw(t)k ≤ w for all t ≥ 0.

10







Q= 






−BCT
D
I|N|
−I|N|
0
0
0
0

0
CB
0
0
0
0
0
0

0
−I|N|
Ag + I|N|
−I|N|
−I|N|
0
0
−BCT

0
I|N|
−I|N|
Al + I|N|
I|N|
0
0
BCT

0
0
0
0
0
0
0
−BCT

Proof of Lemma B.1 . Given fixed σ + , σ − , for all (i, j) ∈
/ σ +,
−
(i, j) ∈
/ σ , we have
V̇2 (w) ≤ k(H(w) − w∗ )T · Γ−1
2 (H(w) − (w − F(w)))
T

− (H(w) − w) Γ2 · Q · Γ2 (H(w) − w)
T

− (H(w) − (w − F(w)))

(B.1a)
(B.1b)

(Γ2 − kΓ−1
2 )(w − H(w))
(B.1c)

− k(w − w∗ )T · Γ−1
2 F(w)

(B.1d)

where Q is a semi-definite positive matrix and Γ2 Q = ∇w F(w),
which given in (B.2). Here, the subscript of I means its
dimension, and |A| means the cardinality of set A.
Given fixed σ + and σ − , F(w) is continuous differentiable.
In this case, (B.1a) and (B.1c) are nonpositive due to as
discussed in the proof of Theorem 4 in [1]. (B.1b) is also
nonpositive as Q is semi-definite positive (see Eq. (B.2)).
Next, we prove that (B.1d) is nonpositive. Similar to the
per-node case, substituting µ(t) ≡ ω(t) into the Lagrangian
L2 (x, ρ) in (11) we obtain a function L̂2 (w) defined as follows.


L̂2 (w) := L2 θ̃ , ω, Pg , Pl , λ , µ, η + , η − , φ̃
µ=ω

1 g T g g
l T l l
T
=
(P ) A P + (P ) A P − ω Dω + zT z
2 



+ λ T Pg − Pl − p −CBφ̃ + ω T Pg − Pl − p −CBθ̃


+ (η − )T θ − φ̃ + (η + )T φ̃ − θ
In addition, denote w1 = (θ̃ , Pg , Pl , φ̃ ), w2 = (λ , ω, η + , η − ).
Then L̂2 is convex in w1 and concave in w2 .
In F(w), [φ̃ −θ ]+
and [θ − φ̃ ]+
have unknown dimensions
η+
η−
+
−
(up to σ and σ ). Fortunately, we have
+
+∗ T
(η + − η +∗ )T [φ̃ − θ ]+
η + ≤ (η − η ) (φ̃ − θ )

= (η + − η +∗ )T ∇η + L̂2
where the inequality holds since ηi+j = 0 ≤ ηi+∗
j and φ̃i j −θ i j <
0 for (i, j) ∈ σ + , i.e., (ηi+j − ηi+∗
θ
)
·
(
φ̃
−
ij
i j ) ≥ 0. Similarly,
j
−

(η − η

−∗ T

) [θ

− φ̃ ]+
η−

−

≤ (η − η

−∗ T

−

−∗ T

= (η − η

) (θ − φ̃ )

) ∇η − L̂2 .

Consequently, it can be verified that
∗ T
(w − w∗ )T Γ−1
2 F(w) ≤ (w − w ) ∇w L̂2 (w) 
∇w1 L̂2
= (w − w∗ )T
(w1 , w2 )
−∇w2 L̂2

0
0
0
0
0
0
0
I|E|−|σ + |

0
0
0
0
0
−CB
0
CB
0
CB
0
−I|E|−|σ + |
0
I|E|−|σ − |
−I|E|−|σ − | BCT CB














(B.2)




∇θ̃ L̂2
∇λ L̂2
 ∇ω L̂2 
∇Pg L̂2 



where, ∇w1 L̂2 := 
 ∇Pl L̂2  and ∇w2 L̂2 := ∇η + L̂2 .
∇φ̃ L̂2
∇η − L̂2
Then we have


∗ T
−1
− k(w − w∗ )T Γ−1
2 F(w) ≤ −k(w − w ) · Γ2 F(w)

= − k(w1 − w∗1 )T ∇w1 L̂2 (w1 , w2 ) + k(w2 − w∗2 )T ∇w2 L̂2 (w1 , w2 )

≤ k L̂2 (w∗1 , w2 ) − L̂2 (w1 , w2 ) + L̂2 (w1 , w2 ) − L̂2 (w1 , w∗2 )


= k L̂2 (w∗1 , w2 ) − L̂2 (w∗1 , w∗2 ) + L̂2 (w∗1 , w∗2 ) − L̂2 (w1 , w∗2 )
|
{z
} |
{z
}
≤0

≤0

≤0

(B.3)

where the first inequality holds because L̂2 is convex in w1
and concave in w2 and the second inequality follows because
(w∗1 , w∗2 ) is a saddle point. Therefore (B.1d) is nonpositive,
proving the first assertion.
To prove the second assertion, we further investigate the
situation that σ + or σ − changes. We only consider the set σ +
since it is the same to σ − . We have the following observations:
• The set σ + is reduced, which only happens when φ̃i j −θ i j
goes through zero, from negative to positive. Hence an
extra term will be added to V2 . As this term is initially
zero, there is no discontinuity of V2 in this case.
+
• The set σ + is enlarged when ηi j goes to zero from
positive while φ̃i j < θ i j . Here V2 will lose a positive term
(γiηj )2 (φ̃i j − θ i j )2 /2, causing discontinuity.
In the context, we conclude that V2 is always nonincreasing
along the trajectory even when σ + or σ − changes and
discontinuity occurs.
To prove that the trajectory w(t) is bounded note that [12,
Theorem 3.1] proves that V̂2 (w) := − (H(w) − w)T F(w) −
1
2
2 ||H(w) − w||2 satisfies V̂2 (w) ≥ 0 over S. Therefore, we have
1
∗
k(w(t) − w∗ )T Γ−2
2 (w(t) − w ) ≤ V2 (w(t)) ≤ V2 (w(0))
2
indicating the trajectory w(t) is bounded.
Lemma B.2. Suppose A1, A4 and A5 hold. Then
1) The trajectory w(t) converges to the largest weakly invariant subset W2∗ contained in W2 := {w ∈ S|V̇2 (w) = 0}.
2) Every point w∗ ∈ W2∗ is an equilibrium point of (16).
Proof of Lemma B.2 . Given an initial point w(0) there is a
compact set Ω0 := Ω(w(0)) ⊂ S such that w(t) ∈ Ω0 for t ≥ 0
and V̇2 (w) ≤ 0 in Ω0 .

11

Invoking the proof of Lemma B.1, V2 is radially unbounded
and positively definite except at equilibrium. As V2 and V̇2 are
nonpathological, we conclude that any trajectory w(t) starting
from Ω0 converges to the largest weakly invariant subset W2∗
contained in W2 = { w ∈ Ω0 | V̇2 (w) = 0 } [14, Proposition 3],
proving the first assertion.
For the second assertion, We fix w(0) ∈ W2∗ and then prove
that w(0) must be an equilibrium point.
From (B.1b), direct computing yields
V̇2 (w(t)) ≤ −ω̇ T Dω̇ − (Ṗg )T Ag Ṗg
− (Ṗl )T Al Ṗl − (CBφ̃˙ − Ṗg + Ṗl )T (CBφ̃˙ − Ṗg + Ṗl )
≤0

(B.4)

Since Ag , Al and D are positively definite diagonal matrices,
V̇2 (w) = 0 holds only when Ṗg = Ṗl = ω̇ = 0. Therefore, for
any w(0) ∈ W2∗ , the trajectory w(t) satisfies
Ṗg (t) = Ṗl (t) = ω̇(t) = 0,
Pg (t),

t ≥0

(B.5)

Pl (t)

Hence
and ω(t) are all constants due to the
boundedness property guaranteed by Lemma B.1.
On the other hand, for V̇2 (w) = 0, both terms in (B.3) have
to be zero, implying that
L̂2 (w1 (t), w∗2 ) = L̂2 (w∗1 , w∗2 )
must hold in W2 . Differentiating with respect to t gives

T
∂
∗
L̂2 (w1 (t), w2 ) · ẇ1 (t) = 0 = −φ̃˙ T (Γφ̃ )−1 φ̃˙ (B.6)
∂ w1
The second equality holds due to Eq. (B.5) and (13h). Then we
can conclude φ̃˙ = 0 immediately, implying φ̃ is also constant
in W ∗ due to its boundedness.
Invoking the close-loop dynamics (13), θ̃˙ (t), η̇ + (t), η̇ − (t)
and λ̇ (t) must be constants in W2∗ as Pg (t), Pl (t), ω(t) and
φ̃ (t) are all constants. Then we conclude that θ̃˙ (t) = η̇ + (t) =
η̇ − (t) = λ̇ (t) = 0 holds for all t ≥ 0 due to the boundedness
property of w(t) (Lemma B.1). This implies that any w(0) ∈
W2∗ must be an equilibrium point, completing the proof.
Proof of Theorem 4 . Fix any initial state w(0) and consider
the trajectory (w(t),t ≥ 0) of the closed-loop system closeloop dynamics (13). As mentioned in the proof of Lemma B.2,
w(t) stays entirely in a compact set Ω0 . Hence there exists an
infinite sequence of time instants tk such that w(tk ) → ŵ∗ as
tk → ∞, for some ŵ∗ ∈ W2∗ . Lemma B.2 guarantees that ŵ∗
is an equilibrium point of the closed-loop system (13), and
hence ŵ∗ = H(ŵ∗ ). Thus, using this specific equlibrium point
ŵ∗ in the definition of V2 , we have
V2∗ = lim V2 (w(t)) = lim V2 (w(tk ))
t→∞

tk →∞

=


lim V2 w(tk ) = V2 (ŵ∗ ) = 0

w(tk )→ŵ∗

Here, the first equality uses the fact that V2 (t) is nonincreasing
in t; the second equality uses the fact that tk is the infinite
sequence of t; the third equality uses the fact that w(t) is
absolutely continuous in t; the fourth equality is due to the
upper semi-continuity of V2 (w), and the last equality holds as
ŵ∗ is an equilibrium point of V2 .

∗
The quadratic term (w− ŵ∗ )T Γ−2
2 (w− ŵ ) in V2 then implies
∗
that w(t) → ŵ as t → ∞, which completes the proof.

