Structure-preserving discretization for port-Hamiltonian
descriptor systems
Volker Mehrmann∗and Riccardo Morandin†
March 25, 2019
Abstract
We extend the modeling framework of port-Hamiltonian descriptor systems to include
under- and over-determined systems and arbitrary differentiable Hamiltonian functions.
This structure is associated with a Dirac structure that encloses its energy balance prop-
erties. In particular, port-Hamiltonian systems are naturally passive and Lyapunov stable,
because the Hamiltonian defines a Lyapunov function. The explicit representation of in-
put and dissipation in the structure make these systems particularly suitable for output
feedback control. It is shown that this structure is invariant under a wide class of non-
linear transformations, and that it can be naturally modularized, making it adequate for
automated modeling. We investigate then the application of time-discretization schemes to
these systems and we show that, under certain assumptions on the Hamiltonian, structure
preservation is achieved for some methods. Relevant examples are provided.
1 Introduction
The port-Hamiltonian (pH) structure [11] is very appealing for modeling, simulation and control
of complex multiphysics systems. PH systems generalize classical Hamiltonian systems and gra-
dient systems by including energy dissipation and interaction with the environment (represented
by ports). Implicit formulations for pH systems enclosing their properties are given through
objects from differential geometry, known as Dirac structures. Similarly as pure Hamiltonian
systems, pH systems also gain from the application of geometric integration methods [3], like the
Gauss-Legendre collocation schemes. With additional requirements imposed on the form of the
Hamiltonian function, structure preseration can be achieved, by a proper discretization of the
Dirac structure [4].
The energy concept can be used as a common language, to interconnect different pH systems
while mantaining the structure, even when they originate in different physical domains, that
may include mechanical, mechatronic, fluidic, thermic, hydraulic, pneumatic, elastic, plastic
or electric components. The explicit incorporation of constraints, often inavoidable when us-
ing modeling packages such as Modelica (https://www.modelica.org), Matlab/Simulink
(https://www.mathworks.org) or Simpack (http://www.simpack.com), produce differential-
algebraic equations (DAEs), also referred to as descriptor systems, which may contain hidden
∗Institut f¨ur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG. Email address:
†Institut f¨ur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG. Email address:
1
constraints, consistency requirements for initial conditions and additional regularity require-
ments. A definition for linear time-varying pH descriptor systems (pHDAEs) has been given in
[1]. While including many interesting examples, that description falls short of fully extending
conventional pH systems, since it is restricted to quadratic Hamiltonian functions and finite di-
mensional states, where the dimension is the same as the number of equations, and the matrix
coefficients are independent from state.
In this paper we extend the definition of pHDAEs to include arbitrary differentiable Hamilto-
nian functions and systems with a different number of states and equations. The new definition
extends to the case of infinite-dimensional states and time- and space-varying coefficients. We
also give a new definition of Dirac structure, so that we can always associate one to our pHDAEs.
We generalize the results in [4] to the case of pHDAEs, so that we can apply structure-preserving
time discretization to port-Hamiltonian descriptor systems.
The paper is organized as follows. In Section II, we introduce our new definition of pHDAEs,
and we investigate its properties. In Section III, we study the application of collocation schemes to
pHDAEs and conditions for which structure preservation is achieved. In Section IV, we illustrate
a simple example of a pHDAE from the electrical circuit domain, we exploit the pH structure
to apply control to the system and we present numerical experiments using Gauss-Legendre
collocation.
2 Port-Hamiltonian descriptor systems
2.1 Formulation
Definition 1 (pHDAE).Let us consider a time interval I⊆R, a state space X ⊆ Rnand the
extended state space S=I× X. A port-Hamiltonian descriptor system (or pHDAE) is a system
of differential (-algebraic) equations of the form
E(t, x) ˙x+r(t, x)=(J(t, x)−R(t, x))z(t, x)+(B(t, x)−P(t, x))u,
y= (B(t, x) + P(t, x))Tz(t, x)+(S(t, x)−N(t, x))u, (1)
associated with an Hamiltonian function H ∈ C1(S,R), where x(t)∈ X is the state, u(t), y(t)∈
Rmare the input and output, E∈ C(S,R`,n) is the flow matrix,r, z ∈ C(S,R`) are the
time-flow and effort functions, J, R ∈ C(S,R`,`) are the structure and dissipation matrices,
B, P ∈ C(S,R`,m) are the port matrices and S, N ∈ C(S,Rm,m) are the feed-through matrices.
Furthermore, the following properties must hold:
1. The extended structure and dissipation matrices Γ, W ∈ C(S,R`+m,`+m), defined as
Γ := J B
−BTN, W := R P
PTS(2)
satisfy Γ = −ΓTand W=WT≥0 pointwise.
2. The gradient of the Hamiltonian satisfies ∂xH=ETzand ∂tH=zTrpointwise.
This definition can be extended to the case of weak solutions and state spaces with infinite
dimension. In particular, this framework can be also used to describe partial-differential-algebraic
equations (PDAEs). For simplicity, the results and examples presented in this paper will be
limited to the finite-dimensional case.
2
Remark 1.The definition of pH system often includes extra conditions on the Hamiltonian,
e.g. that it is a convex function, or that it is always non-negative [1]. While these conditions are
not necessary in general, they are often satisfied, and they strengthen the stability properties
derived from the pH structure.
2.2 Properties
Port-Hamiltonian descriptor systems as defined in (1) satisfy many useful properties. We list
and prove some of these in the following.
2.2.1 Dissipation inequality
Any port-Hamiltonian system must satisfy some kind of dissipation inequality, expressing its
passivity, i.e. energy cannot be created within the system.
Theorem 1 (dissipation inequality).Let us consider a pHDAE of the form (1). Then the power
balance equation (PBE)
d
dtH(t, x(t)) = −z
uT
Wz
u+uTy(3)
holds along any solution x, for any input u. In particular, the dissipation inequality
H(t2, x(t2)) − H(t1, x(t1)) ≤Zt2
t1
u(τ)Ty(τ)dτ(4)
holds.
Proof. The pHDAE (1) can be written as
E˙x+r
0= (Γ −W)z
u+0
y,
in particular
d
dtH(t, x(t)) = ∂tH+∂xHT˙x=zT(E˙x+r) = z
uTE˙x+r
0=
=z
uT(Γ −W)z
u+0
y=−z
uT
Wz
u+uTy,
which is (3). By time integration, since W=WT≥0, we immediately get (4).
Note that, with no input (u= 0), if the Hamiltonian His locally positive-definite in an
equilibrium point x∗(up to shifting it by a constant), then His a Lyapunov function and the
system is stable.
2.2.2 Variable transformations
This class of pHDAEs is closed under many variable transformations:
3
Theorem 2. Consider a pHDAE of the form (1). Let ˜
X ⊆ R˜nbe a second state space, let
˜
S:= Iט
X, let x=ϕ(t, ˜x)∈ C1(˜
S,X)be a local diffeomorphism (with respect to ˜x) and let
U∈ C(˜
S,R`,`)be pointwise invertible. Consider the input-output DAE
˜
E˙
˜x+ ˜r= ( ˜
J−˜
R)˜z+ ( ˜
B−˜
P)u,
y= ( ˜
B+˜
P)T˜z+ (S−N)u, (5)
with ˜
E=UT(E◦ϕ)∂˜xϕ,˜
J=UT(J◦ϕ)U,˜
R=UT(R◦ϕ)U,˜
B=UT(B◦ϕ),˜
P=UT(P◦ϕ),
˜z=U−1(z◦ϕ)and ˜r=UT(r◦ϕ+(E◦ϕ)∂tϕ), where we denote (F◦ϕ)(t, ˜x) = F(t, ϕ(t, ˜x)) for any
F∈ C(S,·), and let ˜
H(t, ˜x) := (H◦ϕ)(t, ˜x). Then (5) is a pHDAE with Hamiltonian function ˜
H,
and to any solution (˜x, u, y)of (5) corresponds a solution (x, u, y)of (1) with x(t) = ϕ(t, ˜x(t)).
Furthermore, if ϕ(t, ·)is a global diffeomorphism for all t∈I, then the two systems are equivalent.
Proof. The transformed DAE system is obtained from (1) by setting x=ϕ(t, ˜x), pre-multiplying
with UTand inserting UU−1in front of zin the first equation. It is then clear that if (˜x, u, y)
is a solution of (5), then (x, u, y) is a solution of the original system. If ϕ(t, ·) is a global diffeo-
morphism, then we can apply the inverse tranformation (ϕ−1, U−1), where ϕ−1is to intended as
ϕ(t, ϕ−1(t, x)) = x, and get a solution (˜x, u, y) for any solution (x, u, y) of the original system,
making the two systems equivalent.
To show that (5) is still a pHDAE, we must check that conditions 1 and 2 in the definition
are satisfied.
1. By substitution, we get
˜
W=˜
R˜
P
˜
PTS=UT(R◦ϕ)U UT(P◦ϕ)
(P◦ϕ)TU S ◦ϕ=
=U0
0IT
(W◦ϕ)U0
0I≥0.
2. By substitution, we get
∂˜
H
∂˜x(t, ˜x) = ∂ϕ
∂˜x
T∂H
∂x ◦ϕ=∂ϕ
∂˜x
T
ETz◦ϕ=
=∂ϕ
∂˜x
T
(ET◦ϕ)UU−1(z◦ϕ) = ˜
ET˜z
and
∂˜
H
∂t (t, ˜x) = ∂H
∂t ◦ϕ+∂H
∂x ◦ϕT∂ϕ
∂t =
=zTr◦ϕ+zTE◦ϕ∂ϕ
∂t =
= (z◦ϕ)Tr◦ϕ+ (E◦ϕ)∂ϕ
∂t =
= ˜zTUTU−T˜r= ˜zT˜r.
This concludes the proof.
4
2.2.3 Autonomous form
Any DAE can be made autonomous by adding time as a state. In particular, any pHDAE can
be made autonomous without destroying the structure:
E r
0 1˙x
˙
t=J−R0
0 0z
0+B−P0
0 1u
1,
y
0=B+P0
0 1Tz
0+S−N0
0 0u
1.
(6)
Note that condition 2 becomes simply ∇˜x˜
H=˜
ET˜z, where the tilde denotes the quantities in
the autonomous system.
2.2.4 Structure-preserving interconnection
Let us consider two autonomous pHDAEs of the form
Ei˙xi= (Ji−Ri)zi+ (Bi−Pi)ui,
yi= (Bi+Pi)Tzi+ (Si−Ni)ui,
with Hamiltonian Hi, for i= 1,2, and assume that the aggregated input u= (u1, u2) and output
y= (y1, y2) satisfy a linear interconnection relation Mu +Ny = 0 for some M, N ∈ C(S,Rk,m).
Then the aggregated system can be written as a pHDAE of the form
E0 0
0 0 0
0 0 0
0 0 0
˙x
˙
ˆu
˙
ˆy
=
Γ−W0 0
Im−MT
0−Im
0M
0−NT
N0
z
ˆu
ˆy
0
+
0
0
Im
0
u,
y= ˆy,
with Hamiltonian H=H1+H2, where we have introduced new state variables ˆu, ˆy∈Rnthat
copy u, y, and we aggregated E= diag(E1, E2), x= (x1, x2), z= (z1, z2), m=m1+m2,
Γ = Π diag(Γ1,Γ2)ΠTand W= Π diag(W1, W2)ΠT, where Π ∈R`+m,`+mis the permutation
matrix
Π =
I`10 0 0
0 0 I`20
0In10 0
0 0 0 In2
.
Note that, in general, we may not be able to reduce the number of inputs and outputs. If we
also assume that the interconnection is energy-preserving (e.g. if Mu +Ny = 0 defines a Dirac
structure for (y, u)), then index reduction [5] and row operations can usually be applied to make
the system smaller.
2.3 Dirac structure
Port-Hamiltonian systems are ofter described through differential geometric structures known as
Dirac structures [11]. We present first the basic definitions.
Definition 2 (linear Dirac structure).Let Fbe a linear space and E:= F∗its dual space. Let
hh·,·ii be the bilinear form on F × E defined as
hh(f1, e1),(f2, e2)ii := he1|f2i+he2|f1i,
5
where h· | ·i denotes the duality pairing. A Dirac structure on F × E is then a linear subspace
D ⊆ F × E, such that D=D⊥⊥.
In particular, in finite dimension, one only needs to prove dim D= dim Fand he|fi= 0 for
all (f, e)∈ D. If (f, e)∈ D, then fand eare called flow and effort, respectively. In [11], the
more general definition of modulated Dirac structure over Xis also presented, as a subbundle
of TX ⊕ T∗X, that denotes the Whitney sum between the tangent and cotangent bundles of X.
To introduce a Dirac structure for pHDAEs, we need to extend this further.
Definition 3. Consider a state space Xand a vector bundle Vover Xwith fibers Vx. A Dirac
structure over Vis a subbundle D ⊆ V ⊕ V∗such that, for all x∈ X,Dx⊆ Vx× V∗
xis a linear
Dirac structure.
Note that modulated Dirac structures are a special case of our definition, where V=TX. To
associate a Dirac structure to our pHDAE system, we first prove the following lemma:
Lemma 1. Let D ⊆ V ⊕ V∗be a vector subbundle with fibers defined by
Dx={(f, e)∈ Vx× V∗
x:f+J(x)e= 0},
where J:X → L(V∗
x,Vx)is a skew-symmetric operator. Then Dis a Dirac structure.
Proof. For generic (f, e)∈ Dxand (f0, e0)∈ Vx× V∗
x,
hh(f, e),(f0, e0)ii =he|f0i+he0|fi=
=he|f0i−he0|Jei=he|f0+Je0i.
We show that (f0, e0)∈ Dxif and only if he|f0+Je0i= 0 for all (f, e)∈ Dx. On one hand, if
(f0, e0)∈ Dx, then he|f0+Je0i= 0 holds for any e∈ E. On the other hand, if (f0, e0)/∈ Dx, then
f0+Je06= 0 and so ∃e∈ E such that he|f0+Je0i= 1, but (f, e)∈ Dxwith f=−Je.
We can now associate a Dirac structure to our pHDAE system in the following way:
Theorem 3. Given an autonomous pHDAE, let us define the flow fiber Vx=Fs
x× Fp
x× Fd
x
for all x∈ X , where Fs
x:= E(x)TxX ⊆ R`is the storage flow fiber, Fp
x:= Rmis the port flow
fiber and Fd
x:= R`+mis the dissipation flow fiber. Let us partition f= (fs, fp, fd)∈ V and
e= (es, ep, ed)∈ V∗. Then the subbundle D ⊆ V ⊕ V∗with
Dx=((f, e)∈ Vx× V∗
x
f+Γ(x)I`+m
−I`+m0e= 0)
is a Dirac structure over V. Furthermore, the system of equations
fs=−E(x) ˙x, es=z(x),
fp=y, ep=u,
ed=−W(x)fd,(f, e)∈ Dx
(7)
is equivalent to the original pHDAE, and he|fi= 0 represents the power balance equation.
Proof. Dis a Dirac structure because of Lemma 1. Note that the pHDAE can be written in
compact form as
E(x) ˙x
−y= (Γ(x)−W(x)) z(x)
u.
6
The condition (f, e)∈ Dxcan be written as
−fs
−fp= Γ(x)es+ep, fd= (es, ep),
that together with the conditions (7) is equivalent to
E(x) ˙x
−y= (Γ(x)−W(x)) z(x)
u,
f=−E(x) ˙x, y, z(x)
u,
e=z(x), u, −W(x)z(x)
u,
which is exactly the compact form of the pHDAE, plus the definition for the flow and effort.
Finally, note that the equation he|fi= 0 can be written as
0 = hz(x)| − E(x) ˙xi+hu|yi+h−W(x)z(x)
u|z(x)
ui=
=d
dtH(x) + yTu−z(x)
uT
W(x)z(x)
u,
which is the power balance equation.
Note that, if we want to retrieve a pHDAE system from a Dirac structure, the additional
conditions (7) and the definition of H(x) are needed. These can also be lifted to a geometric
interpretation, by the means of a Lagrangian submanifold and of a dissipative structure [12].
3 Time discretization
Let us consider a finite-dimensional pHDAE of the form (1). Under some regularity assumptions
[5], many classes of Runge-Kutta methods can be applied to compute time-discretizations of
differential-algebraic equations. An important class of Runge-Kutta methods is the class of
collocation methods. We extend the results from [4] to pHDAEs, proceeding in a similar way.
3.1 Collocation methods for DAEs
Assume that an input function u:I× X → Rm, depending on time and possibly space, is given.
Given a time interval I= [t0, tf] of length h=tf−t0, and a consistent initial condition x0at time
t0, we approximate the solution x(t) of the pHDAE on Iwith a polynomial ˜x(t)∈R[t]s, where
R[t]sdenotes the vector space of polynomials of degree at most s. The polynomial is chosen such
that ˜x(t0) = x0, and that it satisfies the differential-algebraic equation in scollocation points
ti=t0+hγiwith γi∈[0,1] for i= 1, . . . , s.
Let `idenote the i-th Lagrange interpolation polynomial with respect to the nodes γ1, . . . , γs,
i.e.,
`i(τ) :=
s
Y
j=1
j6=i
τ−γj
γi−γj
.
7
Then we can write
˙
˜x(t0+τh) =
s
X
i=1
˙xi`i(τ),
˜x(t0+τh) = x0+h
s
X
j=1
˙xjZτ
0
`j(σ)dσ,
for certain ˙xi=˙
˜x(ti) that we have to compute, and also
xi:= ˜x(ti) = x0+h
s
X
j=1
αij ˙xj,
xf:= ˜x(tf) = x0+h
s
X
j=1
βj˙xj,
where αij := Rγi
0`j(σ)dσand βj:= R1
0`j(σ)dσ. In particular, the constants αij, βj, γi∈Rfor
i, j = 1 . . . s are the coefficients of the Butcher diagram of the associated Runge-Kutta method
[3].
3.2 Dirac structure associated to discretization
Consider for simplicity the autonomous case. The non-autonomous case is similar, with a more
cumbersome notation. Let Dxbe the Dirac structure associated to (6) (as in Theorem 3). We
define the Dirac structure associated to the time-discretization as {Dxi:i= 1, . . . , s}, i.e.,
Dxi=((fi, ei)∈ Vxi× V∗
xi
fi+Γ(ti, xi)I`+m
−I`+m0ei= 0),
with fi= (fi
s, yi, fi
d) and ei= (ei
s, ui, ei
d). Taking (fi, ei)∈ Dxi, together with
xf=x0+h
s
X
i=1
βi˙xi,(8a)
xi=x0+h
s
X
j=1
αij ˙xj,(8b)
fi
s=−E(xi) ˙xi,(8c)
ei
s=z(xi),(8d)
ei
d=−W(xi)fi
d,(8e)
ui=u(xi),(8f)
we get a system that is equivalent to applying the collocation method and computing the discrete
input and output ui, yi, for i= 1, . . . , s. Let us define the collocation flows, efforts, input and
8
output in R[t]s−1as
˜
fs(t0+hτ) =
s
X
i=1
fi
s`i(τ),˜es(t0+hτ) =
s
X
i=1
ei
s`i(τ),
˜
fd(t0+hτ) =
s
X
i=1
fi
d`i(τ),˜ed(t0+hτ) =
s
X
i=1
ei
d`i(τ),
˜y(t0+hτ) =
s
X
i=1
yi`i(τ),˜u(t0+hτ) =
s
X
i=1
ui`i(τ).
Note that, by construction, ( ˜
fs,˜y, ˜
fd; ˜es,˜u, ˜ed)∈ D˜xin all collocation points ti. Let us consider
the evolution of the Hamiltonian Halong the collocation polynomial ˜x(t). In particular, let
H(t) := H(˜x(t)): we can then write
H(t)−H(t0) = Zt
t0
˙
H(s)ds.
In the collocation points, the PBE is satisfied:
˙
H(ti) = ∇H(xi)T˙xi=z(xi)TE(xi) ˙xi=
=−hei
s|fi
si=hei
d|fi
di+hyi|uii,
for i= 1, . . . , s. In particular, if we apply the quadrature rule associated to the collocation
method, we get
H(tf)−H(t0) = h
s
X
j=1
βj˙
H(tj) + O(hp+1) =
=−h
s
X
j=1
βjhej
s|fj
si+O(hp+1) =
=h
s
X
j=1
βjhej
d|fj
di+h
s
X
j=1
βjhyj|uji+O(hp+1),
where p∈Nis the degree of exactness of the quadrature rule. We observe that, for the same
reason,
h
s
X
j=1
βjhej
d|fj
di=Ztf
t0
h˜ed|˜
fdi+O(hp+1),(9a)
h
s
X
j=1
βjhyj|uji=Ztf
t0
h˜y|˜ui+O(hp+1),(9b)
so
H(tf)−H(t0) = Ztf
t0h˜ed(s)|˜
fd(s)i+h˜y(s)|˜u(s)ids+O(hp+1).
We note that, if p≥2s−2, then equations (9a) and (9b) are exact. Furthermore, if βj≥
0 for j= 1, . . . , s (and this is the case for many collocation methods), we can deduce that
hPs
j=1hej
d|fj
di ≤ 0, thus the dissipation component retains its qualitative behaviour.
9
3.3 The quadratic Hamiltonian case
Let us consider the case where the Hamiltonian H(x) is a polynomial of degree (at most) 2, i.e. it
can be written as
H(x) = 1
2xTQx +vTx+c,
for some Q=QT∈Rn,n,v∈Rnand c∈R. Since ˜x∈R[t]n
s, we also have H=H ◦ ˜x∈R[t]2s
and ˙
H∈R[t]2s−1. It is known that the maximum degree of exactness for quadrature rules
with snodes is 2s−1, and that it is attained only with Gaussian quadrature rules, that are
associated with Gauss-Legendre collocation methods. In particular, if we apply such a method,
the integration of ˙
Hwill be exact, i.e.
H(tf)−H(t0) = h
s
X
j=1
βjhej
d|fj
di+h
s
X
j=1
βjhyj|uji=
=Ztf
t0h˜ed(s)|˜
fd(s)i+h˜y(s)|˜u(s)ids.
In particular, since we always have βj≥0 for Gauss-Legendre collocation, the dissipation term
is always non-positive, and we can write the discrete dissipation inequality
H(tf)−H(t0)≤h
s
X
j=1
βjhyj|uji=Ztf
t0
h˜y(s)|˜u(s)ids.
Thus, in the quadratic Hamiltonian case the pH structure is preserved, in the sense that the
PBE and the dissipation inequality have an exact discrete version.
4 Examples
4.1 A basic DC power network example
Let us consider as a toy example the linear electrical circuit represented in Fig. 1, where
RG, RL, RR>0 are resistances, L > 0 is an inductor, C1, C2>0 are capacitors and EGis
a controlled voltage source. This can be interpreted as a basic representation of a DC generator
(EG,RG), connected to a load (RR) with a transmission line (the πmodel given by C1, C2, L, RL).
By means of Kirchhoff’s circuit laws, this system can be naturally written as the following DAE:
−
+
EG
RG
IG
LIRL
RR
IR
C1
V1
I1
C2V2
I2
Figure 1: Basic DC power network example
10
L˙
I=−RLI+V2−V1,
C1˙
V1=I−IG,
C2˙
V2=−I−IR,
0 = −RGIG+V1+EG,
0 = −RRIR+V2.
(10)
The energy in the system can be stored in the inductor and in the two capacitors, giving the
Hamiltonian
H(I, V1, V2) = 1
2LI2+1
2C1V2
1+1
2C2V2
2.(11)
The system can then be written equivalently as an autonomous pHDAE of the form
E˙x= (J−R)x+Bu, (12a)
y=BTx, (12b)
with x= (I, V1, V2, IG, IR), E= diag(L, C1, C2,0,0), B=e4,u=EG,y=IGand
J=
0−1 1 0 0
1 0 0 −1 0
−1000−1
0 1 0 0 0
0 0 1 0 0
, R =
RL0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 RG0
0 0 0 0 RR
.
The power balance equation reads as
˙
H=−RLI2−RGI2
G−RRI2
R+IGEG;
in particular, if we shut down the generator (EG= 0), the Hamiltonian will decrease and converge
to a solution such that ˙
H= 0, that is, I=IG=IR= 0. The only state compatible with (10)
that satisfies that condition is x= 0, so the system will always converge to that asymptotically
stable point.
4.2 Controlling the circuit
Suppose now that RRrepresents a consumer, that requires a fixed amount of power P=RRI2
R
to be delivered to them. We would like then to control the voltage of the generator EG, so that
the state of the system will converge to IR=I∗
R:= −pP/RR. If we assume that a solution
with IR≡I∗
Rexists, then we also get I=IG≡ −I∗
R,V1≡(RR+RL)I∗
R,V2≡RRI∗
Rand
EG≡ −(RR+RL+RG)I∗
R.
This can be interpreted in the following way, exploiting the port-Hamiltonian framework: let
x∗=
I∗
V∗
1
V∗
2
I∗
G
I∗
R
:= rP
RR
1
−RR−RL
−RR
1
−1
denote the desired state. By applying the change of variables ˜x=x−x∗to (12), we get the
equivalent pHDAE
E˙
˜x= (J−R)(˜x+x∗) + e4u, (13a)
y=eT
4(˜x+x∗),(13b)
11
with the same Hamiltonian. Since our goal is having ˜x≡0 as an asymptotically stable solution,
and (J−R)x∗=−pP/RR(RR+RL+RG)e4=: −e4u∗, by construction, we write the equivalent
system
E˙
˜x= (J−R)˜x+e4˜u, (14a)
˜y=eT
4˜x, (14b)
with ˜u=u−u∗and ˜y=y−eT
4x∗=IG−I∗
G, which is again a pHDAE, but with Hamiltonian
˜
H=1
2L(I−I∗)2+1
2C1(V1−V∗
1)2+1
2C2(V2−V∗
2)2.(15)
As before, this shows that choosing ˜u= 0 (i.e. u=−u∗) would make the system converge to the
desired state, that will be asymptotically stable. On the other hand, this is not the only input
that would satisfy this goal. If we want to speed up the convergence, we have to increase the
dissipation. To do so, we can for example apply some linear feedback of the form ˜u=−α˜yfor
some α > 0 (i.e. u=u∗−α(IG−I∗
G)): in this way, the dissipation inequality can be strengthened,
since
˙
˜
H=−RL˜
I2−(RG+α)˜
I2
G−R˜
I2
R≤
≤ −RL˜
I2−RG˜
I2
G−R˜
I2
R.
From another point of view, if we apply a feedback of the form u=u∗−α(y−I∗
G) + ˆu, then
(14) can be written as
E˙
˜x= (J−Rα)˜x+e4ˆu, (16a)
˜y=eT
4˜x, (16b)
with Rα=R+αe4eT
4> R.
4.3 Numerical simulations
We present numerical simulations on the toy example (10). We choose as constants L= 2,
C1= 0.01, C2= 0.02, RL= 0.1, RG= 6 and RR= 3. This choice is not based on real world
data, but it is made to reflect the usual assumption that, for transmission lines, LRLC.
We discretize the system with respect to time, using the midpoint rule, which is the Gauss-
Legendre collocation method with s= 1 stages and order p= 2 for ODEs. Since (10) is a
semi-explicit DAE with index 1, the chosen method will also give convergence order 2 (see [5,
Theorem 5.16]).
First, we apply the time discretization to a system without control (EG= 0), starting from
consistent non-zero initial values (see Fig. 2). The evolution of the state and of the Hamiltonian
show the expected qualitative behaviour: after some time, they all converge to zero. Furthermore,
while in the first half of the simulation the state oscillates, the Hamiltonian always decreases
monotonically, following the discrete dissipation inequality.
As a second example, we apply control to the system, starting from almost-zero consistent
initial values, with the goal of delivering P= 10 power to the consumer. This corresponds to
the desired state
x∗=
I∗
V∗
1
V∗
2
I∗
G
I∗
R
:= r10
3
1
−3.1
−3
1
−1
≈
1.8257
−5.6598
−5.4772
1.8257
−1.8257
,
12
0 0.2 0.4 0.6 0.8 1
-6
-4
-2
0
2
4
6
8
time
state
I
V1
V2
IG
IR
H
Figure 2: Evolution of state and energy with no control. The solid lines represent the state. The
dashed line represents the Hamiltonian H.
and to the control value u∗= 9.1p10/3≈16.614. To simulate the fact that the change of input
does not happen instantly, we actually choose as control
u(t) = u∗arctan 5(t−0.5)+ 0.5,
so that the input increases smoothly and rapidly from 0 to u∗. In Fig. 3 one can see that, after
some initial oscillations, the state quickly converges to the desired values, represented by the
dotted lines. At the same time, the second Hamiltonian ˜
Hdecreases monotonically, converging
to zero, while the first Hamiltonian H, representing actual energy in the system, converges to a
positive value H∗.
5 Conclusions and future works
5.1 Conclusions
We have presented a new definition for port-Hamiltonian descriptor systems, generalizing the
one from [1] to include a larger class of equations. Extension to weak problems and to partial
differential equations is also possible. We verified that this definition satisfies several properties,
in particular a dissipation inequality, invariance under a large class of variable transformations,
and structure-preserving interconnection. We generalized the definition of Dirac structure and
we associated one to every system included in our formulation. We extended the results from
[4] to apply structure-preserving collocation schemes to port-Hamiltonian differential-algebraic
equations. Finally, we illustrated a simple example from the domain of electrical circuits that
can be written with our formulation, and we presented some numerical experiments.
13
0 0.5 1 1.5 2 2.5 3
-15
-10
-5
0
5
time
state
I
V1
V2
IG
IR
H
H1
Figure 3: Evolution of state and energy with control. The solid lines represent the state of the
system. The dotted lines represent the desired state. The cyan and purple dashed lines (H and
H1) represent the first and second Hamiltonian Hand ˜
H, respectively.
5.2 Future Works
Ongoing work is on the analysis of a larger class of numerical methods applied to pHDAEs,
including more general Runge-Kutta schemes and partitioned methods, in particular in the case of
semi-explicit DAEs with index 1. Structure-preserving model reduction and space discretization
for PDEs are also being considered. Future work will include further analysis of the properties
of pHDAEs, in particular the extension of the results from [7] and possibly the characterization
of controllability and observability of pH systems.
6 Acknowledgments
Riccardo Morandin is supported by the Deutsche Forschungsgemeinschaft (DFG) in the subpro-
ject B03, within the SFB Transregio 154. Volker Mehrmann is supported by the German Federal
Ministry of Education and Research (BMBF), in the project EiFer.
References
[1] Christopher Beattie, Volker Mehrmann, Hongguo Xu, and Hans Zwart. Linear port-
Hamiltonian descriptor systems. Mathematics of Control, Signals, and Systems, 30(4):17,
October 2018.
[2] S. Fiaz, D. Zonetti, R. Ortega, J.M.A. Scherpen, and A.J. van der Schaft. A port-
Hamiltonian approach to power network modeling and analysis. European Journal of Con-
trol, 19(6):477–485, December 2013. Relation: https://www.rug.nl/research/itm/ Rights:
University of Groningen, Research Institute of Technology and Management.
14
[3] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric Numerical Integration:
Structure-Preserving Algorithms for Ordinary Differential Equations; 2nd Ed. Springer,
Dordrecht, 2006.
[4] Paul Kotyczka and Laurent Lef`evre. Discrete-time port-Hamiltonian systems based on
Gauss-Legendre collocation. IFAC-PapersOnLine, 51(3):125–130, 2018. 6th IFAC Workshop
on Lagrangian and Hamiltonian Methods for Nonlinear Control LHMNC 2018.
[5] P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and Numerical
Solution. Z¨urich: European Mathematical Society Publishing House, 2006.
[6] C. Mehl, V. Mehrmann, and P. Sharma. Stability radii for linear Hamiltonian systems with
dissipation under structure-preserving perturbations. Bit Numerical Mathematics, 2017.
[7] Christian Mehl, Volker Mehrmann, and Michal Wojtylak. Linear Algebra Properties of
Dissipative Hamiltonian Descriptor Systems. SIAM Journal on Matrix Analysis and Appli-
cations, 39, January 2018.
[8] Volker Mehrmann, Riccardo Morandin, Simona Olmi, and Eckehard Sch¨oll. Qualitative
stability and synchronicity analysis of power network models in port-Hamiltonian form.
Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(10):101102, 2018.
[9] D Portillo, Juan Garcia Orden, and I Romero. Energy-Entropy-Momentum integration
schemes for general discrete non-smooth dissipative problems in thermomechanics. Interna-
tional Journal for Numerical Methods in Engineering, 112, February 2017.
[10] Mark Schiebl and Peter Betsch. Energy-Momentum-Entropy Consistent Numerical Methods
for Thermomechanical Solids Based on the GENERIC Formalism. March 2019.
[11] Arjan van der Schaft and Dimitri Jeltsema. Port-Hamiltonian Systems Theory: An In-
troductory Overview. Foundations and Trends R
in Systems and Control, 1(2-3):173–378,
2014.
[12] Arjan van der Schaft and Bernhard Maschke. Generalized port-Hamiltonian DAE systems.
Systems & Control Letters, 121:31–37, 2018.
15
