Port-Hamiltonian descriptor systems
Christopher Beattie∗and Volker Mehrmann†and Hongguo Xu‡and Hans Zwart §
May 29, 2017
Abstract
The modeling framework of port-Hamiltonian systems is systematically extended to
constrained dynamical systems (descriptor systems, differential-algebraic equations). A
new algebraically and geometrically defined system structure is derived. It is shown that
this structure is invariant under equivalence transformations, and that it is adequate
also for the modeling of high-index descriptor systems. The regularization procedure for
descriptor systems to make them suitable for simulation and control is modified to deal
with the port-Hamiltonian structure. The relevance of the new structure is demonstrated
with several examples.
Keywords: port-Hamiltonian system, descriptor system, differential-algebraic equation, pas-
sivity, stability, system transformation, differentiation-index, strangeness-index, skew-adjoint
operator.
AMS subject classification.: 93A30, 93B17, 93B11.
1 Introduction
Modeling packages such as modelica (https://www.modelica.org/), Matlab/Simulink
(http://www.mathworks.com) or Simpack [42] have come to provide excellent capabilities
for the automated generation of models describing dynamical systems originating in differ-
ent physical domains that may include mechanical, mechatronic, fluidic, thermic, hydraulic,
pneumatic, elastic, plastic, or electric components [1, 16, 20, 40, 41]. Due to the explicit
incorporation of constraints, the resulting systems comprise differential-algebraic equations
(DAEs), also referred to as descriptor systems in the system theory context. Descriptor sys-
tems may contain hidden constraints, consistency requirements for initial conditions, and
unexpected regularity requirements. Therefore, these models usually require further regular-
ization to be suitable for numerical simulation and control, see [11, 27, 30]. Our main focus
will be on linear-time varying descriptor systems, as they may arise from the linearization of
∗Department of Mathematics, Virginia Tech, Blacksburg, VA 24061, USA. [email protected]. Supported by
Einstein Foundation Berlin, through an Einstein Visiting Fellowship.
2Institut f¨ur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG.
[email protected]. Supported by Einstein Foundation Berlin via the Einstein Center ECMath
and by Deutsche Forschungsgemeinschaft via Project A02 within CRC 1029 ’TurbIn’
3Department of Mathematics, University of Kansas, Lawrence, KS 66045, USA. [email protected]. Partially
supported by Alexander von Humboldt Foundation and by Deutsche Forschungsgemeinschaft, through the DFG
Research Center Matheon Mathematics for Key Technologies in Berlin.
4University of Twente, Department of Applied Mathematics, P.O. Box 217, 7500 AE Enschede, The Nether-
lands. [email protected]
1
nonlinear DAE systems along a (non-stationary) reference trajectory, see [10]. These have
the form
E(t) ˙x=A(t)x+B(t)u,
y=C(t)x+D(t)u, (1)
together with an initial condition x(t0) = x0. The coefficient matrices E, A ∈C0(I,Rn,n),
B∈C0(I,Rn,m), C∈C0(I,Rm,n), and D∈C0(I,Rm,m), where we denote by Cj(I,X)
j∈ {0,1,2,3, . . .}the set of j-times continuously differentiable functions from a compact
time interval I= [t0, tf]⊆Rto X=Rn. If it is otherwise clear from the context, the
argument tof the coefficient functions is suppressed.
An important development in recent years has been to employ energy based modeling via
bond graphs [4, 12]. This has been implemented recently in 20-sim (http://www.20sim.com/),
for example. The resulting systems have a port-Hamiltonian (pH) structure, see e. g. [18, 24,
33, 37, 36], that encodes underlying physical principles such as conservation laws directly into
the structure of the system model. The standard form for pH systems appears as
˙x= (J−R)∇xH(x)+(B−P)u,
y= (B+P)T∇xH(x)+(S+N)u, (2)
where the function H(x) is the Hamiltonian which describes the distribution of internal
energy among energy storage elements of the system, J=−JT∈Rn,n is the structure matrix
describing energy flux among energy storage elements within the system; R=RT∈Rn,n is
the dissipation matrix describing energy dissipation/loss in the system; B±P∈Rn,m are
port matrices, describing the manner in which energy enters and exits the system, and S+N,
with S=ST∈Rm,m and N=−NT∈Rm,m, describes the direct feed-through from input to
output. It is necessary that
W="R P
PTS#≥0,(3)
where we write W > 0 (or W≥0) to assert that a real symmetric matrix Wis positive definite
(or positive semi-definite). Port-Hamiltonian systems generalize Hamiltonian systems, in the
sense that the conservation of energy for Hamiltonian systems is replaced by the dissipation
inequality:
H(x(t1)) − H(x(t0)) ≤Zt1
t0
y(t)Tu(t)dt. (4)
In the language of system theory, (4) shows that the dynamical system described in (2) is a
passive system [8]. Furthermore, H(x) defines a Lyapunov function for the unforced system,
so pH systems are implicitly Lyapunov stable [21]. Inequality (4) is an immediate consequence
of (3) and holds even when the coefficient matrices J,R,B,P,S, and Ndepend on xor
explicitly on time t, see [31], or when they are defined as linear operators acting on infinite
dimensional spaces [24, 39].
The physical properties of pH systems are encoded in the algebraic structure of the coeffi-
cient matrices and in geometric structures associated with the flow of the differential equation.
This leads to a remarkably robust modeling paradigm that greatly facilitates the combina-
tion and manipulation of pH systems. Note in particular that the family of pH systems is
closed under power-conserving interconnection (see [25]); model reduction of pH systems via
2
Galerkin projection yields (smaller) pH systems [2, 19, 35]; and conversely, pH systems are
easily extendable in the sense that new state variables can be included while preserving the
structure of (2), and so, the range of application of the model can be increased while ensuring
that the basic conservation principle (4) remains in force.
When state constraints are included in a pH system, the resulting system is a port-
Hamiltonian descriptor system (differential-algebraic equation) (pHDAE). pHDAE systems
arise also in singularly perturbed pH systems when small parameters are set to zero, see [38].
Significantly, there is no systematic way that has yet emerged to describe this problem class
consistently, in a way that reflects both the pH structure and the DAE structure accurately.
The first main topic of this paper is to propose such a systematic approach. This is a challeng-
ing task, in particular when constraints of the DAE are ’hidden’, which is often signaled with
the terminology ’high-index DAE’ [5, 27, 30]. Such DAEs are not well-suited for numerical
simulation and control and so, either a reformulation or a regularization of the model must
first be carried out, [11, 27]. We will briefly summarize the fundamentals of this technique in
Section 4.
It is sometimes stated in the literature, see e. g. [38], that port-Hamiltonian DAEs are
of differentiation-index at most one, i. e., that they do not contain hidden constraints arising
from derivatives. In contrast, we will show that higher-index pHDAEs are actually very com-
mon and so a regularization procedure is necessary. Unfortunately, the usual regularization
strategies do not preserve a given pHDAE structure of the model and so, how one should go
about this task while respecting pHDAE structure is the second main topic of the paper.
The paper is organized as follows. In Section 2 we give a definition of port-Hamiltonian
differential-algebraic systems and demonstrate that this is a relevant class for many appli-
cations. The main properties of this new class of pHDAE systems (such as stability and
dissipativity) are discussed in Section 3. Section 4 extends the definition to the nonlinear
case. The analysis of ‘index at most one’ pHDAEs is discussed in Section 5 while the struc-
tured regularization procedure is discussed in Section 6.
2 Linear port-Hamiltonian Differential-Algebraic Equations
In this section we introduce a new definition of systems of port-Hamiltonian descriptor sys-
tems (pHDAEs). Our new definition is slightly different from the concepts discussed in [38]
and is based on the concept of skew-adjoint differential-algebraic operators, see [29] for the
corresponding self-adjoint case.
Definition 1 A (differential-algebraic) operator
L:= Ed
dt − A : Ω ⊂C1(I,Rn)→C0(I,Rn)
with coefficient functions E ∈ C1(I,Rn,n),A ∈ C(I,Rn,n)is called skew-adjoint, if ET(t) =
E(t)and ˙
E(t) = −(A(t) + AT(t)) for all t∈I.
This definition is motivated by the following observation: starting with vector functions
x1(t), x2(t) that are absolutely continuous on the interval I= (t0, tf) each with square inte-
grable derivative and xi(t0) = xi(tf) = 0 for i= 1,2, and then denoting the usual L2inner
3
product as hx1, x2i=Rtf
t0xT
2x1dt, we have
hx1,L(x2)i=hx1,E˙x2− Ax2i=hx1,d
dt(Ex2)− Ax2−˙
Ex2i
=xT
2Ex1|tf
t0− hET˙x1, x2i−h(AT+˙
ET)x1, x2i
=h−ET˙x1−(AT+˙
ET)x1, x2i=h−E ˙x1+Ax1, x2i.
So formally, the adjoint operator L∗satisfies L∗=−L. Note the boundary terms arising
in partial integration will vanish under a wide variety of conditions replacing the requirement
of zero end conditions on x1(t) and x2(t).
Skew-adjoint operators stay skew-adjoint under time-varying congruence transformations.
Lemma 2 Consider a skew-adjoint differential-algebraic operator
L:= Ed
dt − A : Ω ⊂C1(I,Rn)→C0(I,Rn)
with coefficient functions E ∈ C1(I,Rn,n)and A ∈ C(I,Rn,n). Then for every V ∈ C1(I,Rn,r),
the operator LVdefined by
LV(x) := VTEV ˙x−(VTAV − VTE˙
V)x
is again skew-adjoint, i. e., L∗
V=−LV.
Proof. Since VTEV = (VTEV)T, it remains to consider the coefficient of x. Using ET=E
and ˙
E=−(A+AT), we have
d
dt(VTEV) = ˙
VTEV +VT˙
EV +VTE˙
V
=˙
VTEV − VT(A+AT)V+VTE˙
V
=−(VTAV − VTE˙
V)−(VTAV − VTE˙
V)T.
It should be noted that for any t∈Iand x∈C(I,Rn) we have LV(x(t)) = VT(t)L(V(t)x(t)).
Remark 3 Note that in Lemma 2 we do not need that the transformation matrix Vis
invertible. This implies, in particular, that with a projection matrix
V=Ir0
0 0
the projected system is still skew-adjoint.
Using the definition of skew-adjoint differential-algebraic operators we now present a definition
of pHDAEs.
Definition 4 A linear variable coefficient descriptor system of the form
E˙x= [(J−R)Q−EK]x+ (B−P)u,
y= (B+P)TQx + (S+N)u, (5)
with E, Q ∈C1(I,Rn,n),J, R, K ∈C0(I,Rn,n),B, P ∈C0(I,Rn,m),S=ST, N =−NT∈
C0(I,Rm,m), is called port-Hamiltonian descriptor system (port-Hamiltonian differential-
algebraic system) (pHDAE) if the following properties are satisfied:
4
i) the differential algebraic operator
L:= QTEd
dt −(QTJQ −QTEK) : D⊂C1(I,Rn)→C0(I,Rn) (6)
is skew-adjoint, i. e. we have that QTE∈C1(I,Rn,n)and for all t∈I,
QT(t)E(t) = ET(t)Q(t),and
d
dt(QT(t)E(t)) = QT(t)[E(t)K(t)−J(t)Q(t)] + [E(t)K(t)−J(t)Q(t)]TQ(t);
ii) the matrix function QTEis bounded from below by a constant symmetric matrix H0,
i. e., QT(t)E(t)−H0≥0for all t∈I;
ii) the matrix function
W:= "QTRQ QTP
PTQ S #∈C0(I,Rn+m,n+m) (7)
is positive semidefinite, i. e., W(t) = WT(t)≥0for all t∈I.
The associated Hamiltonian is defined as
H(x) := 1
2xTQTEx :C1(I,Rn)→R.(8)
Besides the matrix function Ein front of the derivative and the different definition of the
Hamiltonian, which gives the option of having singular matrices Eand Q, a major difference
to the definition of standard pH systems is the extra additive term −EKx on the right hand
side of (5), which is needed to accommodate time-varying changes of basis. Note further that
in this definition no further properties of the differential-algebraic operator are assumed, in
particular it is not assumed that it has a certain index as a differential-algebraic equation.
The assumption that the matrix function QTEis bounded by a constant matrix H0from
below implies that the Hamiltonian His bounded from below by a constant. This constant
is irrelevant when the derivative of His considered, but it guarantees that the Hamiltonian
can be interpreted as energy in a real physical system.
Example 5 Consider the model of a simple RLC network, see e. g. [13, 17], given by a linear
constant coefficient DAE
GcCGT
c0 0
0L0
0 0 0
| {z }
:=E
˙
V
˙
Il
˙
Iv
=
−GrR−1
rGT
r−Gl−Gv
GT
l0 0
GT
v0 0
| {z }
:=(J−R)I
V
Il
Iv
,(9)
with real symmetric constant matrices L > 0, C > 0, Rr>0 describing inductances, ca-
pacitances, and resistances, respectively that are present in the network. Here, Gvis of full
column rank, and the subscripts r, c, l, and vrefer to edge quantities corresponding to the re-
sistors, capacitors, inductors, and voltage sources, while V,Idenote the voltage and current,
respectively, on or across the branches of the given RLC network. This model has a pHDAE
5
structure with vanishing B, P, S, N, K, the matrix Qis the identity, E=ET,J=−JT,
QTRQ =R≥0, and
H=
V
Il
Iv
T
E
V
Il
Iv
=V
IlTGcCGT
c0
0LV
Il.
Example 6 In [14, 15] the propagation of pressure waves on the acoustic time scale in a
network of gas pipelines is considered and an infinite-dimensional pHDAE is derived. A
structure preserving mixed finite element space discretization leads to a block-structured
pHDAE system
E˙x= (J−R)Qx +Bu,
y=BTQx, (10)
x(t0) = x0,
with Q=I,P= 0, S+N= 0,
E=
M10 0
0M20
0 0 0
, J =
0−˜
G0
˜
GT0˜
NT
0−˜
N0
, R =
000
0˜
D0
000
B=
0
˜
B2
0
, x =
x1
x2
x3
,
where the vector valued functions x1:R→Rn1, x2:R→Rn2represent the discretized
pressure and flux, respectively, and x3:R→Rn3represents the Lagrange multiplier for
satisfying the space-discretized constraints. The coefficients M1=MT
1,M2=MT
2, and
˜
D=˜
DTare positive definite, and the matrices ˜
Nand ˜
GT˜
NTThave full row rank. The
Hamiltonian is given by H(x) = xTETQx =xT
1M1x1+xT
2M2x2.
Definition 4 brings the pH modeling framework and the DAE framework together in a struc-
tured way. It should be noted, however, that in a DAE we may have hidden constraints that
arise from differentiations, which are not explicitly formulated and the formulation of the DAE
that is used in simulation and control is not unique. One can for example add derivatives
of constraints which leads to an over-determined system, then one can add dummy variables
or Lagrange multipliers to make the number of variables equal to the number of equation
or one can remove some of the dynamical equations to achieve this goal, see [5, 16, 27, 30]
for detailed discussions on this topic. To rewrite these different formulations in the pHDAE
formulation is not always obvious. Let us demonstrate this with an example from multi-body
dynamics.
Example 7 A benchmark example for a nonlinear DAE system is the model of a two-
dimensional three-link mobile manipulator, see [6, 22], which is modeled as
˜
M(Θ)¨
Θ + ˜
C(Θ,˙
Θ) + ˜
G(Θ) = ˜
B1˜u+ ΨTλ,
ψ(Θ) = 0,(11)
where Θ = Θ1Θ2Θ3Tis the vector of joint displacements, ˜uis vector of control torques
at the joints, ˜
Mis mass matrix, ˜
Cis the vector of centrifugal and Coriolis forces, and ˜
Gis the
6
gravity vector. The term ΨTλwith Ψ = ∂ψ
∂Θis the generalized constraint force with Lagrange
multiplier λassociated with the constraint
ψ(Θ) = l1cos(Θ1) + l2cos(Θ1+ Θ2) + l3cos(Θ1+ Θ2+ Θ3)l3−l
Θ1+ Θ2+ Θ3= 0.
Besides the explicit constraint this system contains the first and second time derivative of ψ
as hidden algebraic constraints, see e. g. [16, 27]. There are several regularization procedures
that one can employ to make the system better suited for numerical simulation and control.
One possibility is to replace the original constraint by its time derivative Ψ(Θ) ˙
Θ = 0. In
this case the model equation can easily be written in a pHDAE formulation. Using Cartesian
coordinates for positions p, scaling the constraint equation by −1, and linearizing around a
non-stationary reference solution yields a linear time-varying DAE of the form
˜
Mδ¨p=−˜
Dδ ˙p−˜
Sδp +˜
GTδλ +˜
B1δu,
0 = −˜
Gδ ˙p, (12)
with pointwise symmetric positive definite matrix functions ˜
M, ˜
Sand pointwise symmetric
and positive semidefinite ˜
D. Adding a tracking output of the form y=˜
BT
1δ˙p, see e. g. [23],
and transforming to first order by introducing
x=
x1
x2
x3
:=
δ˙p
δp
δλ
, u =δu,
one obtains a linear time-varying pHDAE system E˙x= (J−R)Qx +Bu,y=BTQx, with
E:=
˜
M0 0
0I0
0 0 0
, R :=
˜
D0 0
0 0 0
0 0 0
, Q :=
I0 0
0˜
S0
0 0 I
,
J:=
0−I˜
GT
I0 0
−˜
G0 0
, B :=
˜
B1
0
0
, P = 0, S +N= 0.
The Hamiltonian in this case is given by H(x) = x1
x2T˜
M0
0˜
Sx1
x2.
Since the Lagrange multipliers in the multibody framework can be interpreted as external
forces, it is also possible to incorporate them in the input (B−P)uto achieve a pHDAE
formulation as in Definition 4, but also other formulations are possible, e. g. we can keep the
original algebraic constraint as well and use an extra Lagrange multiplier for the first time
derivative.
Remark 8 A special case of (5) takes the following form:
E˙x= (J−R)x+ (B−P)u,
y= (B+P)Tx+ (S+N)u, (13)
where E=ET∈C1(I,Rn,n),J=−JT, R =RT, K ∈C0(I,Rn,n),B, P ∈C0(I,Rn,m),
S=ST, N =−NT∈C0(I,Rm,m)as before but now we require,
7
i) the differential algebraic operator
L:= Ed
dt −J:D⊂C1(I,Rn)→C0(I,Rn) (14)
is skew-adjoint, so that we have for all t∈I,
d
dtE(t) = −J(t) + J(t)T;
ii) E(t)is positive semidefinite: E(t)≥0for all t∈I; and
ii) W(t) := "R(t)P(t)
PT(t)S(t)#≥0for all t∈I.
The effective Hamiltonian is now
H(x) := 1
2xTEx :C1(I,Rn)→R.(15)
Notice that in this model description we have merged the roles of Qand E. This is always
possible when Qis pointwise invertible, see Section 3 but this formulation may not be possible
when Qis singular.
3 Properties of pHDAE systems
To analyze the properties of pHDAE systems, we first derive the conservation of energy and
the dissipation inequality.
Theorem 9 A linear time-varying pHDAE system has the following properties:
i) If W≡0in (7) then d
dt H=uTy.
ii) The system satisfies the dissipation inequality (4).
Proof. By Definition 4 we have
d
dtH=1
2˙xT(QTE)x+xTd
dt(QTE)x+xT(QTE) ˙x
=1
2xTd
dt(QTE)x+xTQT(E˙x)
=1
2xTd
dt(QTE)x+xTQT([JQ −RQ −EK]x+Bu −P u)
=1
2xTd
dt(QTE)x+xTQTJQx −xTQTRQx −xTQTEKx +xTQTP u +uTBTQx
=1
2xTd
dt(QTE)x+xTQTJQx −xTQTRQx −xTQTEKx −xTQTP u
+uT(y−PTQx −Su −Nu)
=uTy+1
2xTd
dt(QTE)x+xTQTJQx −xTQTRQx −xTQTEKx
−xTQTPu −uTPTQx −uTSu
=uTy+1
2xTd
dt(QTE)x+xT[QT(JQ −EK)+(JQ −EK)TQ]x
−x
uT
Wx
u.
8
From the skew-adjointness of Lwe have that
d
dtH=uTy−x
uT
Wx
u.
Part i) then follows immediately from the assumption W≡0, while in Part ii) the fact that
W(t)≥0 for all t∈Iimplies that for any t1≥t0,
H(x(t1)) − H(x(t0)) = 1
2Zt1
t0
d
dtHdt ≤Zt1
t0
yTu dt.
An important feature of pHDAE systems is that a change of basis and a scaling with an
invertible matrix function preserves the pHDAE structure and the Hamiltonian.
Theorem 10 Consider a pHDAE system of the form (5) with Hamiltonian (8). Let U∈
C0(I,Rn,n)and V∈C1(I,Rn,n)be pointwise invertible in I. Then the transformed DAE
˜
E˙
˜x= [( ˜
J−˜
R)˜
Q−˜
E˜
K]˜x+ ( ˜
B−˜
P)u,
y= ( ˜
B+˜
P)T˜
Q˜x+ (S+N)u
with
˜
E=UTEV, ˜
Q=U−1QV, ˜
J=UTJU,
˜
R=UTRU, ˜
B=UTB, ˜
P=UTP,
˜
K=V−1KV +V−1˙
V , x =V˜x
is still a pHDAE system with the same Hamiltonian ˜
H(˜x) = 1
2˜xT˜
QT˜
E˜x=H(x).
Proof. The transformed DAE system is obtained from the original DAE system by setting
x=V˜xin (5), by pre-multiplying with UTand by inserting UU−1in front of Q. The
transformed operator corresponding to Lin (14) is
LV:= ˜
QT˜
Ed
dt −˜
QT(˜
J˜
Q−˜
E˜
K).
Because
˜
QT˜
E=VTQTEV, ˜
QT˜
J˜
Q=VTQTJQV, ˜
QT˜
EV −1˙
V=VTQTE˙
V ,
by Lemma 2, LVis skew-adjoint since Ldefined in (14) is skew-adjoint. Hence,
˜
QT˜
E=˜
ET˜
Q,
d
dt(˜
QT˜
E) = −˜
QT(˜
J˜
Q−˜
E˜
K)−(˜
J˜
Q−˜
E˜
K)T˜
Q.
It is straightforward to show that
d
dt ˜
H(˜x) = yTu−˜x
uT
˜
W˜x
u,
9
where
˜
W=˜
QT˜
R˜
Q˜
QT˜
P
˜
PT˜
Q S =VTQTRQV V TQTP
PTQV S
=V0
0IT
WV0
0I,
and Wis defined in (7). Because W(t) is positive semidefinite for all t∈I, so is ˜
W(t).
Therefore, for any t1≥t0,
˜
H(˜x(t1)) −˜
H(˜x(t0)) ≤Zt1
t0
yT(t)u(t)dt,
which establishes the dissipation inequality. Since
˜
QT(t)˜
E(t)−VT(t)H0V(t) = VT(t)(QT(t)E(t)−H0)V(t)≥0 for all t∈I,
and since Vis continuous, and thus VTH0Vis bounded on I, it follows that there exist a
constant symmetric matrix ˜
H0such that ˜
QT(t)˜
E(t)≥˜
H0for all t∈I.
An important point to note is that the Hamiltonian stays invariant under time-varying
changes of basis and the operator LV, the Hamiltonian ˜
H(˜x), and the matrix function ˜
Ware
independent of the choice of the matrix function U.
As we have already pointed out, our definition of pHDAE systems has the extra term
−EKx on the right hand side which is needed to incorporate time-varying changes of basis.
Even if K= 0 in the original system, after the transformation given in Theorem 10 the extra
term −˜
E˜
Kwith ˜
K=V−1˙
Vwill appear. Note that if an orthogonal change of basis is carried
out in a system with K= 0 then the resulting ˜
K=V−1˙
Vis skew-symmetric. Furthermore,
even if K6= 0, this term can be removed via a change of basis transformation which does not
change the Hamiltonian.
Lemma 11 Consider a pHDAE system
˜
E˙
˜x= [( ˜
J−˜
R)˜
Q−˜
E˜
K)]˜x+ ( ˜
B−˜
P)u,
y= ( ˜
B+˜
P)T˜
Q˜x+ (S+N)u
with Hamiltonian ˜
H(˜x) = 1
2˜xT˜
QT˜
E˜x, where ˜
K∈C(I,Rn,n). If V˜
K∈C1(I,Rn,n)is a point-
wise invertible solution of the matrix differential equation ˙
V=˜
KV with initial condition
V(t0) = I, then defining
E=˜
EV −1
K, Q =˜
QV −1
k,
J=˜
J, R =˜
R, B =˜
B,
P=˜
P, ˜x=V−1
Kx,
the system
E˙x= (J−R)Qx + (B−P)u,
y= (B+P)TQx + (S+N)u
is again pHDAE with the same Hamiltonian H(x) = ˜
H(˜x) = 1
2xTQTEx.
10
Proof. For a given matrix function ˜
K, the system ˙
VK=˜
KVKalways has a solution VK
that is pointwise invertible. The remainder of the proof follows by reversing the proof of
Theorem 10 with U=I.
Remark 12 Note that if Kis real and skew-symmetric, then the matrix function VKin
Lemma 11 can be chosen to be pointwise real orthogonal.
Following Theorem 10, if Eis pointwise invertible, then the original system can be trans-
formed into the one with new ˆ
Ebeing the identity, so into a standard port-Hamiltonian
system; and whenever Qis pointwise invertible, then the original system can be transformed
into the one with new ˆ
Qbeing the identity. Which of these formulations is preferable will
depend on the sensitivity (conditioning) of these transformations. In the context of numer-
ical simulation and control methods, these transformations should be avoided if they are
ill-conditioned.
4 Nonlinear DAEs and pHDAEs
In this section we briefly recall the theory of nonlinear DAE systems and then extend these
results to pHDAEs. Consider a general descriptor system of the form
F(t, x, ˙x, u)=0,
x(t0) = x0
y=G(t, x, u).(16)
Assume that F∈C0(I×Dx×D˙x×Du,Rn) and G∈C0(I×Dx×Du,Rm) are sufficiently
smooth, and that Dx,D˙x⊆Rn, and Duare open sets. Note that (in order to deal with
pHDAEs) in contrast to the more general case in [11], we assume square systems with an
equal number of equations and variables and with an equal number of inputs and outputs.
For the analysis and the regularization procedure we make use of the behavior approach
[34], which introduces a descriptor vector v= [xT, uT]T. We could also include the output
vector yin v, but in the context of pHDAEs it is preferable to keep the output equation
separate. The behavior formulation has the form
F(t, v, ˙v)=0,(17)
with F ∈ C0(I×Dv×D˙v,Rn) together with a set of initial conditions c(v(t0)) = v0which
results from the original initial condition. Note that although no initial condition is given for
uin the context of the regularization procedure discussed in [11] such conditions may arise,
so we formally state a condition for v(t0).
To regularize DAEs for numerical simulation and control, see [9, 11, 27], one uses the
behavior system (17) and some or all of its derivatives to produce an equivalent system
with the same solution set (all variables keep their physical interpretation), but where all
explicit and hidden constraints are available. The approach of [11] (adapted for the analysis
of pHDAEs) uses the state equation of (17) to form a derivative array, see [9],
Fµ(t, v, ˙v, . . . , v(µ+1))=0,(18)
which stacks the equation and its time derivatives up to level µinto one large system. We
denote partial derivatives of Fµwith respect to selected variables ζfrom (t, v, ˙v, . . . , v(µ+1))
11
by Fµ;ζ, and the solution set of the nonlinear algebraic equation associated with the derivative
array Fµfor some integer µ(considering the variables as well as their derivatives as algebraic
variables) by Lµ={vµ∈I×Rn+m×. . . ×Rn+m| Fµ(vµ)=0}.
The main assumption for the analysis is that the DAE satisfies the following hypothesis,
which in the linear case can be proved as a Theorem, see [27].
Hypothesis 13 Consider the system of nonlinear DAEs (17). There exist integers µ,r,a,
d, and νsuch that Lµis not empty and such that for every v0
µ= (t0, v0,˙v0, . . . , v(µ+1)
0)∈Lµ
there exists a neighborhood in which the following properties hold.
1. The set Lµ⊆R(µ+2)(n+m)+1 forms a manifold of dimension (µ+ 2)(n+m)+1−r.
2. We have rank Fµ;v, ˙v,...,v(µ+1) =ron Lµ.
3. We have corank Fµ;v, ˙v,...,v(µ+1) −corank Fµ−1;v, ˙v,...,v(µ)=νon Lµ, where the corank is
the dimension of the corange and the convention is used that corank of F−1;vis 0.
4. We have rank Fµ; ˙v,...,v(µ+1) =r−aon Lµsuch that there exist smooth full rank matrix
functions Z2and T2of size (µ+1)n×aand (n+m)×(n+m−a), respectively, satisfying
ZT
2Fµ; ˙v,...,v(µ+1) = 0,rank ZT
2Fµ;v=a, and ZT
2Fµ;vT2= 0 on Lµ.
5. We have rank F˙vT2=d=n−a−νon Lµsuch that there exists a smooth full rank
matrix function Z1of size n×dsatisfying rank ZT
1F˙vT2=d.
The smallest µfor which Hypothesis 13 holds is called the strangeness-index of (17), see
[27]. It generalizes the concept of differentiation-index [5] to over- and under-determined
systems but in contrast to the differentiation-index, ordinary differential equations and purely
algebraic equations have µ= 0 and for other systems the differentiation-index (if defined) is
µ+ 1, see [27]. The quantity νgives the number of trivial equations 0 = 0 in the system.
Of course, these equations can be simply removed and so for our further analysis we assume
that ν= 0.
If Hypothesis 13 holds then, locally (via the implicit function theorem) there exists, see
[26, 27], a system (in the same variables)
ˆ
F1(t, v, ˙v)=0,
ˆ
F2(t, v)=0,(19)
in which the first dequations ˆ
F1=ZT
1Fform a (linear) projection of the original set of
equations representing the dynamics of the system, while the second set ˆ
F2(t, v) = 0 of a
equations contains all explicit and hidden constraints and can be used to parameterize the
solution manifold and to characterize when an initial condition is consistent. Adding again
the output equation and writing (19) in the original variables we obtain the system
ˆ
F1(t, x, ˙x, u)=0,
ˆ
F2(t, x, u)=0,(20)
x(t0) = x0
y=G(t, x, u),
It should be noted that although formally also derivatives of uhave been used to form the
derivative array, no derivatives of uappear in the regularized system (20). This has been
12
shown in various contexts [7, 27, 28] and is due to the fact that only derivatives of equations
where Fµ;u≡0 are needed to generate (20). This means, in particular, that the equations in
ˆ
F2(t, x, u) = 0 can be partitioned further into equations that arise from the original system,
which include those algebraic equations in the original system which are explicit constraints
(in the behavior sense) so that the system can be made to be of differentiation index at
most one by a feedback (the part that is impulse controllable or controllable at infinity), and
implicit hidden constraints arising from differentiations of equations for which Fµ,u ≡0 in the
derivative array (the parts that are not impulse controllable).
Using these observations, the regularized system can be (locally in the nonlinear case)
written as
ˆ
E1˙x=ˆ
A1x+B1u,
0 = ˆ
A2x+B2u, (21)
0 = ˆ
A3x,
x(t0) = x0,
y=Cx +Du,
where the third equation that is representing all the hidden algebraic constraints of differenti-
ation index larger that one. Performing an appropriate (local) change of basis one can identify
some (transformed variables) which vanish and the remaining system consisting of the first
two equations is of index at most one in the behavior sense, see [7, 27, 28] for details. For the
first two equations in (20) and (21) one can always find an initial feedback u=k(x) + ˜uso
that the resulting system is strangeness-free (of differentiation-index one) as a system with
input ˜u= 0, see [3, 11] for a detailed analysis and regularization procedures. In the following
we assume that this reinterpretation has been done, so that the n×nmatrix (functions)
ˆ
E1
ˆ
A2
ˆ
A3
,
(ˆ
F1)˙x(t, x, ˙x, u)
(ˆ
F2)x(t, x, u)
(ˆ
F3)x(t, x)
,(22)
respectively, are locally invertible, see [27]. Furthermore there exists a (local) partitioning of
the variables so that the strangeness-free formulation takes the form
ˆ
E11 ˆ
E12 ˆ
E13
0 0 0
0 0 0
˙x1
˙x2
˙x3
=
ˆ
A11 ˆ
A12 ˆ
A13
ˆ
A21 ˆ
A22 ˆ
A23
ˆ
A31 ˆ
A32 ˆ
A33
x1
x2
x3
+
ˆ
B1
ˆ
B2
0
u(23)
with the property that ˆ
A33 is invertible and the reduced system obtained by solving for x3is
strangeness-free (of differentiation index at most one) when setting u= 0.
The described regularization procedure holds for general DAEs but it does not reflect
an available port-Hamiltonian structure. We will now modify this approach for nonlinear
systems with a pHDAE structure, which (based on the linear time-varying formulation) we
define as follows.
Definition 14 Consider a general DAE model in the form (16) and a Hamiltonian H(x)with
the property that for a given input u(t)and associated trajectory x(t)the Hessian Yloc(t) =
Hxx(x(t)) can be expressed locally as Eloc(t)TQloc(t), where Eloc(t) = F˙x(t),Fx(t) = (Jloc(t)−
13
Rloc(t))Qloc(t)−Eloc(t)Kloc(t),Fu(t) = Bloc(t)−Ploc(t),Gx(t) = (Bloc(t) + Ploc(t))TQloc(t),
Gu(t) = Sloc(t) + Nloc(t), with Eloc(t), Aloc(t), Qloc(t), Rloc(t) = RT
loc(t), Kloc(t)∈C0(I,Rn,n),
Bloc(t), Ploc(t)∈C0(I,Rn,m),Sloc(t) = ST
loc(t), Nloc(t) = −NT
loc(t)∈C0(I,Rm,m). Then the
system is called pHDAE system if the following properties are satisfied:
i) the (local) differential-algebraic operator
Lloc := QT
loc(t)Eloc(t)d
dt −QT
loc(t)(Jloc(t)Qloc(t)−Eloc(t)Kloc(t)) (24)
is skew-adjoint,
ii) locally there exists a constant matrix H0such that Hessian Hxx(t)− H0is positive
semidefinite;
ii) locally
Wloc(t) = "Qloc(t)TRloc(t)Qloc(t)Qloc(t)TPloc(t)
PT
loc(t)Qloc(t)Sloc(t)#≥0 for all t∈I.(25)
Clearly standard pH systems of the form (2) and linear time-varying systems as in Definition 4
directly fit in this framework. The same holds for multibody systems as in Example 7 with
the derivative of the constraint
˜
M(Θ)¨
Θ + ˜
C(Θ,˙
Θ) + ˜
G(Θ) = ˜u+ ΨTλ,
Ψ(Θ) ˙
Θ=0.(26)
Remark 15 There is a lot of choice in the local matrices Qloc and Eloc when factoring
the Hessian. In some cases we can just choose Qloc to be the identity (see Remark 8), so
that Eloc =ET
loc defines the Hessian. In other cases one chooses the block-diagonal matrix
Eloc = diag(Id,0) and obtains a semi-explicit formulation of the pHDAE. However, in general,
this freedom should be chosen to make the system robust to perturbations for simulation and
control methods.
There are multiple reasons why constraints may arise in pH systems. A typical example arises
as a limiting situation in a singularly perturbed problem which has pH structure. Typical
examples are mechanical multibody systems where small masses are ignored.
Example 16 Finite element modeling of the acoustic field in the interior of a car, see e. g.
[32], leads (after several simplifications) to a large scale constant coefficient differential-
algebraic equation system of the form
M¨p+D˙p+Kp =B1u,
where pis the coefficient vector associated with the pressure in the air and the displacements
of the structure, B1uis an external force, Mis a symmetric positive semidefinite mass matrix,
Dis a symmetric positive semidefinite matrix, and Kis a symmetric positive definite stiffness
matrix. Here Mis only semidefinite since small masses were set to zero, so Mis a perturbation
14
of a positive definite matrix. Performing a first order formulation we obtain the state equation
of a pHDAE system E˙z= (J−R)Qz +Bu, where
E:= M0
0I, J := 0−I
I0, R := D
0, z := ˙p
p,
Q:= I0
0K, B := B1
0, P := 0,
and the Hamiltonian is
H=1
2(zTETQz) = 1
2( ˙pTM˙p+pTKp).
Note that this model is nonlinear originally, but the simplifications carried out in the modeling
process, e. g. linearization and omitting nonlinear terms with small coefficients leads to a
linear model.
The other class of examples are systems such as as Example 7, where the dynamics is con-
strained to a manifold. If the system like in this example has hidden constraints, then the
formulation as pHDAE system is not unique because different formulations of the equations
and the constraints can be be made. We will come back to this question in Section 6.
As mentioned in the introduction, it is sometimes claimed that port-Hamiltonian DAEs
are of differentiation-index at most one (i. e., satisfy Hypothesis 13 with µ= 0). If this would
be the case then in the derivative array with µ= 0 the matrix F˙xlocally has constant rank
dand if ZT
2is a maximal rank matrix such that (locally) ZT
2F˙x= 0 and Z1is such that it
completes Z2to an invertible matrix Z= [Z1, Z2], then the matrix ¯
E:= ZT
1F˙x
ZT
2Fxis locally
invertible.
Let us check this for some of the examples. In Example 5 we have ZT
2=0 0 Iand
obtain
¯
E=
GcCGT
c0 0
0L0
−GT
v0 0
which is clearly not invertible, except if the last row and column is empty. The same matrix
Z2can be used in Example 6 and yields
¯
E=
M10 0
0M20
0N0
which is also not invertible except if the last row and column is empty. Actually due to
the special structure it can be shown that both systems have µ= 1, i. e., differentiation-
index two, when the input is chosen to be 0. The analysis of Example 7 with the original
constraint Ψ(Θ) has µ= 2 (differentiation-index three) and the formulation as pHDAE is not
straightforward, but using as constraint its derivative yields µ= 1 (differentiation-index two)
if ˜
G˜
GTis invertible. This replacement corresponds to an index reduction. How to carry out
such a regularization for pHDAEs will be discussed in Section 6. But let us first (in the next
section) analyze in detail the case of differentiation-index one pHDAEs.
15
5 PHDAEs of differentiation-index at most one
In this section we characterize pHDAE systems of differentiation-index at most one (µ= 0)
and we first study linear time-varying pHDAE systems. In this case Hypothesis 13 implies
that the matrix function E(t) has constant rank. Then, see e. g. Theorem 3.9 in [27], there
exist pointwise orthogonal matrix functions ˜
Uand Vsuch that
˜
UTEV =E11 0
0 0 ,
where E11 is pointwise invertible. Because QTEis real symmetric, in
˜
UTQV =Q11 Q12
Q21 Q22 ,
one has QT
11E11 =ET
11Q11 and also Q12 = 0. Partition in the same way
˜
UTJ˜
U=J11 J12
J21 J22 ,˜
UTR˜
U=R11 R12
RT
12 R22 ˜x=VTx=x1
x2,
˜
UTB P =B1P1
B2P2,˜
K=K+VT˙
V=K11 K12
K21 K22 ,
˜
UT(J−R)˜
U=J11 J12
J21 J22 −R11 R12
RT
12 R22 =: L11 L12
L21 L22 ,
so that the transformed pHDAE system has the form
E11 0
0 0 ˙x1
˙x2=L11 L12
L21 L22 Q11 0
Q21 Q22 −E11K11 E11K12
0 0 x1
x2
+B1−P1
B2−P2u,
y=(B1+P1)T(B2+P2)TQ11 0
Q21 Q22 x1
x2+ (S+N)u.
Since the system has differentiation-index at most one, the block L22Q22 either does not
occur (in this case we have an implicitly defined standard pH system) or it must be pointwise
invertible, see [27], i. e., both L22 and Q22 are pointwise invertible. Let U=˜
UT, where
T:= I0
T21 I, T21 =−L−T
22 (L12 −E11K12Q−1
22 )T.
Then a transformation of the original pHDAE with Uand Vyields the pHDAE system
˜
E˙x1
˙x2= [( ˜
J−˜
R)˜
Q−˜
E˜
K)] x1
x2+ ( ˜
B−˜
P)u,
y= ( ˜
B+˜
P)T˜
Qx1
x2+ ( ˜
S+˜
N)u, (27)
16
where
˜
E=TTE11 0
0 0 =E11 0
0 0 ,˜
K=K, ˜
S=S, ˜
N=N,
˜
Q=T−1Q11 0
Q21 Q22 =Q11 0
Q21 −T21Q11 Q22 =: Q11 0
˜
Q21 Q22 ,
˜
J=TTJ11 J12
J21 J22 T=J11 +TT
21J21 +J12T21 +TT
21J22T21 J12 +TT
21J22
J21 +J22T21 J22
=: ˜
J11 ˜
J12
˜
J21 J22 ,
˜
R=TTR11 R12
RT
12 R22 T=R11 +TT
21RT
12 +J12R12 +TT
21R22T21 R12 +TT
21R22
RT
12 +R22T21 R22
=: ˜
R11 ˜
R12
˜
RT
12 R22 ,
˜
B=TTB1
B2=B1+TT
21B2
B2=: ˜
B1
B2,˜
P=TTP=P1+TT
21P2
P2=: ˜
P1
P2.
Following Theorem 10 this transformation will not change the Hamiltonian, and based on the
construction of T,
(˜
J−˜
R)˜
Q−˜
E˜
K=T((J−R)Q−EK)
=(˜
J11 −˜
R11)Q11 −E11(K11 −K12Q−1
22 ˜
Q21) 0
(˜
J21 −˜
RT
12)Q11 + (J22 −R22)˜
Q21 (J22 −R22)Q22 .
Note that these transformations should not be performed in a numerical integration or control
design technique, since the inversion of the matrices Q22 and L22 may be highly ill-conditioned.
However, from an analytic point of view have the following theorem.
Theorem 17 Suppose that the pHDAE system (5) is of differentiation-index at most one
(i. e. satisfies Hypothesis 13 with µ= 0) for ν= 0, and that E(t)has constant rank. Let
U, V and ˜
E,˜
Q,˜
J,˜
R,˜
B,˜
Pbe as in (27) and let VTx=xT
1xT
2T. Then system (5) can
be reduced to the implicit pHDAE system (for the state x1)
ˆ
E˙x1= [( ˆ
J−ˆ
R)ˆ
Q−ˆ
Eˆ
K]x1+ ( ˆ
B−ˆ
P)u
y= ( ˆ
B+ˆ
P)Tˆ
Qx1+ ( ˆ
S+ˆ
N)u(28)
with Hamiltonian ˆ
H(x1(t)) = 1
2xT
1ˆ
QTˆ
Ex1=H(x), and coefficients
ˆ
E=E11,ˆ
Q=Q11,ˆ
J=˜
J11,ˆ
R=˜
R11,ˆ
K=K11 −K12Q−1
22 ˜
Q21,
ˆ
B=˜
B1−1
2(˜
JT
21 −˜
R12)L−T
22 (B2+P2),ˆ
P=˜
P1−1
2(˜
JT
21 −˜
R12)L−T
22 (B2+P2),
ˆ
S=S−1
2[(B2+P2)TL−1
22 (B2−P2)+(B2−P2)TL−T
22 (B2+P2)],
ˆ
N=N−1
2[(B2+P2)TL−1
22 (B2−P2)−(B2−P2)TL−T
22 (B2+P2)],
together with the explicit algebraic constraint
L22Q22x2=−[( ˜
J21 −˜
RT
12)Q11 +L22 ˜
Q21]x1−(B2−P2)u, (29)
17
for the state x2, which also gives a consistency constraint for the initial condition
L22(t0)Q22(t0)x2(t0) = −[(L21(t0))Q11(t0) + L22(t0)˜
Q21(t0)]x1(t0)−(B2(t0)−P2(t0))u(t0).
Proof. Equations (28) and (29) follow directly from (27). The output equation is obtained
directly by substituting (29). It remains to prove that (28) is port-Hamiltonian. From
d
dt ˜
QT˜
E=−˜
QT(˜
J+˜
JT)˜
Q+˜
QT˜
E˜
K+˜
KT˜
ET˜
Q,
and since Q22 is invertible, one obtains
0 = J22 +JT
22,
0 = −QT
11(˜
J12 +˜
JT
21)−˜
QT
21(J22 +JT
22) + QT
11E11K12Q−1
22 ,
d
dtQT
11E11 =−QT
11(˜
J11 +˜
JT
11)Q11 −QT
11(˜
J12 +˜
JT
21)˜
Q21
−˜
QT
21(˜
J12 +˜
JT
21)TQ11 −˜
QT
21(J22 +JT
22)˜
Q21 +QT
11E11K11 +KT
11ET
11Q11,
which leads to
J22 +JT
22 = 0,(30)
QT
11(˜
J12 +˜
JT
21) = QT
11E11K12Q−1
22 ,(31)
and
d
dtQT
11E11 =−QT
11(˜
J11 +˜
JT
11)Q11 +QT
11E11(K11 −K12Q−1
22 ˜
Q21)+(K11 −K12Q−1
22 ˜
Q21)TET
11Q11,
and this last equation is just
d
dt ˆ
QTˆ
E=−ˆ
QT(ˆ
J+ˆ
JT)ˆ
Q+ˆ
QTˆ
Eˆ
K+ˆ
KTˆ
ETˆ
Q.
Since ˆ
QTˆ
E=ˆ
ETˆ
Q, the operator ˆ
QTˆ
Ed
dt −ˆ
QT(ˆ
Jˆ
Q−ˆ
Eˆ
K) is skew-adjoint.
The invariance of the Hamiltonian follows directly, since
ˆ
H(x1) = 1
2xT
1ˆ
QTˆ
Ex1=1
2x1
x2T
˜
QT˜
Ex1
x2=1
2xTQTEx =H(x).
It remains to prove the dissipation inequality. We have that
d
dt ˆ
H(x1) = d
dtH(x) = yTu−x
uT
Wx
u=yTu−
x1
x2
u
T
˜
W
x1
x2
u
where
˜
W=
Q11 0 0
˜
Q21 Q22 0
0 0 I
T
˜
R11 ˜
R12 ˜
P1
˜
RT
12 R22 P2
˜
PT
1PT
2S
Q11 0 0
˜
Q21 Q22 0
0 0 I
.
18
Eliminating x2by using (29), we obtain
x1
x2
u
T
˜
W
x1
x2
u
=x1
uT
ˆ
Wx1
u,
where ˆ
W=XT˜
WX with
X=
I0
−Q−1
22 (L−1
22 (˜
J21 −˜
RT
12)Q11 +˜
Q21)−Q−1
22 L−1
22 (B2−P2)
0I
,
Note that L22 =J22 −R22 and by (30) we have J22 =−JT
22, and thus R22 =−1
2(L22 +LT
22).
Also, from (31) and the formulas of ˜
J21,˜
R12, T21, it follows that QT
11(˜
JT
21 +˜
R12) = 0. Then,
by straightforward calculations we obtain
ˆ
W=QT
11 ˆ
RQ11 QT
11 ˆ
P
ˆ
PTQ11 ˆ
S.
Hence
d
dt ˆ
H=yTu−x1
uT
ˆ
Wx1
u.
Since Wis symmetric positive semidefinite, so is ˆ
W, and hence the reduced system in x1is
still port-Hamiltonian with Hamiltonian ˆ
H(x1).
Note that for the numerical integration or in the control context, as for general DAEs, it
is sufficient to carry out the transformation with ˜
Upointwise from the left and the insertion
of I=˜
U˜
UTbefore Q. In this way a differentiation of a computed transformation matrix can
be avoided and the pHDAE structure is preserved nonetheless.
Remark 18 For nonlinear pHDAE systems with differentiation index at most one (µ= 0),
the corresponding local result follows directly via the implicit function theorem and applica-
tion of Theorem 17 to the linearization as in Definition 14.
6 Regularization of higher index pHDAE systems
In this section we discuss how to modify the regularization procedure discussed for general
DAEs in Section 4 to preserve the pHDAE structure. Let us first consider the linear time-
varying case (5) and set L=J−R. Suppose that the state equation with u= 0 already
satisfies Hypothesis 13, i. e., as discussed in Section 4, no reinterpretation of variables or initial
feedbacks are necessary. It has been shown in [7] that the extra constraint equations (hidden
constraints) that arise from derivatives are uncontrollable, because otherwise the index re-
duction could have been done via feedback. This means that these extra constraint equations
are of the form ˆ
A3x= 0 which corresponds to ˆ
F3(t, x) = 0 in the nonlinear case. We add just
these constraint equations to our original pHDAE and obtain an overdetermined strangeness-
free system. Note again that under our assumptions the explicit algebraic constraints are
included in the first two equations in (20), resp. (21).
Let us make the weak assumption that E(t) has constant rank. This is a restriction
that however holds in all examples that we have encountered so far, and it can be removed
19
by considering the system in a piecewise fashion, see [27]. Then there exist real orthogonal
matrix functions U, V1∈C1(I,Rn,n) such that
UT
1EV1=˜
E11 0
0 0
with pointwise invertible ˜
E11. Perform a transformation of the pHDAE (5) as in Theorem 10
and also form ˆ
A3V=ˆ
A31 ˆ
A32 partitioned accordingly. By the property that ˆ
A3contains
all the high index constraints it follows that ˆ
A23 has full row rank for all t∈I, and hence
there exists a real orthogonal matrix function V2such that ˆ
A32V2=0A33 with A33
pointwise invertible. Performing a change of variables of the pHDAE with
V:= V1I0
0V2
I0 0
0I0
−ˆ
A31A−1
33 0I
we obtain a pHDAE of the form
˜
E11 0 0
0 0 0
0 0 0
˙x1
˙x2
˙x3
=˜
L
˜
Q11 ˜
Q12 ˜
Q13
˜
Q21 ˜
Q22 ˜
Q23
˜
Q31 ˜
Q32 ˜
Q33
x1
x2
x3
−
˜
E11 0 0
0 0 0
0 0 0
K11 K12 K13
K21 K22 K23
K31 K32 K33
x1
x2
x3
+
˜
B1−˜
P1
˜
B2−˜
P2
˜
B3−˜
P3
u, (32)
y=(˜
B1+˜
P1)T(˜
B2+˜
P2)T(˜
B3+˜
P3)T
˜
Q11 ˜
Q12 ˜
Q13
˜
Q21 ˜
Q22 ˜
Q23
˜
Q31 ˜
Q32 ˜
Q33
x1
x2
x3
+ (S+N)u,
where ˜
K= (VTKV +˙
V), ˜
L=LV , together with the constraint 0 = A33x3, i. e. x3= 0.
Assuming further that the matrix function
˜
Q11 ˜
Q12
˜
Q21 ˜
Q22
˜
Q31 ˜
Q32
has constant rank, there exists a pointwise real orthogonal matrix function U2such that
UT
2
˜
Q11 ˜
Q12 ˜
Q13
˜
Q21 ˜
Q22 ˜
Q23
˜
Q31 ˜
Q32 ˜
Q33
=
Q11 Q12 Q13
Q21 Q22 Q23
0 0 Q33
20
Transforming the pHDAE (34) with UT
2we get a pHDAE of the form
E11 0 0
E21 0 0
E31 0 0
˙x1
˙x2
˙x3
=
L11 L12 L13
L21 L22 L23
L31 L32 L33
Q11 Q12 Q13
Q21 Q22 Q23
0 0 Q33
x1
x2
x3
−
˜
E11 0 0
˜
E21 0 0
˜
E31 0 0
K11 K12 K13
K21 K22 K23
K31 K32 K33
x1
x2
x3
+
˜
B1−˜
P1
˜
B2−˜
P2
˜
B3−˜
P3
u, (33)
y=(˜
B1+˜
P1)T(˜
B2+˜
P2)T(˜
B3+˜
P3)T
˜
Q11 ˜
Q12 ˜
Q13
˜
Q21 ˜
Q22 ˜
Q23
˜
Q31 ˜
Q32 ˜
Q33
x1
x2
x3
+ (S+N)u,
together with the constraint 0 = x3.
By Theorem 10, system (34) is still a pHDAE system, and the solution of the overdeter-
mined system (34) together with x3= 0 is the same as that of (34) and the Hamiltonian is
unchanged. Since the resulting system is still port-Hamiltonian, using that x3= 0, we have
that the subsystem given by the first two block rows together with output equation is an
index at most one phDAE which has the form
E11 0
E21 0 ˙x1
˙x2=L11 L12
L21 L22 Q11 Q12
Q21 Q22 x1
x2
−E11 0
E21 0K11 K12
K21 K22 x1
x2+˜
B1−˜
P1
˜
B2−˜
P2u, (34)
y=(˜
B1+˜
P1)T(˜
B2+˜
P2)T˜
Q11 ˜
Q12
˜
Q21 ˜
Q22 x1
x2
+ (S+N)u,
To this system we can apply the results of the previous section and obtain that the system
can be further reduced to an implicit standard pH system.
Example 19 Consider again the semidiscretized Example 6. It has been shown in [15] that
for a (permuted) singular value decomposition (SVD) of NT
N>=U>
N0
ΣVN,
with real orthogonal matrices UN, VNand a nonsingular diagonal matrix Σ ∈Rn3,n3. Scal-
ing the second row of (10) with UNand setting x2=VNx>
2,2x>
2,3>, as well as x0
2=
VNhx0
2,2
>x0
2,3
>i>we obtain a transformed system
M10 0 0
0M2,2M2,30
0M>
2,3M3,30
0 0 0 0
d
dt
x1
x2,2
x2,3
x3
+
0G1,2G1,30
−G>
1,2D2,2D2,30
−G>
1,3D>
2,3D3,3−Σ
0 0 Σ 0
x1
x2,2
x2,3
x3
=
0
B2,2
B3,2
0
u. (35)
21
It follows immediately that x2,3= 0, which is the uncontrollable (index two) constraint in
the DAE that in particular the initial condition x0
2,3has to satisfy. The vectors x1, x2,2are
solutions of the implicit ordinary pH system
M10
0M2,2d
dt x1
x2,2+0G1,2
−G>
1,2D2,2x1
x2,2=0
B2,2u, (36)
with initial conditions x1(0) = x0
1,x2,2(0) = x0
2,2, so they are well-defined continuously differ-
entiable functions for any piecewise continuous uand any choice of the initial conditions.
Finally we get the component x3(the Lagrange multiplier) via
x3= Σ−1(M>
2,3
d
dtx2,2−G>
1,3x1+D>
2,3x2,2−B3,2u),(37)
and this is the implicit index one constraint in the DAE. Since both type of (the explicit
and the hidden) constraints have to be satisfied for the initial condition, it means that the
transformed initial condition also has to satisfy the consistency condition
x3(0) = Σ−1(M>
2,3
d
dtx2,2(0) −G>
1,3x1(0) + D>
2,3x2,2(0) −B3,2u(0)) (38)
Condition (38) leads to a relationship between the input uand the state at t= 0, which
is a constraint that has to be satisfied to have a classical solution. Furthermore, we see
immediately that to obtain a continuous x3the function B3,2uhas to be continuous and
uhas to be such that B3,2uleads to a continuous M>
2,3d
dt x2,2. The implicit ordinary pH
system (36) describes the dynamics of the system, while the other two equations describe the
constraints.
Remark 20 For nonlinear pHDAE systems satisfying Hypothesis 13 with µ > 0, the corre-
sponding local result follows directly via linearization and the implicit function theorem.
Conclusion
A new definition of port-Hamiltonian descriptor systems has been derived. It has been shown
that this formulation is valid also for DAEs of differentiation-index larger than one, and it has
been demonstrated that under some (local) constant rank assumption any such pHDAE can
be reformulated as an implicitly defined standard PH system plus an algebraic constraint that
describes the manifold where the dynamics of the system takes place and that also describes
the consistent initial conditions. As for standard DAEs the reformulated system is well suited
for numerical integration and control, since all constraints are available.
Acknowledgments
We acknowledge many interesting discussions with Robert Altmann and Philipp Schulze from
TU Berlin and Arjan Van der Schaft from RU Groningen. The first author has been supported
by Einstein Foundation Berlin, through an Einstein Visiting Fellowship. The second author
has been supported by Deutsche Forschungsgemeinschaft for Research support via Project
A02 in CRC 1029 TurbIn and by Einstein Foundation Berlin within the Einstein Center
ECMath.
22
References
[1] J. Bals, G. Hofer, A. Pfeiffer, and C. Schallert. Virtual iron bird – A multidisciplinary
modelling and simulation platform for new aircraft system architectures. In Deutscher
Luft- und Raumfahrtkongress, Friedrichshafen, Germany, 2005.
[2] C. Beattie and S. Gugercin. Structure-preserving model reduction for nonlinear port-
Hamiltonian systems. In 50th IEEE Conference on Decision and Control and European
Control Conference (CDC-ECC), 2011, pages 6564–6569. IEEE, 2011.
[3] A. Binder, V. Mehrmann, A. Miedlar, and P. Schulze. A Matlab toolbox for the regu-
larization of descriptor systems arising from generalized realization procedures. Preprint
24–2015, Institut f¨ur Mathematik, TU Berlin, 2015.
[4] P. C. Breedveld. Modeling and Simulation of Dynamic Systems using Bond Graphs,
pages 128–173. EOLSS Publishers Co. Ltd./UNESCO, Oxford, UK, 2008.
[5] K. E. Brenan, S. L. Campbell, and L. R. Petzold. Numerical Solution of Initial-Value
Problems in Differential Algebraic Equations. SIAM Publications, Philadelphia, PA, 2nd
edition, 1996.
[6] A. Bunse-Gerstner, R. Byers, V. Mehrmann, and N. K. Nichols. Feedback design for
regularizing descriptor systems. Linear Algebra Appl., 299:119–151, 1999.
[7] R. Byers, P. Kunkel, and V. Mehrmann. Regularization of linear descriptor systems with
variable coefficients. SIAM J. Control Optim., 35(1):117–133, 1997.
[8] C. I. Byrnes, A. Isidori, and J. C. Willems. Passivity, feedback equivalence, and the
global stabilization of minimum phase nonlinear systems. IEEE Trans. Autom. Control,
36:1228–1240, 1991.
[9] S. L. Campbell. A general form for solvable linear time varying singular systems of
differential equations. 18:1101–1115, 1987.
[10] S. L. Campbell. Linearization of DAEs along trajectories. Z. Angew. Math. Phys.,
46:70–84, 1995.
[11] S. L. Campbell, P. Kunkel, and V. Mehrmann. Regularization of linear and nonlinear
descriptor systems. In L. T. Biegler, S. L. Campbell, and V. Mehrmann, editors, Con-
trol and Optimization with Differential-Algebraic Constraints, Advances in Control and
Design, chapter 2, pages 17–36. 2012.
[12] F. Couenne, C. Jallut, B. M. Maschke, M. Tayakout, and P. C. Breedveld. Bond graph
for dynamic modelling in chemical engineering. Chemical engineering and processing,
Elsevier, Amsterdam, 47:1994–2003, 2008.
[13] L. Dai. Singular Control Systems, volume 118 of Lecture Notes in Control and Inform.
Sci. Springer-Verlag, Berlin/Heidelberg, 1989.
[14] H. Egger and T. Kugler. Damped wave systems on networks: Exponential stability and
uniform approximations. Technical report, 2016. http://arxiv.org/abs/1605.03066.
23
[15] H. Egger, T. Kugler, B. Liljegren-Sailer, N. Marheineke, and V. Mehrmann. Structure
preserving model reduction for damped wave propagation in transport networks. Preprint
06-2017, Institut f¨ur Mathematik, TU Berlin, 2017. http://arxiv.org/abs/1704.03206.
[16] E. Eich-Soellner and C. F¨uhrer. Numerical Methods in Multibody Dynamics. Teubner,
Stuttgart, 1998.
[17] R. W Freund. The sprim algorithm for structure-preserving order reduction of general
rcl circuits. In Model reduction for circuit simulation, pages 25–52. Springer, 2011.
[18] G. Golo, A. J. van der Schaft, P. C. Breedveld, and B. M. Maschke. Hamiltonian
formulation of bond graphs. In A. Rantzer R. Johansson, editor, Nonlinear and Hybrid
Systems in Automotive Control, pages 351–372. Springer, Heidelberg, 2003.
[19] S. Gugercin, R. V. Polyuga, C. Beattie, and A. J. van der Schaft. Structure-preserving
tangential interpolation for model reduction of port-Hamiltonian systems. Automatica,
48:1963–1974, 2012.
[20] M. Hiller and K. Hirsch. Multibody system dynamics and mechatronics.
Z. Angew. Math. Mech., 86(2):87–109, 2006.
[21] D. Hinrichsen and A. J. Pritchard. Mathematical Systems Theory I. Modelling, State
Space Analysis, Stability and Robustness. Springer-Verlag, New York, NY, 2005.
[22] M. Hou. A three–link planar manipulator model. Sicherheitstechnische Regelungs- und
Meßtechnik, Bergische Universit¨at–GH Wuppertal, Germany, May 1994.
[23] M. Hou and P. C. M¨uller. LQ and tracking control of descriptor systems with application
to constrained manipulator. Technical report, Sicherheitstechnische Regelungs- und Meß-
technik, Universit¨at Wuppertal, Gauß-Straße 20, D-5600 Wuppertal 1, Germany, 1994.
[24] B. Jacob and H. Zwart. Linear port-Hamiltonian systems on infinite-dimensional spaces.
Operator Theory: Advances and Applications, 223. Birkh¨auser/Springer Basel AG, Basel
CH, 2012.
[25] C. Kleijn. 20-sim 4C 2.1 Reference manual. Controlab Products B.V., 2013.
[26] P. Kunkel and V. Mehrmann. Analysis of over- and underdetermined nonlinear
differential-algebraic systems with application to nonlinear control problems. Math. Con-
trol Signals Systems, 14:233–256, 2001.
[27] P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and Numerical
Solution. EMS Publishing House, Z¨urich, Switzerland, 2006.
[28] P. Kunkel, V. Mehrmann, and W. Rath. Analysis and numerical solution of control
problems in descriptor form. Math. Control Signals Systems, 14:29–61, 2001.
[29] P. Kunkel, V. Mehrmann, and L. Scholz. Self-adjoint differential-algebraic equations.
Math. Control Signals Syst., 26:47–76, 2014.
[30] R. Lamour, R. M¨arz, and C. Tischendorf. Differential-algebraic equations: a projector
based analysis. Springer Science & Business Media, 2013.
24
[31] B. M. Maschke, A. J. van der Schaft, and P. C. Breedveld. An intrinsic Hamiltonian
formulation of network dynamics: non-standard poisson structures and gyrators. J.
Franklin Inst., 329:923–966, 1992.
[32] V. Mehrmann and C. Schr¨oder. Nonlinear eigenvalue and frequency response problems
in industrial practice. J. Math. Ind., 1:18, 2011.
[33] R. Ortega, A. J. van der Schaft, Y. Mareels, and B. M. Maschke. Putting energy back
in control. Control Syst. Mag., 21:18–33, 2001.
[34] J. W. Polderman and J. C. Willems. Introduction to Mathematical Systems Theory: A
Behavioural Approach. Springer-Verlag, New York, NY, 1998.
[35] R. V. Polyuga and A. J. van der Schaft. Structure preserving model reduction of port-
Hamiltonian systems by moment matching at infinity. Automatica, 46:665–672, 2010.
[36] A. J. van der Schaft. Port-Hamiltonian systems: an introductory survey. In J. L. Verona
M. Sanz-Sole and J. Verdura, editors, Proc. of the International Congress of Mathemati-
cians, vol. III, Invited Lectures, pages 1339–1365, Madrid, Spain.
[37] A. J. van der Schaft. Port-Hamiltonian systems: network modeling and control of non-
linear physical systems. In Advanced Dynamics and Control of Structures and Machines,
CISM Courses and Lectures, Vol. 444. Springer Verlag, New York, N.Y., 2004.
[38] A. J. van der Schaft. Port-Hamiltonian differential-algebraic systems. In Surveys in
Differential-Algebraic Equations I, pages 173–226. Springer-Verlag, 2013.
[39] A. J. van der Schaft and B. M. Maschke. Hamiltonian formulation of distributed-
parameter systems with boundary energy flow. J. Geom. Phys., 42:166–194, 2002.
[40] W. Schiehlen. Advanced multibody system dynamics. Kluwer Academic Publishers,
Stuttgart, Germany, 1993.
[41] K. Schlacher and A. Kugi. Automatic control of mechatronic systems.
Int. J. Appl. Math. Comput. Sci., 11(1):131–164, 2001.
[42] W. Trautenberg. Simpack 8.9. Manual, INTEC GmbH, Angelrieder Feld 13, 82234
Wessling.
25
