Optimal power dispatch in networks of high-dimensional models of
synchronous machines

arXiv:1603.06688v1 [math.OC] 22 Mar 2016

Tjerk Stegink and Claudio De Persis and Arjan van der Schaft
Abstract— This paper investigates the problem of optimal
frequency regulation of multi-machine power networks where
each synchronous machine is described by a sixth order model.
By analyzing the physical energy stored in the network and
the generators, a port-Hamiltonian representation of the multimachine system is obtained. Moreover, it is shown that the openloop system is passive with respect to its steady states which
implies that passive controllers can be used to control the multimachine network. As a special case, a distributed consensus
based controller is designed that regulates the frequency and
minimizes a global quadratic generation cost in the presence
of a constant unknown demand. In addition, the proposed
controller allows freedom in choosing any desired connected
undirected weighted communication graph.

I. INTRODUCTION
The control of power networks has become increasingly
challenging over the last decades. As renewable energy
sources penetrate the grid, the conventional power plants
have more difficulty in keeping the frequency around the
nominal value, e.g. 50 Hz, leading to an increased chance of
a network failure of even a blackout.
The current developments require that more advanced
models for the power network must be established as the
grid is operating more often near its capacity constraints.
Considering high-order models of, for example, synchronous
machines, that better approximate the reality allows us to establish results on the control and stability of power networks
that are more reliable and accurate.
At the same time, incorporating economic considerations
in the power grid has become more difficult. As the scale of
the grid expands, computing the optimal power production
allocation in a centralized manner as conventionally is done
is computationally expensive, making distributed control far
more desirable compared to centralized control. In addition,
often exact knowledge of the power demand is required for
computing the optimal power dispatch, which is unrealistic
in practical applications. As a result, there is an increased
desire for distributed real-time controllers which are able to
compensate for the uncertainty of the demand.
In this paper, we propose an energy-based approach for the
modeling, analysis and control of the power grid, both for
This work is supported by the NWO (Netherlands Organisation for
Scientific Research) programme Uncertainty Reduction in Smart Energy
Systems (URSES) under the auspices of the project ENBARK.
T.W. Stegink and C. De Persis are with the Engineering and Technology
institute Groningen (ENTEG), University of 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 physical network as well as for the distributed controller
design. Since energy is the main quantity of interest, the portHamiltonian framework is a natural approach to deal with the
problem. Moreover, the port-Hamiltonian framework lends
itself to deal with complex large-scale nonlinear systems like
power networks [5], [12], [13].
The emphasis in the present paper lies on the modeling and
control of (networked) synchronous machines as they play an
important role in the power network since they are the most
flexible and have to compensate for the increased fluctuation
of power supply and demand. However, the full-order model
of the synchronous machine as derived in many power
engineering books like [2], [6], [8] is difficult to analyze,
see e.g. [5] for a port-Hamiltonian approach, especially when
considering multi-machine networks [4], [9]. Moreover, it is
not necessary to consider the full-order model when studying
electromechanical dynamics [8].
On the other hand of the spectrum, many of the recent optimal controllers in power grids that deal with optimal power
dispatch problems rely on the second-order (non)linear swing
equations as the model for the power network [7], [11], [16],
[17], or the third-order model as e.g. in [14]. However, the
swing equations are inaccurate and only valid on a specific
time scale up to the order of a few seconds so that asymptotic
stability results are often invalid for the actual system [2],
[6], [8].
Hence, it is appropriate to make simplifying assumptions
for the full-order model and to focus on multi-machine
models with intermediate complexity which provide a more
accurate description of the network compared to the secondand third-order models [2], [6], [8]. However, for the resulting intermediate-order multi-machine models the stability
analysis is often carried out for the linearized system, see
[1], [6], [8]. Consequently, the stability results are only valid
around a specific operating point.
Our approach is different as the nonlinear nature of the
power network is preserved. More specifically, in this paper
we consider a nonlinear sixth-order reduced model of the
synchronous machine that enables a quite accurate description of the power network while allowing us to perform a
rigorous analysis.
In particular, we show that the port-Hamiltonian framework is very convenient when representing the dynamics
of the multi-machine network and for the stability analysis.
Based on the physical energy stored in the generators and
the transmission lines, a port-Hamiltonian representation of
the multi-machine power network can be derived. More
specifically, while the system dynamics is complex, the in-

terconnection and damping structure of the port-Hamiltonian
system is sparse and, importantly, state-independent.
The latter property implies shifted passivity of the system
[15] which respect to its steady states which allows the
usage of passive controllers that steer the system to a desired
steady state. As a specific case, we design a distributed realtime controller that regulates the frequency and minimizes
the global generation cost without requiring any information about the unknown demand. In addition, the proposed
controller design allows us to choose any desired undirected
weighted communication graph as long as the underlying
topology is connected.
The main contribution of this paper is to combine distributed optimal frequency controllers with a high-order
nonlinear model of the power network, which is much
more accurate compared to the existing literature, to prove
asymptotic stability to the set of optimal points by using
Lyapunov function based techniques.
The rest of the paper is organized as follows. In Section
II the preliminaries are stated and a sixth order model of
a single synchronous machine is given. Next, the multimachine model is derived in Section III. Then the energy
functions of the system are derived in Section IV, which
are used to represent the multi-machine system in portHamiltonian form, see Section V. In Section VI the design of
the distributed controller is given and asymptotic stability to
the set of optimal points is proven. Finally, the conclusions
and possibilities for future research are discussed in Section
VII.
II. PRELIMINARIES
Consider a power grid consisting of n buses. The network
is represented by a connected and undirected graph G =
(V, E), where the set of nodes, V = {1, ..., n}, is the set of
buses and the set of edges, E = {1, ..., m} ⊂ V × V, is the
set of transmission lines connecting the buses. The ends of
edge l ∈ E are arbitrary labeled with a ‘+’ and a ‘-’, so that
the incidence matrix D of the network is given by


+1 if i is the positive end of l
Dil = −1 if i is the negative end of l
(1)


0
otherwise.
Each bus represents a synchronous machine and is assumed to have controllable mechanical power injection and
a constant unknown power load. The dynamics of each
synchronous machine i ∈ V is assumed to be given by [8]
Mi ω̇i = Pmi − Pdi − Vdi Idi − Vqi Iqi

δ̇i
0
0
Tdi Ėqi
0
0
Tqi
Ėdi
00 00
Tdi
Ėqi
00 00
Tqi Ėdi
see also Table I.

0
0
= −Edi
− (Xqi − Xqi
)Iqi

0
00
0
00
= Eqi
− Eqi
+ (Xdi
− Xdi
)Idi

0
00
0
00
= Edi
− Edi
− (Xqi
− Xqi
)Iqi ,

rotor angle w.r.t. synchronous reference frame
frequency deviation
mechanical power injection
power demand
moment of inertia
synchronous reactances
transient reactances
subtransient reactances
exciter emf/voltage
internal bus transient emfs/voltages
internal bus subtransient emfs/voltages
external bus voltages
generator currents
open-loop transient time-scales
open-loop subtransient time-scales

TABLE I
M ODEL PARAMETERS AND VARIABLES .

Assumption 1: When using model (2), we make the following simplifying assumptions [8]:
•
•
•

•

The frequency of each machine is operating around the
synchronous frequency.
The stator winding resistances are zero.
The excitation voltage Ef i is constant for all i ∈ V.
00
The subtransient saliency is negligible, i.e. Xdi
=
00
Xqi , ∀i ∈ V.

The latter assumption is valid for synchronous machines with
damper windings in both the d and q axes, which is the case
for most synchronous machines [8].
It is standard in the power system literature to represent
the equivalent synchronous machine circuits along the dqaxes as in Figure 1, [6], [8]. Here we use the conventional
00
00
00
00
00
00
+jEdi
where E qi :=
phasor notation E i = E qi +E di = Eqi
00
00
00
Eqi
, E di := jEdi
, and the phasors I i , V i are defined likewise
[8], [10]. Remark that internal voltages Eq0 , Ed0 , Eq00 , Ed00 as
depicted in Figure 1 are not necessarily at steady state but are
governed by (2), where it should be noted that, by definition,
the reactances of a round rotor synchronous machine satisfy
0
00
0
00
Xdi > Xdi
> Xdi
> 0, Xqi > Xqi
> Xqi
> 0 for all i ∈ V
[6], [8].
By Assumption 1 the stator winding resistances are negligible so that synchronous machine i can be represented
by a subtransient emf behind a subtransient reactance, see
Figure 2 [6], [8]. As illustrated in this figure, the internal
and external voltages are related to each other by [8]
00

00
E i = V i + jXdi
I i,

i ∈ V.

(3)

III. MULTI-MACHINE MODEL

= ωi

0
0
= Ef i − Eqi
+ (Xdi − Xdi
)Idi

δi
ωi
Pmi
Pdi
Mi
Xqi , Xdi
0 , X0
Xqi
di
00 , X 00
Xdi
qi
Ef i
0 , E0
Eqi
di
00 , E 00
Eqi
di
Vqi , Edi
Iqi , Idi
0 ,T0
Tqi
di
00 , T 00
Tqi
di

(2)

Consider n synchronous machines which are interconnected by RL-transmission lines and assume that the network
is operating at steady state. As the currents and voltages of
each synchronous machine is expressed w.r.t. its local dqreference frame, the network equations are written as [10]
00

I = diag(e−jδi )Y diag(ejδi )E .

(4)

j(Xd − Xd0 )

Td0

Ef

j(Xd0 − Xd00 )
0

Eq

jXd00

00

Td00

Id

Vq

Eq

To simplify the analysis further, we assume that the
network resistances are negligible so that G = 0. By equating
the real and imaginary part of (4) we obtain the following
expressions for the dq-currents entering generator i ∈ V:
X 

00
00
00
Bik (Edk
Idi = Bii Eqi
−
sin δik + Eqk
cos δik ) ,
k∈Ni

j(Xq − Xq0 )

Tq0

j(Xq0 − Xq00 )
0

Ed

jXq00

00

Tq00

Ed

Iqi =
Iq

Vd

Fig. 1: Generator equivalent circuits for both dq-axes [8].
For aesthetic reasons the subscript i is dropped.

00
jXdi

00

Ei

Vi

that D is the incidence matrix of the network defined by (1).

00

Ei

jXT

Vi

Il

X 

00
00
Bik (Eqk
sin δik − Edk
cos δik ) ,

k∈Ni

IV. ENERGY FUNCTIONS

−1 T
Here the admittance matrix1 Y := D(R + jX)P
D satisfies
YP
ik = −Gik − jBik and Yii = Gii + jBii =
k∈Ni Gik +
j k∈Ni Bik where G denotes the conductance and B ∈
Rn×n
denotes the susceptance of the network [10]. In
≤0
addition, Ni denotes the set of neighbors of node i.
Remark 1: As the electrical circuit depicted in Figure 2 is
00
in steady state (3), the reactance Xdi
can also be considered
as part of the network (an additional inductive line) and
is therefore implicitly included into the network admittance
matrix Y, see also Figure 3.

00
jXdi

−

(5)
where δik := δi − δk . By substituting (5) and (3) into (2)
we obtain after some rewriting a sixth-order multi-machine
model given by equation (6), illustrated at the top of the next
page.
Remark 2: Since the transmission lines are purely inductive by assumption, there are no energy losses in the
transmission lines implying
that the following energy conserP
∗
vation law holds: i∈V Pei = 0 where Pei = Re(E i I i ) =
00
00
Edi Idi + Eqi Iqi is the electrical power produced by synchronous machine i.

Ii

Fig. 2: Subtransient emf behind a subtransient reactance.

1 Recall

00
−Bii Edi

00
jXdk

00

V k Ek

Fig. 3: Interconnection of two synchronous machines by a
purely inductive transmission line with reactance XT .

When analyzing the stability of the multi-machine system
one often searches for a suitable Lyapunov function. A
natural starting point is to consider the physical energy as a
candidate Lyapunov function. Moreover, when we have an
expression for the energy, a port-Hamiltonian representation
of the associated multi-machine model (6) can be derived,
see Section V.
Remark 3: It is convenient in the definition of the Hamiltonian to multiply the energy stored in the synchronous
machine and the transmission lines by the synchronous
frequency ωs since a factor ωs−1 appears in each of the energy
functions. As a result, the Hamiltonian has the dimension of
power instead of energy. Nevertheless, we still refer to the
Hamiltonian as the energy function in the sequel.
In the remainder of this section we will first identify the
electrical and mechanical energy stored in each synchronous
machine. Next, we identify the energy stored in the transmission lines.
A. Synchronous Machine
1) Electrical Energy: Note that, at steady state, the energy
(see Remark 3) stored in the first two reactances2 of generator
i as illustrated in Figure 1 is given by
!
0
0
00 2
− Ef i )2
(Eqi
− Eqi
)
1 (Eqi
Hedi =
+
0
0 − X 00
2
Xdi − Xdi
Xdi
di
!
(7)
0 2
0
00 2
1
(Edi )
(Edi − Edi )
Heqi =
.
0 + X 0 − X 00
2 Xqi − Xqi
qi
qi
Remark 4: The energy stored in the third (subtransient)
reactance will be considered as part of the energy stored in
the transmission lines, see also Remark 1 and Section IV-B.
2 In

both the d- and the q-axes.

Mi ∆ω̇i = Pmi − Pdi +

X
k∈Ni

h
i
00 00
00 00
00 00
00 00
Bik (Edi
Edk + Eqi
Eqk ) sin δik + (Edi
Eqk − Eqi
Edk ) cos δik

δ̇i = ∆ωi
0
0
0
0
00
Tdi
Ėqi
= Ef i − Eqi
+ (Xdi − Xdi
)(Bii Eqi
−
0
0
Tqi
Ėdi

=

0
−Edi

+ (Xqi −

0
00
Xqi
)(Bii Edi

=

0
Edi

−

00
Edi

+

0
(Xqi

−

00
00
Xqi
)(Bii Edi

B. Inductive Transmission Lines
Consider an interconnection between two SG’s with a
purely inductive transmission line (with reactance XT ) at
steady state, see Figure 3. When expressed in the local dqreference frame of generator i, we observe from Figure 3
that at steady state one obtains3
00

jXl I l = E i − e−jδik E k ,

(8)

where the total reactance between the internal buses of
00
00
generator i and k is given by Xl := Xdi
+ XT + Xdk
.
Note that at steady state the modified energy of the inductive
transmission line l between nodes i and k is given by
∗
Hl = 21 Xl I l I l , which by (8) can be rewritten as


1
00 00
00
00
Hl = − Bik 2 Edi
Eqk − Edk
Eqi
sin δik
2

00 00
00 00
(9)
−2 Edi Edk + Eqi Eqk cos δik

002
002
002
002
+Edi
+ Edk
+ Eqi
+ Eqk
,
where the line susceptance satisfies Bik = − X1l < 0 [10].
C. Total Energy
The total physical energy of the multi-machine system is
equal to the sum of the individual energy functions:
X
X
Hp =
(Hdei + Hqei + Hmi ) +
Hl .
(10)
i∈V

(6)

k∈Ni

2) Mechanical Energy: The kinetic energy of synchronous machine i is given by
1
1
Hmi = Mi ωi2 = Mi−1 p2i ,
2
2
where pi = Mi ωi is the angular momentum of synchronous
machine i with respect to the synchronous rotating reference
frame.

00

k∈N

X  i

00
00
Bik (Edk
−
cos δik − Eqk
sin δik ) )

00 00
0
00
0
00
00
Tdi
Ėqi = Eqi
− Eqi
+ (Xdi
− Xdi
)(Bii Eqi
−
00 00
Tqi
Ėdi

X 

00
00
Bik (Edk
sin δik + Eqk
cos δik ) )

l∈E

V. PORT-HAMILTONIAN REPRESENTATION
Using the energy functions from the previous section, the
multi-machine model (6) can be put into a port-Hamiltonian
form. To this end, we derive expressions for the gradient of
each energy function.
3 The mapping from dq-reference frame k to dq-reference frame i in the
phasor domain is done by multiplication of e−jδik [10].

−

X 

00
00
Bik (Edk
sin δik + Eqk
cos δik ) )
k∈Ni

X 

00
00
Bik (Edk
cos δik − Eqk
sin δik ) )

k∈Ni

A. Transmission Line Energy
Recall that the energy stored in transmission line l between
internal buses i and k is given by (9). It can
P be verified that
the gradient of the total energy HL := l∈E Hl stored in
the transmission lines takes the form
 ∂HL  
  
00
E 00 Idi + Eqi
Iqi
Pei
∂δi
 ∂H00L   di
 =  Idi  ,
Idi
 ∂Eqi  =
∂HL
Iqi
Iqi
00
∂Edi

where Idi , Iqi are given by P
(5). Here it is used that the selfsusceptances satisfy Bii = k∈Ni Bik for all i ∈ V.
1) State transformation: In the sequel, it is more convenient to consider a different set of variable describing
the voltage angle differences. Define for each edge l ∈ E
ηl := δik where i, k are respectively the positive and negative
ends of l. In vector form we obtain η = DT δ ∈ Rm , and
observe that this implies
D

∂HL
∂Hp
=D
= Pe .
∂η
∂η

B. Electrical Energy SG
Further, notice that the electrical energy stored in the
equivalent circuits along the d- and q-axis of generator i
is given by (7) and satisfies
 " ∂Hedi #  0


0
0
0
Eqi − Ef i
Xdi − Xdi
Xdi − Xdi
∂Eqi
=
0
00
00
0
∂Hedi
0
Xdi
− Xdi
Eqi
− Eqi
00
∂Eqi
"
#

 ∂Heqi


0
0
0
0
Xqi − Xqi
Xqi − Xqi
Edi
∂Edi
.
0
00
00
0
∂Heqi =
0
Xqi
− Xqi
Edi
− Edi
∂E 00
di

By the previous observations, and by aggregating the states,
the dynamics of the multi-machine system can now be
written in the form (11) where the Hamiltonian is given
0
0
0
00
by (10) and X̂di := Xdi − Xdi
, X̂di
:= Xdi
− Xdi
, X̂d =
0
0
diagi∈V {X̂di } and X̂d , X̂q , X̂q are defined likewise. In addi0
tion, Td0 = diagi∈V {Tdi
} and Td0 , Tq , Tq0 are defined similarly.
Observe that the multi-machine system (11) is of the form
ẋ = (J − R)∇H(x) + gu
y = g T ∇H(x)

(12)

−D
0
0
0
0
0
0
−(Td0 )−1 X̂d
0
0
0
−(Tq0 )−1 X̂q
0
0
0
0
0
0
0

T
g= I 0 0 0 0 0 .

 
0
ṗ
 η̇  DT
 0 
 Ėq   0
 
ẋp = 
 Ė 0  =  0
 d 
Ė 00   0


q
Ėd00

y = g T ∇Hp ,

where J = −J T , R = RT are respectively the antisymmetric and symmetric part of the matrix depicted in (11).
Notice that the dissipation matrix of the electrical part is
positive definite (which implies R ≥ 0) if


0
X −X 0
Xdi −Xdi
2 diT 0 di
0
Tdi
di

∀i ∈ V,
0
0
00  > 0,
Xdi
−Xdi
Xdi −Xdi
2
T0
T 00
di

di

which, by invoking the Schur complement, holds if and only
if
0
00
0
0
00
4(Xdi
− Xdi
)Tdi
− (Xdi − Xdi
)Tdi
> 0,

∀i ∈ V. (13)

Note that a similar condition holds for the q-axis.
Proposition 1: Suppose that for all i ∈ V the following
holds:
0
00
0
4(Xdi
− Xdi
)Tdi
− (Xdi
0
00
0
4(Xqi − Xqi )Tqi
− (Xqi

0
00
− Xdi
)Tdi
0
00
− Xqi )Tqi

>0
> 0.

(14)

Then (11) is a port-Hamiltonian representation of the multimachine network (6).
It should be stressed that (14) is not a restrictive assumption
00
0
00
, Tqi

 Tdi
as it holds for a typical generator since Tdi
0
Tqi , see also Table 4.2 of [6] and Table 4.3 of [8].
Because the interconnection and damping structure J − R
of (11) is state-independent, the shifted Hamiltonian
H̄(x) = H(x) − (x − x̄)T ∇H(x̄) − H(x̄)

(15)

acts as a local storage function for proving passivity in a
neighborhood of a steady state x̄ of (12), provided that the
Hessian of H evaluated at x̄ (denoted as ∇2 H(x̄)) is positive
definite4 .
Proposition 2: Let ū be a constant input and suppose
there exists a corresponding steady state x̄ to (12) such
that ∇2 H(x̄) > 0. Then the system (12) is passive in a
neighborhood of x̄ with respect to the shifted external portvariables ũ := u − ū, ỹ := y − ȳ where ȳ := g T ∇H(x̄).
Proof: Define the shifted Hamiltonian by (15), then we
obtain
ẋ = (J − R)∇H(x) + gu

= (J − R)(∇H̄(x) + ∇H(x̄)) + gu
= (J − R)∇H̄(x) + g(u − ū)

= (J − R)∇H̄(x) + gũ

ỹ = y − ȳ = g T (∇H(x) − ∇H(x̄)) = g T ∇H̄(x).
4 Observe

that ∇2 H(x) = ∇2 H̄(x) for all x.

(16)

0
0
−(Td0 )−1 X̂d
0
00 −1 0
−(Td ) X̂d
0


0

0


0
 ∇Hp + g(Pm − Pd ),
−(Tq0 )−1 X̂q 


0
00 −1 0
−(Tq ) X̂q

(11)

As ∇2 H(x̄) > 0 we have that H̄(x̄) = 0 and H̄(x) > 0
for all x 6= x̄ in a sufficiently small neighborhood around x̄.
Hence, by (16) the passivity property automatically follows
where H̄ acts as a local storage function.
VI. MINIMIZING GENERATION COSTS
The objective is to minimize the total quadratic generation
cost while achieving zero frequency deviation. By analyzing
the steady states of (6), it follows that a necessary condition
for zero frequency deviation is 1T Pm = 1T Pd , i.e., the total
supply must match the total demand. Therefore, consider the
following convex minimization problem:
1 T
P QPm
2 m
s.t. 1T Pm = 1T Pd ,

min
Pm

(17)

where Q = QT > 0 and Pd is a constant unknown power
load.
Remark 5: Note that that minimization problem (17) is
easily extended to quadratic cost functions of the form
1 T
T
n
2 Pm QPm + b Pm for some b ∈ R . Due to space limitations, this extension is omitted.
As the minimization problem (17) is convex, it follows
that Pm is an optimal solution if and only if the KarushKuhn-Tucker conditions are satisfied [3]. Hence, the optimal
points of (17) are characterized by
∗
Pm
= Q−1 1λ∗ ,

λ∗ =

1 T Pd
T
1 Q−1 1

(18)

Next, based on the design of [14], consider a distributed
controller of the form
T θ̇ = −Lc θ − Q−1 ω

Pm = Q−1 θ − Kω

(19)

where T = diagi∈V {Ti } > 0, K = diagi∈V {ki } > 0 are
controller parameters. In addition, Lc is the Laplacian matrix of some connected undirected weighted communication
graph.
The controller (19) consists of three parts. Firstly, the term
−Kω corresponds to a primary controller and adds damping
into the system. The term −Q−1 ω corresponds to secondary
control for guaranteeing zero frequency deviation on longer
time-scales. Finally, the term −Lc θ corresponds to tertiary
control for achieving optimal production allocation over the
network.

Note that (19) admits the port-Hamiltonian representation
ϑ̇ = −Lc ∇Hc − Q−1 ω
Pm = Q−1 ∇Hc − Kω,

Hc =

1 T −1
ϑ T ϑ,
2

(20)

where ϑ := T θ. By interconnecting the controller (20) with
(11), the closed-loop system amounts to
  

 
ẋp
g
J − R − RK GT
=
∇H −
P
0 d
−G
−Lc
ϑ̇
(21)
 −1

0 0 0 0 0
G= Q
where J − R is given as in (11), H := Hp + Hc , and RK =
blockdiag(0, K, 0, 0, 0, 0). Define the set of steady states of
(21) by Ω and observe that any x := (xp , ϑ) ∈ Ω satisfies
the optimality conditions (18) and ω = 0.
Assumption 2: Ω 6= ∅ and there exists x̄ ∈ Ω such that
∇2 H(x̄) > 0.
Remark 6: While the Hessian condition of Assumption
2 is required for proving local asymptotic stability of (21),
guaranteeing that this condition holds can be bothersome.
However, while we omit the details, it can be shown that
∇2 H(x̄) > 0 if
the generator reactances are small compared to the
transmission line reactances.
• the subtransient voltage differences are small.
• the rotor angle differences are small.
Theorem 1: Suppose Pd is constant and there exists x̄ ∈ Ω
such that Assumption 2 is satisfied. Then the trajectories of
the closed-loop system (21) initialized in a sufficiently small
neighborhood around x̄ converge to the set of optimal points
Ω.
Proof: Observe by (16) that the shifted Hamiltonian
defined by (15) satisfies
•

H̄˙ = −(∇H̄)T blockdiag(R + RK , Lc )∇H̄ ≤ 0
where equality holds if and only if ω = 0, T −1 ϑ = θ = 1θ∗
for some θ∗ ∈ R, and ∇E H̄(x) = ∇E H(x) = 0. Here
∇E H(x) is the gradient of H with respect to the internal
voltages Eq0 , Ed0 , Eq00 , Ed00 . By Assumption 2 there exists a
compact neighborhood Υ around x̄ which is forward invariant. By invoking LaSalle’s invariance principle, trajectories
initialized in Υ converge to the largest invariant set where
H̄˙ = 0. On this set ω, η, θ, Eq0 , Ed0 , Eq00 , Ed00 are constant and,
T
more specifically, ω = 0, θ = 1λ∗ = 1 1T1QP−1d 1 corresponds
to an optimal point of (17) as Pm = Q−1 1λ∗ where λ∗
is defined in (18). We conclude that the trajectories of the
closed-loop system (21) initialized in a sufficiently small
neighborhood around x̄ converge to the set of optimal points
Ω.
Remark 7: While by Theorem 1 the trajectories of the
closed-loop system (21) converge to the set of optimal points,
it may not necessarily converge to a unique steady state as the
closed-loop system (21) may have multiple (isolated) steady
states.

VII. CONCLUSIONS
We have shown that a much more advanced multi-machine
model than conventionally used can be analyzed using the
port-Hamiltonian framework. Based on the energy functions
of the system, a port-Hamiltonian representation of the
model is obtained. Moreover, the system is proven to be
incrementally passive which allows the use of a passive
controller that regulates the frequency in an optimal manner,
even in the presence of an unknown constant demand.
The results established in this paper can be extended in
many possible ways. Current research has shown that the
third, fourth and fifth order model as given in [8] admit a
similar port-Hamiltonian structure as (11). It is expected that
the same controller as designed in this paper can also be used
in these lower order models.
While the focus in this paper is about (optimal) frequency
regulation, further effort is required to investigate the possibilities of (optimal) voltage control using passive controllers.
Another extension is to include transmission line resistances
of the network. Finally, one could look at the possibility to
extend the results to the case where inverters and frequency
dependent loads are included into the network as well.
R EFERENCES
[1] F. Alvarado, J. Meng, C. DeMarco, and W. Mota. Stability analysis of
interconnected power systems coupled with market dynamics. IEEE
Transactions on Power Systems, 16(4):695–701, November 2001.
[2] P. Anderson and A. Fouad. Power System Control and Stability. The
Iowa State Univsersity Press, first edition, 1977.
[3] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge
University Press, first edition, 2004.
[4] S. Caliskan and P. Tabuada. Compositional transient stability analysis
of multimachine power networks. IEEE Transactions on Control of
Network systems, 1(1):4–14, 2014.
[5] S. Fiaz, D. Zonetti, R. Ortega, J. Scherpen, and A. van der Schaft. A
port-Hamiltonian approach to power network modeling and analysis.
European Journal of Control, 19(6):477–485, December 2013.
[6] P. Kundur. Power System Stability and Control. Mc-Graw-Hill
Engineering, 1993.
[7] 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, pages 735–740. IEEE, 2014.
[8] J. Machowski, J. Bialek, and J. Bumby. Power System Dynamics:
Stability and Control. John Wiley & Sons, Ltd, second edition, 2008.
[9] R. Ortega, M. Galaz, A. Astolfi, Y. Sun, and T. Shen. Transient
stabilization of multimachine power systems with nontrivial transfer
conductances. Automatic Control, IEEE Transactions on, 50(1):60–75,
2005.
[10] J. Schiffer, D. Zonetti, R. Ortega, A. Stankovic, T. Sezi, and J. Raisch.
Modeling of microgrids-from fundamental physics to phasors and
voltage sources. arXiv preprint arXiv:1505.00136, 2015.
[11] 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.
[12] T. Stegink, C. De Persis, and A. van der Schaft. A port-Hamiltonian
approach to optimal frequency regulation in power grids. arXiv
preprint arXiv:1509.07318, 2015.
[13] T. Stegink, C. De Persis, and A. van der Schaft. Port-Hamiltonian
formulation of the gradient method applied to smart grids. In 5th
IFAC Workshop on Lagrangian and Hamiltonian Methods for Non
Linear Control, Lyon, France, July 2015.
[14] 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, 64:240–253, 2016.
[15] A. van der Schaft and D. Jeltsema. Port-Hamiltonian systems theory:
An introductory overview. Foundations and Trends in Systems and
Control, 1(2-3):173–378, 2014.

[16] X. Zhang and A. Papachristodoulou. A real-time control framework for
smart power networks: Design methodology and stability. Automatica,
58:43–50, 2015.
[17] 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), pages 1–6.
IEEE, 2015.

