1

A unifying energy-based approach to stability of
power grids with market dynamics

arXiv:1604.05200v2 [math.OC] 27 Jul 2016

Tjerk Stegink, Claudio De Persis, Member, IEEE, and Arjan van der Schaft, Fellow, IEEE

Abstract—In this paper a unifying energy-based approach
is provided to the modeling and stability analysis of power
systems coupled with market dynamics. We consider a standard
model of the power network with a third-order model for the
synchronous generators involving voltage dynamics. By applying
the primal-dual gradient method to a social welfare optimization,
a distributed dynamic pricing algorithm is obtained, which can
be naturally formulated in port-Hamiltonian form. By interconnection with the physical model a closed-loop port-Hamiltonian
system is obtained, whose properties are exploited to prove
asymptotic stability to the set of optimal points. This result is
extended to the case that also general nodal power constraints are
included into the social welfare problem. Additionally, the case of
line congestion and power transmission costs in acyclic networks
is covered. Finally, a dynamic pricing algorithm is proposed that
does not require knowledge about the power supply and demand.
Index Terms—port-Hamiltonian, frequency regulation, optimal
power dispatch, dynamic pricing, social welfare, distributed
control, convex optimization.

I. I NTRODUCTION

P

ROVISIONING energy has become increasingly complicated due to several reasons, including the increased share
of renewables. As a result, the generators operate more often
near their capacity limits and transmission line congestion
occurs more frequently.
One effective approach to alleviate some of these challenges
is to use real-time dynamic pricing as a control method
[1]. This feedback mechanism can be used to encourage the
consumers to change their usage when in some parts of the
grid (control areas) it is difficult for the generators and the
network to match the demand.
Real-time dynamic pricing also allows producers and consumers to fairly share utilities and costs associated with the
generation and consumption of energy among the different
control areas. The challenge of achieving this in an optimal
manner while the grid operates within its capacity limits, is
called the social welfare problem [2], [3].
Many of the existing dynamic pricing algorithms focus on
the economic part of optimal supply-demand matching [2],
[4]. However, if market mechanisms are used to determine
the optimal power dispatch (with near real-time updates of
This work is supported by the NWO (Netherlands Organisation for Scientific Research) programme Uncertainty Reduction in Smart Energy Systems
(URSES).
T.W. Stegink and C. De Persis are with the Engineering and Technology
institute Groningen (ENTEG), University of Groningen, 9747 AG Groningen,
the Netherlands. {t.w.stegink, c.de.persis}@rug.nl
A.J. van der Schaft is with the Johann Bernoulli Institute for Mathematics
and Computer Science, University of Groningen, Nijenborgh 9, 9747 AG
Groningen, the Netherlands. a.j.van.der.schaft@rug.nl

the dispatch commands) dynamic coupling occurs between
the market update process and the physical response of the
physical power network dynamics [5].
Consequently, under the assumption of market-based dispatch, it is essential to consider the stability of the coupled
system incorporating both market operation and electromechanical power system dynamics simultaneously.
While on this subject a vast literature is already available,
the aim of this paper is to provide a rigorous and unifying
passivity-based stability analysis. We focus on a more accurate
and higher order model for the physical power network than
conventionally used in the literature. In particular, we use a
third-order model for the synchronous generators including
voltage dynamics. As a result, market dynamics, frequency dynamics and voltage dynamics are considered simultaneously.
Finally, we propose variations of the basic controller design
that, among other things, allow the incorporation of capacity
constraints on the generation and demand of power and on
the transmission lines, and enhance the transient dynamics of
the closed-loop system. The approach taken in this paper is to
model both the dynamic pricing controller as well as the physical network in a port-Hamiltonian way, emphasizing energy
storage and power flow. This provides a unified framework for
the modeling, analysis and control of power networks with
market dynamics, with possible extensions to more refined
models of the physical power network, including for example
turbine dynamics.
A. Literature review
The coupling between a high-order dynamic power network
and market dynamics has been studied before in [5]. Here
a fourth-order model of the synchronous generator is used
in conjunction with turbine and exciter dynamics, which is
coupled to a simple model describing the market dynamics.
The results established in [5] are based on an eigenvalue
analysis of the linearized system.
It is shown in [6] that the third-order model (often called the
flux-decay model) for describing the power network, as used
in the present paper, admits a useful passivity property that
allows for a rigorous stability analysis of the interconnection
with optimal power dispatch controllers, even in the presence
of time-varying demand.
A common way to solve a general optimization problem
like the social welfare problem is by applying the primaldual gradient method [7], [8], [9]. Also in power grids this
is a commonly used approach to design optimal distributed
controllers, see e.g. [10], [11], [12], [13], [14], [15]. The

2

problem formulations vary in these papers, with the focus
being on either the generation side [10], [12], the load side
[11], [16], [17] or both [13], [14], [18], [15]. We will elaborate
on these references in the following two paragraphs.
A vast literature focuses on linear power system models
coupled with gradient-method-based controllers [10], [12],
[11], [17], [13], [19]. In these references the property that
the linear power system dynamics can be formulated as a
gradient method applied to a certain optimization problem is
exploited. This is commonly referred to as reverse-engineering
of the power system dynamics [13], [10], [12]. However, this
approach falls short in dealing with models involving nonlinear
power flows.
Nevertheless, [14], [18], [15], [16] show the possibility to
achieve optimal power dispatch in power networks with nonlinear power flows using gradient-method-based controllers.
On the other hand, the controllers proposed in [14], [18], [15]
have restrictions in assigning the controller parameters and in
addition require that the topology of the physical network is
a tree.
B. Main contributions
The contribution of this paper is to propose a novel energybased approach to the problem that differs substantially from
the aforementioned works. We proceed along the lines of
[20], [21], where a port-Hamiltonian approach to the design
of gradient-method-based controllers in power networks is
proposed. In those papers it is shown that both the power
network as well as the controller designs admit a portHamiltonian representation which are then interconnected to
obtain a closed-loop port-Hamiltonian system.
After showing that the third-order dynamical model describing the power network admits a port-Hamiltonian representation, we provide a systematic method to design gradientmethod-based controllers that is able to balance power supply
and demand while maximizing the social welfare at steady
state. This design is carried out first by establishing the optimality conditions associated with the social welfare problem.
Then the continuous-time gradient method is applied to obtain
the port-Hamiltonian form of the dynamic pricing controller.
Then, following [20], [21], the market dynamics is coupled to
the physical power network in a power-preserving manner so
that all the trajectories of the closed-loop system converge
to the desired synchronous solution and to optimal power
dispatch.
Although the proposed controllers share similarities with
others presented in the literature, the way in which they are
interconnected to the physical network, which is based on
passivity, is to the best of our knowledge new. Moreover, they
show several advantages.
1) Physical model: Since our approach is based on passivity and does not require to reverse-engineer the power system
dynamics as a primal-dual gradient dynamics, it allows to deal
with more complex nonlinear models of the power network.
More specifically, the physical model for describing the power
network in this paper admits nonlinear power flows and timevarying voltages, and is more accurate and reliable than the
classical second-order model [22], [23], [24].

In addition, most of the results that are established in the
present paper are valid for the case of nonlinear power flows
and cyclic networks, in contrast to e.g. [13], [10], [12], [19],
[19], where the power flows are linearized and e.g. [15], [18],
[14] where the physical network topology is a tree. Moreover,
in the aforementioned references the voltages are assumed
to be constant. While the third-order model for the power
network as considered in this paper has been studied before
using passivity based techniques [6], the combination with
gradient method based controllers is novel. In addition, the
stability analysis does not rely on linearization and is based on
energy functions which allow us to establish rigorous stability
results.
2) State transformation: In [10], [12] it is shown how a
state transformation of the closed-loop system can be used to
eliminate the information about the demand from the controller
dynamics, which improves implementation of the resulting
controller. We pursue this idea and show that the same kind
of state transformation can also be used for more complex
physical models as considered in this paper. This avoids the
requirement of knowing the demand to determine the market
price.
3) Controller parameters: In the present paper we show
that both the physical power network as well as the dynamic
pricing controllers admit a port-Hamiltonian representation,
and in particular are passive systems. As a result, the interconnection between the controller and the nonlinear power
system is power-preserving, implying passivity of the closedloop system as well. Consequently, we do not have to impose
any condition on controller design parameters for guaranteeing
asymptotic stability, contrary to [14], [18], [15].
4) Port-Hamiltonian framework: Because of the use of the
port-Hamiltonian framework, the proposed controller designs
have the potential to deal with more complex models for
the power network compared to the model described in this
paper. As long as the more complex model remains portHamiltonian, the results continue to be valid. This may lead
to inclusion of, for example, turbine dynamics or automatic
voltage regulators in the analysis, although this is beyond the
scope of the present paper. Furthermore, higher order models
for the synchronous generator could be considered.
In addition, we propose various extensions to the basic controller design that have not been investigated in the aforementioned references.
5) Transmission costs: In addition to nodal power constraints and line congestion, we also consider the possibility
of including power transmission costs into the social welfare
problem. Including such costs may in particular be useful in
reducing energy losses or the risk of a breakdown of certain
transmission lines.
6) Non-strict convex objective functions: By relaxing the
conditions on the objective function, we show that also nonstrict convex/concave cost/utility functions can be considered
respectively. In addition, the proposed technique allows to add
damping in the gradient method based controller which may
improve the convergence rate of the closed-loop system.

3

7) Barrier functions: We highlight the possibility to use
barrier functions to enforce the trajectories to stay withing the
feasible region, which allows operation within the capacity
constraints for all time, even during transients. This permits a
more realistic application of the proposed controller design.
The remainder of this paper is organized as follows. In
Section II the preliminaries are stated. Then the basic dynamic
pricing algorithm is discussed in Section III and convergence
of the closed-loop system is proven. Variations of the basic
controller design are discussed in Section IV where in Section
IV-A nodal power congestion is included into the social
welfare problem, and in Section IV-B the case line congestion
for the acyclic power networks is discussed. A dynamic pricing
algorithm is proposed in Section IV-C which does not require
knowledge about the power supply and demand. In Section
IV-D the possibility to relax the convexity assumption and
to improve the transient dynamics of the basic controller is
discussed. Finally, the conclusions and suggestions for future
research are discussed in Section V.

II. P RELIMINARIES
A. Notation
Given a symmetric matrix A ∈ Rn×n , we write A > 0 (A ≥
0) to indicate that A is a positive (semi-)definite matrix. The
set of positive real numbers is denoted by R>0 and likewise
the set of vectors in Rn whose elements are positive by Rn>0 .
We write A∨B to indicate that either condition A or condition
B holds, e.g., a > 0 ∨ b > 0 means that either a > 0 or b > 0
holds. Given v ∈ Rn then v  0 denotes the element-wise
inequality. The notation 1 ∈ Rn is used for the vector whose
elements are equal to 1. Given a twice-differentiable function
f : Rn → Rn then the Hessian of f evaluated at x is denoted
by ∇2 f (x). Given a vector η ∈ Rm , we denote by sin η ∈ Rm
the element-wise sine function. Given a differentiable function
f (x1 , x2 ), x1 ∈ Rn1 , x2 ∈ Rn2 then ∇f (x1 , x2 ) denotes the
gradient of f with respect to x1 , x2 evaluated at x1 , x2 and
likewise ∇x1 f (x1 , x2 ) denotes the gradient of f with respect
to x1 . Given a solution x of ẋ = f (x), where f : Rn → Rn
is a Lebesgue measurable function and locally bounded, the
omega-limit set Ω(x) is defined as [25]
n
Ω(x) := x̄ ∈ Rn | ∃{tk }∞
k=1 ⊂ [0, ∞)
o
with lim tk = ∞ and lim x(tk ) = x̄ .
k→∞

are arbitrary labeled with a ‘+’ and a ‘-’, so that the incidence
matrix D of the resulting directed graph is given by

Dik

Each bus represents a control area and is assumed to have a
controllable power supply and demand. The dynamics at each
bus is assumed to be given by [22], [6]
δ̇i = ωi
Mi ω̇i = −

0
0
Bij Eqi
Eqj
sin δij − Ai ωi + Pgi − Pdi

0
0
= Ef i − (1 − (Xdi − Xdi
)Bii )Eqi
X
0
0
− (Xdi − Xdi
)
Bij Eqj
cos δij ,

(1)

j∈Ni

which is commonly referred to as the flux-decay model. Here
we use a similar notation as used in established literature on
power systems [22], [23], [26], [24]. See Table I for a list of
symbols used in the model (1) and throughout the paper.
δi
ωib
ωn
ωi
0
Eqi
Ef i
Pdi
Pgi
Mi
Ni
Ai
Bij
Bii
Xdi
0
Xdi
0
Tdi

Voltage angle
Frequency
Nominal frequency
Frequency deviation ωi := ωib − ω n
q-axis transient internal voltage
Excitation voltage
Power demand
Power generation
Moment of inertia
Set of buses connected to bus i
Asynchronous damping constant
Negative of the susceptance of transmission line (i, j)
Self-susceptance
d-axis synchronous reactance of generator i
d-axis transient reactance of generator i
d-axis open-circuit transient time constant
Table I
PARAMETERS AND STATE VARIABLES OF MODEL (1).

Assumption 1. By using the power network model (1) the
following assumptions are made, which are standard in a broad
range of literature on power network dynamics [22].
•

•

Consider a power grid consisting of n buses. The network is
represented by a connected and undirected graph G = (V, E),
where the nodes, V = {1, ..., n}, is the set of buses and the
edges, E = {1, . . . , m} ⊂ V × V, is the set of transmission
lines connecting the buses. The k-th edge connecting nodes i
and j is denoted as k = (i, j) = (j, i). The ends of edge k

X
j∈Ni

0
0
Tdi
Ėqi

k→∞

B. Power network model



+1 if i is the positive end of edge k
= −1 if i is the negative end of edge k


0
otherwise.

•

Lines are purely inductive, i.e., the conductance is zero.
This assumption is generally valid for the case of high
voltage lines connecting different control areas.
The grid is operating around the synchronous frequency
which implies ωib ≈ ω n for each i ∈ V.
In addition, we assume for simplicity that the excitation
voltage Ef i is constant for all i ∈ V.

Define the voltage angle differences between the buses by
η = DT δ. Further define the angular momenta by p := M ω,
where ω = ω b −1ω n are the (aggregated) frequency deviations
and M = diagi∈V {Mi } are the moments of inertia. Let
0
0
0
0
Γ(Eq0 ) = diagk∈E {γk } and γk = Bij Eqi
Eqj
= Bji Eqi
Eqj

4

where k corresponds to the edge between node i and j. Then
we can write (1) more compactly as [6]
η̇ = DT ω
M ω̇ = −DΓ(Eq0 ) sin η − Aω + Pg − Pd
Td0 Ėq0

=

−F (η)Eq0

(2)

+ Ef

where A = diagi∈V {Ai }, Pg = coli∈V {Pgi }, Pd =
0
0
coli∈V {Pdi }, Td0 = diagi∈V {Tdi
}, Eq0 = coli∈V {Eqi
}, Ef =
coli∈V {Ef i }. For a given η, the components of the matrix
F (η) ∈ Rn×n are defined as
1
0 + Bii ,
Xdi − Xdi
Fij (η) = −Bij cos ηk = Fji (η),
Fii (η) =

i∈V

(3)

k = (i, j) ∈ E

and Fij (η) = 0 otherwise.P
Since for realistic power networks
0
Xdi > Xdi
, and Bii =
j∈Ni Bij > 0 for all i ∈ V, it
follows that F (η) > 0 for all η [22], [23].
Considering the physical energy1 stored in the generator and
the transmission lines respectively, we define the Hamiltonian
as
!
0
2
(E
−
E
)
1X
f
i
qi
Hp =
Mi−1 p2i +
0
2
Xdi − Xdi
i∈V
(4)

1X
0 2
0 2
0
0
+
Bij (Eqi
) + (Eqj
) − 2Eqi
Eqj
cos ηk
2

The assumption η̄ ∈ (−π/2, π/2)m is standard in studies on
power grid stability and is also referred to as a security constraint [6]. In addition, the Hessian condition guarantees the
existence of a local storage function around the equilibrium.
The following result, which establishes decentralized conditions for checking the positive definiteness of the Hessian,
was proven in [28]:
0
Proposition 1. Let Ēqi
∈ Rn>0 and η̄ ∈ (−π/2, π/2)m . If for
all i ∈ V we have
0
X
sin2 η̄k
Eqj
1
+
B
+
B
ii
ij
0
0 cos η̄
Xdi − Xdi
Eqi
k
k=(i,j)∈E


X
Ēqi
>
Bij cos η̄k 1 +
tan2 η̄k > 0
Ēqj
k=(i,j)∈E

then ∇2 Hp (x̄p ) > 0.
It can be verified that the condition stated in Proposition 1 is
satisfied if the following holds [28]:
• the generator reactances are small compared to the transmission line reactances
• the voltage (angle) differences are small.
Remarkably, these conditions hold for a typical operation point
in power transmission networks.

k=(i,j)∈E

where ηk = δi − δj . The first term of the Hamiltonian Hp
represents the shifted kinetic energy stored in the rotors of the
generators and the second term corresponds to the magnetic
energy stored in the generator circuits. Finally, the last term of
Hp corresponds to the magnetic energy stored in the inductive
transmission lines.
By (4), the system (2) can be written in port-Hamiltonian
form [27] as




0 0
0
DT
0
 ∇Hp + I −I  up
0
ẋp = −D −A
0 0
0
0
−Rq
(5)


 
0 I 0
ω
yp =
∇Hp =
0 −I 0
−ω
where xp = col(η, p, Eq0 ), up = col(Pg , Pd ) and
Rq = (Td0 )−1 (Xd − Xd0 ) > 0,
0
Td0 = diagi∈V {Tdi
} > 0,
0
Xd − Xd0 = diagi∈V {Xdi − Xdi
} > 0.

For a study on the stability and equilibria of the flux-decay
model (5), based on the Hamiltonian function (4), we refer to
[6]. The stability results established in [6] rely on the following
assumption.
Assumption 2. Given a constant input up = ūp . There exists
an equilibrium (η̄, p̄, Ēq0 ) of (5) that satisfies η̄ ∈ im DT , η̄ ∈
(−π/2, π/2)m and ∇2 H(η̄, p̄, Ēq0 ) > 0.
aesthetic reasons we define the Hamiltonian Hp as ω n times the
physical energy as the factor 1/ω n appears in each of the energy functions.
As a result, Hp has the dimension of power instead of energy.
1 For

C. Social welfare problem
We define the social welfare by S(Pg , Pd ) := U (Pd ) −
C(Pg ), which consists of a utility function U (Pd ) of the power
consumption Pd and the cost C(Pg ) associated to the power
production Pg . We assume that C(Pg ), U (Pd ) are strictly
convex and strictly concave functions respectively.
Remark 1. It is also possible to include mutual costs and
utilities among the different control areas, provided that the
convexity/concavity assumption is satisfied.
The objective is to maximize the social welfare while
achieving zero frequency deviation. By analyzing the equilibria of (1), it follows that a necessary condition for zero
frequency deviation is 1T Pd = 1T Pg [6], i.e., the total supply
must match the total demand. It can be noted that (Pg , Pd ) is
a solution to the latter equation if and only if there exists a
v ∈ Rmc satisfying Dc v −Pg +Pd = 0 where Dc ∈ Rn×mc is
the incidence matrix of some connected communication graph
with mc edges and n nodes. Because of the latter equivalence,
we consider the following convex minimization problem:
min

Pg ,Pd ,v

− S(Pg , Pd ) = C(Pg ) − U (Pd )

s.t. Dc v − Pg + Pd = 0.

(6a)
(6b)

Remark 2. Although this paper focuses on optimal active
power sharing, we stress that it is also possible to consider
(optimal) reactive power sharing simultaneously, see e.g. [28]
for more details.
The Lagrangian corresponding to (6) is given by
L = C(Pg ) − U (Pd ) + λT (Dc v − Pg + Pd )

(7)

5

with Lagrange multipliers λ ∈ Rn . The resulting first-order
optimality conditions are given by the Karush–Kuhn–Tucker
(KKT) conditions
∇C(P̄g ) − λ̄ = 0,
−∇U (P̄d ) + λ̄ = 0,
DcT λ̄ = 0,

(8)

Dc v̄ − P̄g + P̄d = 0.
Since the minimization problem is convex, strong duality holds
and it follows that (P̄g , P̄d , v̄) is an optimal solution to (6) if
and only if there exists an λ̄ ∈ Rn that satisfies (8) [29].
III. BASIC PRIMAL - DUAL GRADIENT CONTROLLER
In this section we design the basic dynamic pricing algorithm which will be used as the starting point for the
controllers designs discussed in Section IV. Its dynamics is
obtained by applying the primal-dual gradient method [7],
[10], [14] to the minimization problem (6), resulting in
τg Ṗg = −∇C(Pg ) + λ + ugc
τd Ṗd = ∇U (Pd ) − λ +

udc

we interchangeably write the system dynamics in terms of
both xc and zc for ease of notation. In these new variables the
dynamics (9) admits a natural port-Hamiltonian representation
[20], which is given by


0 0 0
I
0 0 0
−I 

ẋc = 
 0 0 0 −DcT  ∇Hc (xc ) + ∇S(zc )
−I I Dc
0


I 0
0 I 

u
+
(10)
0 0 c
0 0


 
I 0 0 0
P
yc =
∇Hc (xc ) = g ,
0 I 0 0
Pd
1
(11)
Hc (xc ) = xTc τc−1 xc .
2

(9a)

Note that the system (10) is indeed a port-Hamiltonian system2
since S is concave and therefore satisfies the incremental
passivity property

(9b)

(z1 − z2 )T (∇S(z1 ) − ∇S(z2 )) ≤ 0, ∀z1 , z2 ∈ R3n+mc .

τv v̇ = −DcT λ

(9c)

τλ λ̇ = Dc v − Pg + Pd .

(9d)

Here we introduce additional inputs uc = col(ugc , udc )
which are to be specified later on, and τc
:=
blockdiag(τg , τd , τv , τλ ) > 0 are controller design parameters.
Recall from Section II-C that there is freedom in choosing a
communication network and the associated incidence matrix.
Depending on the application, one may prefer all-to-all
communication where the underlying graph is complete, or
communication networks where its associated graph is a star,
line or cycle graph. In addition, τc determines the converge
rate of the dynamics (9); a large τc gives a slow convergence
rate whereas a small τc gives a fast convergence rate.
Observe that the dynamics (9) has a clear economic interpretation [1], [2], [5]: each power producer aims at maximizing their own profit, which occurs whenever their individual
marginal cost is equal to the local price λi + ugci . At the same
time, each consumer maximizes its own utility but is penalized
by the local price λi − udci .
The equations (9c), (9d) represent the distributed dynamic
pricing mechanism where the quantity v represents a virtual
power flow along the edges of the communication graph
with incidence matrix Dc . We emphasize virtual, since v
may not correspond to the real physical power flow as the
communication graph may be different than the physical
network topology. Equation (9d) shows that the local price
λi rises if the power demand plus power outflow at node
i ∈ V is greater than the local power supply plus power
inflow of power at node i and vice versa. The inputs ugc , udc are
interpreted as additional penalties or prices that are assigned to
the power producers and consumers respectively. These inputs
can be chosen appropriately to compensate for the frequency
deviation in the physical power network as we will show now.
To this end, define the variables xc = (xg , xd , xv , xλ ) =
(τg Pg , τd Pd , τv v, τλ λ) = τc zc and note that, in the sequel,

The port-Hamiltonian controller (10) is interconnected to the
physical network (5) by taking uc = −yp , up = yc . Define the
extended vectors of variables by
 

η
I 0 0 0 0 0 0
0 M 0 0 0 0 0   ω 
 

0 0 I 0 0 0 0  Eq0 
 

 
x := 
0 0 0 τg 0 0 0   Pg  =: τ z. (12)
0 0 0 0 τd 0 0   Pd 
 

0 0 0 0 0 τv 0   v 
0 0 0 0 0 0 τλ
λ
Then the closed-loop port-Hamiltonian system takes the form


0
DT
0
0
0
0
0
−D −A
0
I −I 0
0 


 0
0
−Rq 0
0
0
0 


−I
0
0
0
0
I 
ẋ = 
 ∇H(x)
 0

 0
I
0
0
0
0
−I


T
 0
0
0
0
0
0 −Dc
0
0
0
−I I Dc
0
+ ∇S(z),
(13)
where H = Hp +Hc is equal to the sum of the energy function
(4) corresponding to the physical model, and the controller
Hamiltonian (11). In the sequel we write (13) more compactly
as
ẋ = (J − R)∇H(x) + ∇S(z),
where R = RT ≥ 0, J = −J T . We define the equilibrium set
of (13), expressed in the variable z, by
Z1 = {z̄ | z̄ is an equilibrium of (13)}.

(14)

Note that each z̄ ∈ Z1 satisfies the optimality conditions
(8) and simultaneously the zero frequency constraints of the
2 Strictly

speaking, (10) is an incremental port-Hamiltonian system [27].

6

physical network (5) given by ω̄ = 0. Hence, Z1 corresponds
to the desired equilibria, and the next theorem states the
convergence to this set of optimal points.
Theorem 1. For every z̄ ∈ Z1 satisfying Assumption 2 there
exists a neighborhood Υ around z̄ where all trajectories z
satisfying (13) with initial conditions in Υ converge to the set
Z1 . In addition, the convergence of each such trajectory is to
a point.
Proof. Let z̄ ∈ Z1 and define the shifted Hamiltonian H̄
around x̄ := τ z̄ as [27], [6]
H̄(x) = H(x) − (x − x̄)T ∇H(x̄) − H(x̄).

(15)

After rewriting, the closed-loop port-Hamiltonian system (13)
is equivalently described by

− (∇Eq0 H̄)T Rq ∇Eq0 H̄ ≤ 0,

A. Including nodal power constraints
The results of Section III can be extended to the case where
nodal constraints on the power production and consumption
are included into the optimization problem (6). To this end,
consider the social welfare problem
min

Pg ,Pd ,v

− S(Pg , Pd ) := C(Pg ) − U (Pd )

s.t. Dc v − Pg + Pd = 0,
g(Pg , Pd )  0

(17a)
(17b)
(17c)

where g : R2n → Rl is a convex function.

ẋ = (J − R)∇H̄(x) + ∇S(z) − ∇S(z̄).
The shifted Hamiltonian H̄ satisfies
H̄˙ = −ω T Aω + (z − z̄)T (∇S(z) − ∇S(z̄))

include, among other things, the possibility to incorporate
nodal power constraints, and line congestion in conjunction
with transmission costs into the social welfare problem.

(16)

where equality holds if and only if ∇Eq0 H̄(x) = ∇Eq0 H(x) =
0, ω = 0, Pg = P̄g , Pd = P̄d since S(z) is strictly concave in
Pg and Pd . Bearing in mind Assumption 2, it is observed
that ∇2 H(x) = ∇2 H̄(x) > 0 for all x in a sufficiently
small open neighborhood around x̄. Hence, as H̄˙ ≤ 0, there
exists a compact sublevel set Υ of H̄ around z̄ contained in
such neighborhood, which is forward invariant. By LaSalle’s
invariance principle, each the solution with initial conditions
in Υ converges to the largest invariant set S contained in
Υ ∩ {z | ∇Eq0 H(x) = 0, ω = 0, Pg = P̄g , Pd = P̄d }.
On such invariant set λ = λ̄ and η, v, Eq0 are constant. Hence,
z converges to S ⊂ Z1 as t → ∞.
Finally, we prove that the convergence of each solution of
(13) initializing in Υ is to a point. This is equivalent to proving
that its omega-limit set Ω(x) is a singleton. Since the solution
x is bounded, Ω(x) 6= ∅ by the Bolzano-Weierstrass theorem
[30]. By contradiction, suppose now that there exist two
distinct point in Ω(x), say x̄1 , x̄2 ∈ Ω(x), x̄1 6= x̄2 . Then there
exists H̄1 (x), H̄2 (x) defined by (15) with respect to x̄1 , x̄2
respectively and scalars c1 , c2 ∈ R>0 such that H̄1−1 (≤ c1 ) :=
{x | H̄1 (x) ≤ c1 }, H̄2−1 (≤ c2 ) := {x | H̄2 (x) ≤ c2 } are
disjoint and compact as the Hessian of H̄1 , H̄2 is positive definite in the neighborhood Υ. Since each trajectory z converges
to Z1 as proven above, it follows that τ −1 x̄1 , τ −1 x̄2 ∈ Z1 .
Together with x̄1 ∈ Ω(x), this implies that there exists a finite
time t1 > 0 such that x(t) ∈ H̄1−1 (≤ c1 ) for all t ≥ t1 as
the set H̄1−1 (≤ c1 ) is invariant by the dissipation inequality
(16). Similarly, there exists a finite time t2 > 0 such that
x(t) ∈ H̄2−1 (≤ c2 ) for all t ≥ t2 . This implies that the
solution x(t) satisfies x(t) ∈ H̄1−1 (≤ c1 ) ∩ H̄1−1 (≤ c1 ) = ∅
for t ≥ max(t1 , t2 ) which is a contradiction. This concludes
the proof.

Remark 3. Note that (17c) captures the convex inequality
constraints considered in the existing literature. For example,
by choosing g as

 

Pg − Pgmax
g1 (Pg , Pd )
g2 (Pg , Pd )  Pgmin − Pg 

 
g(Pg , Pd ) = 
g3 (Pg , Pd ) = Pd − P max  ,
d
g4 (Pg , Pd )
Pdmin − Pd
the resulting inequality constraints (17c) become Pgmin 
Pg  Pgmax , Pdmin  Pd  Pdmax which, among others, are
used in [13], [14], [18].
In the sequel, we assume that (17) satisfies Slater’s condition
[29]. As a result, (P̄g , P̄d , v̄) is an optimal solution to (17)
if and only if there exists λ̄ ∈ Rn , µ̄ ∈ Rl≥0 satisfying the
following KKT optimality conditions:
∂g
(P̄g , P̄d )µ̄ = 0,
∂Pg
∂g
(P̄g , P̄d )µ̄ = 0,
−∇U (P̄d ) + λ̄ +
∂Pd
Dc v̄ − P̄g + P̄d = 0,
DcT λ̄ = 0,
∇C(P̄g ) − λ̄ +

g(P̄g , P̄d )  0,

µ̄  0,

(18)

µ̄T g(P̄g , P̄d ) = 0.

Next, we introduce the following subsystems [20], [9]
(
gi (wi )
if µi > 0
+
ẋµi = (gi (wi ))µi :=
max{0, gi (w)} if µi = 0
(19)
1 T −1
yµi = ∇gi (wi )∇Hµi (xµi ), Hµi (xµi ) = xµi τµi xµi
2
with state xµi := τµi µi ∈ R≥0 , outputs yµi ∈ Rl , inputs
wi ∈ R2n , and i ∈ I := {1, . . . , l}. Here gi (.) is the i’th entry
of the vector-valued function g(.) = coli∈V {gi (.)}. Note that,
for a given i ∈ I and for a constant input w̄i , the equilibrium
set Zµi of (19) is characterized by all (µ̄i , w̄i ) satisfying
gi (w̄i ) ≤ 0,

µ̄i ≥ 0,

µ̄i = 0 ∨ gi (w̄i ) = 0.

(20)

IV. VARIATIONS IN THE CONTROLLER DESIGN

More formally, for i ∈ I the equilibrium set Zµi of (19) is
given by

In this section we propose several variations and extensions
of the controller designed in the previous section. These

Zµi := {(µ̄i , w̄i ) | (µ̄i , w̄i ) satisfies (20)}.

7

Remark 4. In case the inequality constraints of Remark 3 (e.g.
Pg  Pgmax ) are considered, the subsystems (19) take the
decentralized form
(
max
Pgi − Pgi
if µi > 0
max +
ẋµi = (Pgi − Pgi
)µi =
max
max{0, Pgi − Pgi } if µi = 0
1
yµi = ∇Hµi (xµi ), Hµi (xµi ) = xTµi τµ−1
xµi , i ∈ V,
i
2
(21)
and similar expressions can be given for the remaining inequalities Pgmin  Pg , Pdmin  Pd  Pdmax .
The subsystems (19) have the following passivity property
[20].

states, the closed-loop system inherits the same property
provided that an equilibrium of (24) exists.

Proposition 2 ([20]). Let i ∈ I, (µ̄i , w̄i ) ∈ Zµi and define
ȳµi := ∇gi (w̄i )µ̄i . Then (19) is passive with respect to the
shifted external port-variables w̃i := wi − w̄i , ỹµi := yµi − ȳµi .
Additionally, (µi , wi ) → Zµi as t → ∞ for (µi , wi ), wi = w̄i
satisfying (19).

where x̃µ := xµ − x̄µ and H̄ is defined by (23). By Proposition
2 and Remark 5, the time-derivative of H̄e satisfies
X
H̄˙ e ≤ ũT ỹ + w̃T
ỹµi = ũT ỹ − ũT ỹ = 0

Consider again system (13)
ẋ = (J − R)∇H(x) + ∇S(z) + GT u
 

P
0 0 I
y = G∇H(x) = g , G =
Pd
0 0 0

0
I

0 0
0 0



(22)

with the state x defined by (12), where we introduce an
additional input u ∈ R2n and output y ∈ R2n .
Remark 5. Note that for any steady state (x̄, ū) of (22), the
latter system is passive with respect to the shifted external
port-variables ũ := u − ū, ỹ = y − ȳ, ȳ := G∇H(x̄), using
the storage function
H̄(x) := H(x) − (x − x̄)T ∇H(x̄) − H(x̄).

(23)

We interconnect the subsystems (19) to (22) in a powerpreserving way by
X
wi = w = y ∀i ∈ I,
u=−
yµi
i∈I

to obtain the closed-loop system
η̇ = DT ω
M ω̇ =
Td0 Ėq0

=

−DΓ(Eq0 ) sin η −
−F (η)Eq0 + Ef

(24a)
Aω + Pg − Pd

(24c)

∂g
(Pg , Pd )µ − ω
τg Ṗg = −∇C(Pg ) + λ −
∂Pg
∂g
τd Ṗd = ∇U (Pd ) − λ −
(Pg , Pd )µ + ω
∂Pd
τv v̇ = −∇CT (v) − DT λ
τλ λ̇ = Dv − Pg + Pd
τµi µ̇i = (gi (Pg , Pd ))+
µi ,

(24b)

(24d)
(24e)
(24f)
(24g)

i ∈ I.

(24h)

Observe that the equilibrium set Z2 of (24), if expressed in the
co-energy variables, is characterized by all (z̄, µ̄) that satisfy
(18) in addition to ω̄ = 0, −DΓ(Ēq0 ) sin η̄ + P̄g − P̄d =
0, −F (η̄)Ēq0 + Ef = 0, and therefore corresponds to the
desired operation points.
Since both the subsystems (19) and the system (13) admit
an incrementally passivity property with respect to their steady

Theorem 2. For every (z̄, µ̄) ∈ Z2 satisfying Assumption 2
there exists a neighborhood Υ of (z̄, µ̄) where all trajectories
z satisfying (24) with initial conditions in Υ converge to the
set Z2 and the convergence of each such trajectory is to a
point.
Proof. Let (z̄, µ̄) ∈ Z2 and consider the shifted Hamiltonian
H̄e around (x̄, x̄µ ) = (τ z̄, τµ µ̄) defined by
X
1
x̃µ
H̄e (x, xµ ) := H̄(x) +
H̄µi (xµi ) = H̄(x) + x̃Tµ τx−1
µ
2
i∈I

i∈I

where equality holds only if Pg = P̄g , Pd = P̄d , ω =
0, ∇Eq0 H(x) = 0. On the largest invariant set where H̄˙ e = 0 it
follows by the second statement of Proposition 2 that µ = µ̄.
As a result, λ = λ̄ and v, η, Eq0 are constant on this invariant
set. Since the right-hand side of (19) is discontinuous and takes
the same form as in [25], we can apply the invariance principle
for discontinuous Caratheodory systems [25, Proposition 2.1]
to conclude that (z, µ) → Z2 as t → ∞. By following
the same line of arguments as in the proof of Theorem 1,
convergence of each trajectory to a point is proven.
Remark 6. Theorem 2 uses the Caratheodory variant of the
Invariance Principle which requires that the Caratheodory
solution of (24) is unique and that its omega-limit set is
invariant [25]. These requirements are indeed satisfied by
extending Lemmas 4.1-4.4 of [25] to the case where equality
constraints and nonstrict convex/concave (utility) functions are
considered in the optimization problem [25, equation (3)],
noting that these lemmas only require convexity/concavity
instead of their strict versions. In particular, by adding a
quadratic function of the Lagrange multipliers associated with
the equality constraints to the Lyapunov function, it can be
proven that monotonicity of the primal-dual dynamics with
respect to primal-dual optimizers as stated in [25, Lemma 4.1]
holds for this more general case as well, see also [20], [31].
Remark 7. Instead of using the hybrid dynamics (19) for
dealing with the inequality constraints (17c), we can instead introduce the so called barrier functions Bi =
−ν log(−gi (Pg , Pd )) that are added to the objective function
[29]. Simultaneously, the corresponding inequalities in the
social welfare problem (17) are removed to obtain the modified
convex optimization problem
X
min − S(Pg , Pd ) − ν
log(−gi (Pg , Pd ))
Pg ,Pd ,v
(25)
i∈V
s.t. Dc v − Pg + Pd = 0.
Here ν > 0 is called the barrier parameter and is usually
chosen small. By applying the primal-dual gradient method to

8

(25) it can be shown that, if the system is initialized in the
interior of the feasible region, i.e. where (17c) holds, then the
trajectories of the resulting gradient dynamics remain within
the feasible region and the system converges to a suboptimal
value of the social welfare [7], [29], [32]. However, if Slater’s
condition holds, this suboptimal value which depends on ν
converges to the optimal value of the social welfare problem as
ν → 0 [29]. The particular advantage of using barrier functions
is to avoid the use of an hybrid controller and to enforce that
the trajectories remain within the feasible region for all future
time.

the physical model (2), we obtain the following closed-loop
system:
η̇ = DT ω
M ω̇ =
Td0 Ėq0

=

min

Pg ,Pd ,v

− S(Pg , Pd , v) := C(Pg ) + CT (v) − U (Pd ) (26a)

s.t. Dv − Pg + Pd = 0

(26b)

− κ  v  κ,

(26c)

where κ ∈ Rm satisfies the element-wise inequality κ  0.
Note that in this case the communication graph is chosen to
be identical with the topology of the physical network, i.e.,
Dc = D. As a result, the additional constraints (26c) bound
the (virtual) power flow along the transmission lines as |vk | ≤
κk , k ∈ E. The corresponding Lagrangian is given by
L = C(Pg ) + CT (v) − U (Pd ) + λT (Dv − Pg + Pd )
+

µT+ (v

− κ) +

µT− (−κ

with Lagrange multipliers λ ∈ R , µ+ , µ+ ∈
resulting KKT optimality conditions are given by
∇C(P̄g ) − λ̄ = 0,

Rm
≥0 .

(28c)

τd Ṗd = ∇U (Pd ) − λ + ω

(28e)

T

τv v̇ = −∇CT (v) − D λ − µ+ + µ−

(28f)

τλ λ̇ = Dv − Pg + Pd

(28g)

(v − κ)+
µ+
(−κ − v)+
µ− .

(28h)
(28i)

The latter system can partially be put into a port-Hamiltonian
form, since equations (28a)-(28g) can be rewritten as
ẋ = (J − R)∇H(x) + ∇S(z) + N µ

T
0 0 0 0 −I 0
N=
,
0 0 0 0 I 0

(29)

where the variables x, z and the Hamiltonian H are respectively defined by (12) and (13) as before, and µ =
col(µ+ , µ− ).
Since the network topology is a tree (i.e. ker D = {0}),
the equilibrium of (28) satisfies v̄ = Γ(Ēq0 ) sin η̄. Hence, the
controller variable v corresponds to the physical power flow
of the network if the closed-loop system is at steady state.
Consequently, the constraints and costs on v correspond to
constraints and costs of the physical power flow if the system
converges to an equilibrium.
Theorem 3. Let the network topology be acyclic and let (z̄, µ̄)
be an (isolated) equilibrium of (28) satisfying Assumption 2.
Then all trajectories (z, µ) of (28) initialized in a sufficiently
small neighborhood around (z̄, µ̄) converge asymptotically to
(z̄, µ̄)
Proof. Let (z̄, µ̄) be the equilibrium of (28). By defining the
shifted Hamiltonian H̄(x) around x̄ := τ z̄ by

− v)
n

(28b)
(28d)

τ− µ̇− =

The previous section shows how to include nodal power
constraints into the social welfare problem. In case the network
is acyclic, line congestion and power transmission costs can
be incorporated into the optimization problem as well.
To this end, define the (modified) social welfare by U (Pd )−
C(Pg )−CT (v) where the convex function CT (v) corresponds
to the power transmission cost. If security constraints on
the transmission lines are included as well, the optimization
problem (6) modifies to

Aω + Pg − Pd

τg Ṗg = −∇C(Pg ) + λ − ω

τ+ µ̇+ =
B. Including line congestion and transmission costs

(28a)

−DΓ(Eq0 ) sin η −
−F (η)Eq0 + Ef

H̄(x) = H(x) − (x − x̄)T ∇H(x̄) − H(x̄)
The
one can rewrite (29) as
ẋ = (J − R)∇H̄(x) + ∇S(z) − ∇S(z̄) + N µ̃

−∇U (P̄d ) + λ̄ = 0,

(30)

T

∇CT (v̄) + D λ̄ + µ̄+ − µ̄− = 0,
−κ  v̄  κ,
µ̄+ , µ̄−  0,

µ̄T+ (v̄ − κ) = 0,

Dv̄ − P̄g + P̄d = 0,

(27)

µ̄T− (−κ − v̄) = 0.

Suppose that Slater’s condition holds. Then, since the optimization problem (26) is convex, it follows that (P̄g , P̄d , v̄)
is an optimal solution to (26) if and only if there exists
λ̄ ∈ Rn , µ̄ = col(µ̄+ , µ̄− ) ∈ R2m
≥0 satisfying (27) [29].
By applying the gradient method to (26) in a similar
manner as before and connecting the resulting controller with

where µ̃ := µ − µ̄. Consider candidate Lyapunov function
1
1
V (x, µ) = H̄(x) + µ̃+ τµ+ µ̃+ + µ̃T− τµ− µ̃−
2
2
and observe that
T
T
T
µ̃T+ (v − κ)+
µ+ ≤ µ̃+ (v − κ) = µ̃+ (v̄ − κ + ṽ) ≤ µ̃+ ṽ
T
µ̃T− (−κ − v)+
µ− ≤ µ̃− (−κ − v)

= µ̃T− (−κ − v̄ − ṽ) ≤ −µ̃T− ṽ.

9

Bearing in mind (30), the time-derivative of V amounts to

1
1
Ĥ(x̂) = Hp + xTg τg−1 xg + xTd τd−1 xd
2
2
1
1
+ xTv τv−1 xv + (xθ − p)τλ−1 (xθ − p).
2
2
By writing the system of differential equations (31) explicitly
we obtain

V̇ = −ω T Aω − (∇Eq0 H(x))T Rq ∇Eq0 H(x)
+ (z − z̄)T (∇S(z) − ∇S(z̄))
T
+
− ṽ T µ̃+ + ṽ T µ̃− + µ̃T+ (v − κ)+
µ+ + µ̃− (−κ − v)µ−

≤ −ω T Aω + (z − z̄)T (∇S(z) − ∇S(z̄))
− (∇Eq0 H(x))T Rq ∇Eq0 H(x) ≤ 0

η̇ = DT ω

where equality holds only if ∇Eq0 H(x) = 0, ω = 0, Pg =
P̄g , Pd = P̄d . On the largest invariant set S where ∇Eq0 H(x) =
0, ω = 0, Pg = P̄g , Pd = P̄d it follows that, since the graph
contains no cycles λ = λ̄, v = v̄, µ = µ̄ and that η, Eq0 are
constant, which corresponds to an equilibrium. In particular
∇V (x, µ) = 0 for all (z, µ) ∈ S and (z̄, µ̄) ∈ S. Since
by Assumption 2 we have ∇2 V (x̄, µ̄) > 0, it follows that
(z̄, µ̄) is isolated. By the invariance principle for discontinuous Caratheodory systems [25] all trajectories (z, µ) of (28)
initializing in a sufficienly small neighborhood around (z̄, µ̄)
satisfy µ → µ̄, z → z̄ as t → ∞.

Remark 8. It is possible to include nodal power constraints,
line congestion and transmission costs simultaneously. However, as the results in this section are only valid for acyclic
graphs, it should also be assumed for the more general case
that the physical network is a tree.

C. State transformation
Consider again the minimization problem (6). As shown
before, by applying the gradient method to the social welfare
problem, the closed-loop system (13) is obtained.
Note that in the λ-dynamics the demand Pd appears, which
in practice is often uncertain. A possibility to eliminate the demand from the controller dynamics is by a state transformation
[10], [12]. To this end, define the new variables
  

 
η
η
I 0 0 0 0 0 0
p
 p  0 I 0 0 0 0 0
 0 

 0
Eq  0 0 I 0 0 0 0
Eq 
  

 



 Pg  =: τ̂ ẑ,
0
0
0
I
0
0
0
x̂ := 
=
x
=
τ̂
x
 g 

 
 xd   0 0 0 0 I 0 0 
 Pd 
  

 
 xv   0 0 0 0 0 I 0 
v
0 I 0 0 0 0 I
xθ
θ
i.e., xθ := τθ θ = p + xλ . Then
(13) transforms to

0
DT
0
0
−D −A
0
I

 0
0
−Rq 0

˙x̂ =  0
−I
0
0

 0
I
0
0

 0
0
0
0
−D −A
0
0

the port-Hamiltonian system
0
−I
0
0
0
0
0

0
0
0
0
0
0
Dc

with Hamiltonian


DT
−A 

0 

0 
 ∇Ĥ(x̂)
0 

−DcT 
−A

+ ∇S(ẑ)
(31)

M ω̇ = −DΓ(Eq0 ) sin η − Aω + Pg − Pd
Td0 Ėq0 = −F (η)Eq0 + Ef
τg Ṗg = −∇C(Pg ) + τλ−1 (τθ θ − M ω) − ω
τd Ṗd = ∇U (Pd ) −
τv v̇ =

τλ−1 (τθ θ

−DcT τλ−1 (τθ θ

(32)

− M ω) + ω

− M ω)

τθ θ̇ = Dc v − DΓ sin η − Aω.
Define Z4 as the set of all ẑ ∗ := (η̄, ω̄, P̄g , P̄d , v̄, θ̄) that are
an equilibrium of (32). Using the previous established tools
we can prove asymptotic stability to the set of optimal points
Z4 .
Theorem 4. For every ẑ ∗ ∈ Z4 satisfying Assumption 2 there
exists a neighborhood Υ around ẑ ∗ where all trajectories ẑ
satisfying (32) (or equivalently (31)) and initializing in Υ
converge to Z4 . In addition, the convergence of each such
trajectory is to a point.
Proof. We proceed along the same lines as in the proof of
Theorem 1. Since the stability result of Theorem 1 is preserved
after a state transformation, the proof is concluded.
Note that the latter result holds for all τg , τd , τv , τλ , τθ > 0.
The controller appearing in (32) can be simplified by choosing
τλ = τθ = M . As a result, the controller dynamics is described
by
τg Ṗg = −∇C(Pg ) + θ − 2ω

(33a)

τd Ṗd = ∇U (Pd ) − θ + 2ω

(33b)

τv v̇ =

−DcT (θ

M θ̇ = Dc v −

− ω)

DΓ(Eq0 ) sin η

(33c)
− Aω.

(33d)

The main advantage of controller design (33) is that no
information about the power supply and demand is required
in the dynamic pricing algorithm (33c), (33d), where we
observe that the quantity θ − 2ω acts here as the electricity
price for the producers and consumers. Another benefit of the
proposed dynamic pricing algorithm is that, contrary to [16],
no information is required about ω̇.
On the other hand, knowledge about the physical power
flows and the power system parameters M, A is required.
Determining the radius of uncertainty of these parameters
under which asymptotic stability is preserved remains an open
question [10]; see [17] for results in a similar setting where
only the damping term A is assumed to be uncertain.
D. Relaxing the strict convexity assumption
By making a minor modification to the social welfare
problem (6), it is possible to relax the condition that the

10

functions C, U are strictly convex and concave respectively.
To this end, consider the optimization problem
1
C(Pg ) − U (Pd ) + ρ||Dc v − Pg + Pd ||2
Pg ,Pd ,v
2
s.t. Dc v − Pg + Pd = 0,
min

(34a)
(34b)

where ρ > 0, C(Pg ) is convex and U (Pd ) is concave, which
makes the optimization problem (34) convex. Suppose that
there exists a feasible solution to the minimization problem,
then the set of optimal points of (34) is identical with the
set of optimal points of (6) which is characterized by set of
points satisfying the KKT conditions (8). The corresponding
augmented Lagrangian of (34) is given by
Lp = C(Pg ) − U (Pd ) − λT (Dc v + Pg − Pd )
1
+ ρ||Dc v + Pg − Pd ||2 .
2
Consequenctly, the distributed dynamics of the primal-dual
gradient method applied to (34) amounts to
τg Ṗg = −∇C(Pg ) + λ − ρ(Dc v + Pg − Pd )
τd Ṗd = ∇U (Pd ) − λ + ρ(Dc v + Pg − Pd )
τv v̇ = DcT λ − ρDcT (Dc v + Pg − Pd )

(35)

τλ λ̇ = −Dc v − Pg + Pd ,
which can be written in the same port-Hamiltonian form as
(13) where in this case
1
S(Pg , Pd , v) = U (Pd ) − C(Pg ) − ρ||Dc v − Pg + Pd ||2 .
2
(36)
This leads to the following result.
Theorem 5. Consider the system (13) where S is given by
(36) and suppose that C, U are convex and concave functions
respectively. Then for every z̄ ∈ Z1 satisfying Assumption 2,
where Z1 is defined by (14), there exists a neighborhood Υ
around z̄ wherein each trajectory z satisfying (13) converges
to a point in Z1 .
Proof. Let z̄ ∈ Z1 . By the proof of Theorem 1 it follows that
H̄˙ = −ω T Aω + (z − z̄)T (∇S(z) − ∇S(z̄))
− (∇Eq0 H̄)T Rq ∇Eq0 H̄,

S ⊂ Z1 . By continuing along the same lines as the proof
of Theorem 1, convergence of each trajectory to a point is
proven.
Remark 9. Adding the quadratic term in the social welfare
problem as done in (34a) provides an additional advantage.
As this introduces more damping in the resulting gradientmethod-based controller, see (37), it may improve the convergence properties of the closed-loop dynamics [33], [34].
Moreover, the amount of damping injected into the system
depends on parameter ρ, which can be chosen freely.
V. C ONCLUSIONS AND FUTURE RESEARCH
In this paper a unifying and systematic energy-based approach in modeling and stability analysis of power networks
has been established. Convergence of the closed-loop system
to the set of optimal points using gradient-method-based
controllers have been proven using passivity based arguments.
This result is extended to the case where nodal power constraints are included into the problem as well. However,
for line congestion and power transmission cost the power
network is required to be acyclic to prove asymptotic stability
to the set of optimal points.
The results established in this paper lend themselves to
many possible extensions. One possibility is to design an
additional (passive) controller that regulates the voltages to the
desired values or achieves alternative objectives like (optimal)
reactive power sharing. This could for example be realized by
continuing along the lines of [28].
Recent observations, see [35], suggest that the portHamiltonian framework also lends itself to consider higherdimensional models for the synchronous generator than the
third-order model used in this paper, while the same controllers
as designed in the present paper can be used in this case
as well. In addition, current research includes extending the
results of the present paper to network-preserving models
where a distinction is made between generator and load nodes.
One of the remaining open questions is how to deal with
line congestion and power transmission costs in cyclic power
networks with nonlinear power flows. In addition, all of
the results established for the nonlinear power network only
provide local asymptotic stability to the set of optimal points.
Future research includes determining the region of attraction.

where the second term can be written as
P̃dT (∇U (Pd ) − ∇U (P̄d )) − P̃gT (∇C(Pg ) − ∇C(P̄g ))
 T 
 
−I
I
−Dc
P̃g
P̃g
−I
Dc  P̃d  ≤ 0
− ρ P̃d   I
(37)
−DcT DcT −DcT Dc
ṽ
ṽ
where P̃g = Pg − P̄g , P̃d = Pd − P̄d , ṽ = v − v̄. Hence,
we obtain that H̄˙ ≤ 0 where equality holds only if ω =
0, ∇Eq0 H̄(x) = 0 and Dc ṽ + P̃g − P̃d = Dc v + Pg − Pd = 0.
On the largest invariant set S where H̄˙ = 0 we have ω = 0
c

and η, Eq0
optimality
invariance
where all

are constant and (Pg , Pd , v, λ) satisfy the KKT
conditions (8). Therefore S ⊂ Z1 and by LaSalle’s
principle there exists a neighborhood Υ around z̄
trajectories z satisfying (13) converge to the set

R EFERENCES
[1] F. Alvarado, “The stability of power system markets,” IEEE Transactions
on Power Systems, vol. 14, no. 2, pp. 505–511, May 1999.
[2] A. Kiani and A. Annaswamy, “The effect of a smart meter on congestion
and stability in a power market,” in 49th IEEE Conference on Decision
and Control, Atlanta, USA, December, 2010.
[3] ——, “A hierarchical transactive control architecture for renewables
integration in smart grids,” in Decision and Control (CDC), 2012 IEEE
51st Annual Conference on. IEEE, 2012, pp. 4985–4990.
[4] M. Roozbehani, M. Dahleh, and S. Mitter, “On the stability of wholesale
electricity markets under real-time pricing,” in 49th IEEE Conference on
Decision and Control (CDC), 2010, pp. 1911–1918.
[5] F. Alvarado, J. Meng, C. DeMarco, and W. Mota, “Stability analysis
of interconnected power systems coupled with market dynamics,” IEEE
Transactions on Power Systems, vol. 16, no. 4, pp. 695–701, 2001.
[6] S. Trip, M. Bürger, and C. De Persis, “An internal model approach
to (optimal) frequency regulation in power grids with time-varying
voltages,” Automatica, vol. 64, pp. 240–253, 2016.

11

[7] J. Arrow, L. Hurwicz, H. Uzawa, and H. Chenery, Studies in linear and
non-linear programming. Stanford University Press, 1958.
[8] D. Feijer and F. Paganini, “Stability of primal–dual gradient dynamics
and applications to network optimization,” Automatica, vol. 46, pp.
1974–1981, September 2010.
[9] A. Jokić, M. Lazar, and P. P. Van den Bosch, “On constrained steadystate regulation: dynamic KKT controllers,” Automatic Control, IEEE
Transactions on, vol. 54, no. 9, pp. 2250–2254, 2009.
[10] N. Li, L. Chen, C. Zhao, and S. H. Low, “Connecting automatic
generation control and economic dispatch from an optimization view,”
in American Control Conference. IEEE, 2014, pp. 735–740.
[11] E. Mallada and S. Low, “Distributed frequency-preserving optimal load
control,” in IFAC World Congress, 2014.
[12] Y. Seungil and C. Lijun, “Reverse and forward engineering of frequency
control in power networks,” in Proc. of IEEE Conference on Decision
and Control, Los Angeles, CA, USA, 2014.
[13] X. Zhang, N. Li, and A. Papachristodoulou, “Achieving real-time economic dispatch in power networks via a saddle point design approach,”
in Power & Energy Society General Meeting, 2015. IEEE, 2015, pp.
1–5.
[14] X. Zhang and A. Papachristodoulou, “A real-time control framework
for smart power networks with star topology,” in American Control
Conference (ACC), 2013, pp. 5062–5067.
[15] ——, “A real-time control framework for smart power networks: Design
methodology and stability,” Automatica, vol. 58, pp. 43–50, 2015.
[16] C. Zhao, E. Mallada, and S. Low, “Distributed generator and loadside secondary frequency control in power networks,” in 49th Annual
Conference on Information Sciences and Systems (CISS). IEEE, 2015,
pp. 1–6.
[17] E. Mallada, C. Zhao, and S. Low, “Optimal load-side control for
frequency regulation in smart grids,” in Communication, Control, and
Computing (Allerton), 2014 52nd Annual Allerton Conference on.
IEEE, 2014, pp. 731–738.
[18] X. Zhang and A. Papachristodoulou, “Distributed dynamic feedback
control for smart power networks with tree topology,” in American
Control Conference (ACC), 2014, pp. 1156–1161.
[19] C. Zhao, U. Topcu, N. Li, and S. Low, “Design and stability of load-side
primary frequency control in power systems,” Automatic Control, IEEE
Transactions on, vol. 59, no. 5, pp. 1177–1189, 2014.
[20] T. W. Stegink, C. De Persis, and A. J. van der Schaft, “Port-Hamiltonian
formulation of the gradient method applied to smart grids,” IFACPapersOnLine, vol. 48, no. 13, pp. 13–18, 2015.
[21] ——, “A port-Hamiltonian approach to optimal frequency regulation in
power grids,” arXiv preprint arXiv:1509.07318, 2015.
[22] J. Machowski, J. Bialek, and J. Bumby, Power System Dynamics:
Stability and Control, 2nd ed. Ltd: John Wiley & Sons, 2008.
[23] P. Kundur, Power System Stability and Control. Mc-Graw-Hill Engineering, 1993.
[24] P. Sauer and M. Pai, Power system dynamics and stability. PrenticeHall, 1998.
[25] A. Cherukuri, E. Mallada, and J. Cortés, “Asymptotic convergence of
constrained primal–dual dynamics,” Systems & Control Letters, vol. 87,
pp. 10–15, 2016.
[26] P. Anderson and A. Fouad, Power System Control and Stability, 1st ed.
The Iowa State Univsersity Press, 1977.
[27] A. van der Schaft and D. Jeltsema, “Port-Hamiltonian systems theory:
An introductory overview,” Foundations and Trends in Systems and
Control, vol. 1, no. 2-3, pp. 173–378, 2014.
[28] C. De Persis and N. Monshizadeh, “A modular design of incremental
Lyapunov functions for microgrid control with power sharing,” arXiv
preprint arXiv:1510.05811, 2015.
[29] S. Boyd and L. Vandenberghe, Convex Optimization, 1st ed. Cambridge
University Press, 2004.
[30] W. Rudin, Principles of mathematical analysis. McGraw-Hill New
York, 1964, vol. 3.
[31] A. Cherukuri, B. Gharesifard, and J. Cortes, “Saddle-point dynamics:
conditions for asymptotic stability of saddle points,” arXiv preprint
arXiv:1510.02145, 2015.
[32] J. Wang and N. Elia, “A control perspective for centralized and distributed convex optimization,” in Decision and Control and European
Control Conference (CDC-ECC), 2011 50th IEEE Conference on.
IEEE, 2011, pp. 3800–3805.
[33] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method
of multipliers,” Foundations and Trends in Machine Learning, vol. 3,
no. 1, pp. 1–122, 2011.

[34] R. T. Rockafellar, “The multiplier method of Hestenes and Powell
applied to convex programming,” Journal of Optimization Theory and
applications, vol. 12, no. 6, pp. 555–562, 1973.
[35] T. W. Stegink, C. De Persis, and A. J. van der Schaft, “Optimal
power dispatch in networks of high-dimensional models of synchronous
machines,” arXiv preprint arXiv:1603.06688, 2016.

Tjerk Stegink is a Ph.D. candidate at the Engineering and Technology Institute, Faculty of Mathematics and Natural Sciences, University of Groningen,
the Netherlands. He received his B.Sc. (2012) in
Applied Mathematics and M.Sc. (2014, cum laude)
in Systems, Control and Optimization from the
same university. His main research interests are in
modeling and nonlinear distributed control of power
systems.

Claudio De Persis received the Laurea degree (cum
laude) in electrical engineering in 1996 and the
Ph.D. degree in system engineering in 2000, both
from Sapienza University of Rome, Rome, Italy. He
is currently a Professor at the Engineering and Technology Institute, Faculty of Mathematics and Natural
Sciences, University of Groningen, the Netherlands.
He is also affiliated with the Jan Willems Center for Systems and Control. Previously he was
with the Department of Mechanical Automation and
Mechatronics, University of Twente and with the
Department of Computer, Control, and Management Engineering, Sapienza
University of Rome. He was a Research Associate at the Department of
Systems Science and Mathematics, Washington University, St. Louis, MO,
USA, in 2000–2001, and at the Department of Electrical Engineering, Yale
University, New Haven, CT, USA, in 2001–2002. His main research interest
is in control theory, and his recent research focuses on dynamical networks,
cyberphysical systems, smart grids and resilient control. He was an Editor
of the International Journal of Robust and Nonlinear Control (2006–2013),
and is currently an Associate Editor of the IEEE Transactions On Control
Systems Technology (2010–2015), of the IEEE Transactions On Automatic
Control (2012–2015), and of Automatica (2013–present).

Arjan van der Schaft received the undergraduate
(cum laude) and Ph.D. degrees in Mathematics from
the University of Groningen, The Netherlands. In
1982 he joined the Department of Applied Mathematics, University of Twente, where he was appointed as full professor in Mathematical Systems
and Control Theory in 2000. In September 2005 he
returned to his Alma Mater as a full professor in
Mathematics.
Arjan van der Schaft is Fellow of the Institute
of Electrical and Electronics Engineers (IEEE), and
Fellow of the International Federation of Automatic Control (IFAC). He was
Invited Speaker at the International Congress of Mathematicians, Madrid,
2006. He was the 2013 recipient of the 3-yearly awarded Certificate of
Excellent Achievements of the IFAC Technical Committee on Nonlinear
Systems.
He is (co-)author of the following books: System Theoretic Descriptions
of Physical Systems (1984), Variational and Hamiltonian Control Systems
(1987, with P.E. Crouch), Nonlinear Dynamical Control Systems (1990,
with H. Nijmeijer), L2-Gain and Passivity Techniques in Nonlinear Control
(1996, 2000), An Introduction to Hybrid Dynamical Systems (2000, with J.M.
Schumacher), Port-Hamiltonian Systems Theory: An Introductory Overview
(2014, with D. Jeltsema).

