Structure-preserving discretization for port-Hamiltonian descriptor systems

# Structure-preserving discretization for port-Hamiltonian descriptor systems

Volker Mehrmann111 Institut für Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG. Email address: mehrmann@math.tu-berlin.de.  and Riccardo Morandin222 Institut für Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG. Email address: morandin@math.tu-berlin.de.
###### 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 properties. In particular, port-Hamiltonian systems are naturally passive and Lyapunov stable, because the Hamiltonian defines a Lyapunov function. The explicit representation of input 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 nonlinear 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 gradient 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 using 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 constraints, consistency requirements for initial conditions and additional regularity requirements. 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 dimensional 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 Hamiltonian 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 , a state space and the extended state space . 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 , where is the state, are the input and output, is the flow matrix, are the time-flow and effort functions, are the structure and dissipation matrices, are the port matrices and are the feed-through matrices. Furthermore, the following properties must hold:

1. The extended structure and dissipation matrices , defined as

 Γ:=[JB−BTN],W:=[RPPTS] (2)

satisfy and pointwise.

2. The gradient of the Hamiltonian satisfies and pointwise.

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.

###### 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)

 \textupd\textupdtH(t,x(t))=−[zu]TW[zu]+uTy (3)

holds along any solution , for any input . In particular, the dissipation inequality

 H(t2,x(t2))−H(t1,x(t1))≤∫t2t1u(τ)Ty(τ)\textupdτ (4)

holds.

###### Proof.

The pHDAE (1) can be written as

 [E˙x+r0]=(Γ−W)[zu]+[0y],

in particular

 \textupd\textupdtH(t,x(t)) =∂tH+∂xHT˙x=zT(E˙x+r)=[zu]T[E˙x+r0]= =[zu]T((Γ−W)[zu]+[0y])=−[zu]TW[zu]+uTy,

which is (3). By time integration, since , we immediately get (4). ∎

Note that, with no input (), if the Hamiltonian is locally positive-definite in an equilibrium point (up to shifting it by a constant), then is a Lyapunov function and the system is stable.

#### 2.2.2 Variable transformations

This class of pHDAEs is closed under many variable transformations:

###### Theorem 2.

Consider a pHDAE of the form (1). Let be a second state space, let , let be a local diffeomorphism (with respect to ) and let 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 , , , , , and , where we denote for any , and let . Then (5) is a pHDAE with Hamiltonian function , and to any solution of (5) corresponds a solution of (1) with . Furthermore, if is a global diffeomorphism for all , then the two systems are equivalent.

###### Proof.

The transformed DAE system is obtained from (1) by setting , pre-multiplying with and inserting in front of in the first equation. It is then clear that if is a solution of (5), then is a solution of the original system. If is a global diffeomorphism, then we can apply the inverse tranformation , where is to intended as , and get a solution for any solution 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∘φ)UUT(P∘φ)(P∘φ)TUS∘φ]==[U00I]T(W∘φ)[U00I]≥0.
2. By substitution, we get

 ∂~H∂~x(t,~x)=∂φ∂~xT(∂H∂x∘φ)=∂φ∂~xT(ETz∘φ)==∂φ∂~xT(ET∘φ)UU−1(z∘φ)=~ET~z

and

 ∂~H∂t(t,~x)=∂H∂t∘φ+(∂H∂x∘φ)T∂φ∂t==zTr∘φ+(zTE∘φ)∂φ∂t==(z∘φ)T(r∘φ+(E∘φ)∂φ∂t)==~zTUTU−T~r=~zT~r.

This concludes the proof. ∎

#### 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:

 [Er01][˙x˙t]=[J−R000][z0]+[B−P001][u1],[y0]=[B+P001]T[z0]+[S−N000][u1]. (6)

Note that condition 2 becomes simply , 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 , for , and assume that the aggregated input and output satisfy a linear interconnection relation for some . Then the aggregated system can be written as a pHDAE of the form

 ⎡⎢ ⎢ ⎢⎣E00000000000⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢⎣˙x˙^u˙^y⎤⎥ ⎥⎦ =⎡⎢ ⎢ ⎢ ⎢⎣Γ−W00Im−MT0−Im0M0−NTN0⎤⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣z^u^y0⎤⎥ ⎥ ⎥⎦+⎡⎢ ⎢ ⎢⎣00Im0⎤⎥ ⎥ ⎥⎦u, y =^y,

with Hamiltonian , where we have introduced new state variables that copy , and we aggregated , , , , and , where is the permutation matrix

 Π=⎡⎢ ⎢ ⎢ ⎢⎣Iℓ100000Iℓ200In100000In2⎤⎥ ⎥ ⎥ ⎥⎦.

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 defines a Dirac structure for ), 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 be a linear space and its dual space. Let be the bilinear form on defined as

 $$⟨$$$$⟨$$(f1,e1),(f2,e2)$$⟩$$$$⟩$$:=⟨e1|f2⟩+⟨e2|f1⟩,

where denotes the duality pairing. A Dirac structure on is then a linear subspace , such that .

In particular, in finite dimension, one only needs to prove and for all . If , then and are called flow and effort, respectively. In [11], the more general definition of modulated Dirac structure over is also presented, as a subbundle of , that denotes the Whitney sum between the tangent and cotangent bundles of . To introduce a Dirac structure for pHDAEs, we need to extend this further.

###### Definition 3.

Consider a state space and a vector bundle over with fibers . A Dirac structure over is a subbundle such that, for all , is a linear Dirac structure.

Note that modulated Dirac structures are a special case of our definition, where . To associate a Dirac structure to our pHDAE system, we first prove the following lemma:

###### Lemma 1.

Let be a vector subbundle with fibers defined by

 Dx={(f,e)∈Vx×V∗x:f+J(x)e=0},

where is a skew-symmetric operator. Then is a Dirac structure.

###### Proof.

For generic and ,

 $$⟨$$$$⟨$$(f,e),(f′,e′)$$⟩$$$$⟩$$ =⟨e|f′⟩+⟨e′|f⟩= =⟨e|f′⟩−⟨e′|Je⟩=⟨e|f′+Je′⟩.

We show that if and only if for all . On one hand, if , then holds for any . On the other hand, if , then and so such that , but with . ∎

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 for all , where is the storage flow fiber, is the port flow fiber and is the dissipation flow fiber. Let us partition and . Then the subbundle with

 Dx={(f,e)∈Vx×V∗x∣∣∣f+[Γ(x)Iℓ+m−Iℓ+m0]e=0}

is a Dirac structure over . Furthermore, the system of equations

 fs =−E(x)˙x, es =z(x), (7) fp =y, ep =u, ed =−W(x)fd, (f,e) ∈Dx

is equivalent to the original pHDAE, and represents the power balance equation.

###### Proof.

is 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].

The condition can 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 can be written as

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 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 , depending on time and possibly space, is given. Given a time interval of length , and a consistent initial condition at time , we approximate the solution of the pHDAE on with a polynomial , where denotes the vector space of polynomials of degree at most . The polynomial is chosen such that , and that it satisfies the differential-algebraic equation in collocation points with for .

Let denote the -th Lagrange interpolation polynomial with respect to the nodes , i.e.,

 ℓi(τ):=s∏j=1j≠iτ−γjγi−γj.

Then we can write

 ˙~x(t0+τh)=s∑i=1˙xiℓi(τ),~x(t0+τh)=x0+hs∑j=1˙xj∫τ0ℓj(σ)\textupdσ,

for certain that we have to compute, and also

 xi :=~x(ti)=x0+hs∑j=1αij˙xj, xf :=~x(tf)=x0+hs∑j=1βj˙xj,

where and . In particular, the constants for 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 be the Dirac structure associated to (6) (as in Theorem 3). We define the Dirac structure associated to the time-discretization as , i.e.,

 Dxi ={(fi,ei)∈Vxi×V∗xi∣∣∣fi+[Γ(ti,xi)Iℓ+m−Iℓ+m0]ei=0},

with and . Taking , together with

 xf =x0+hs∑i=1βi˙xi, (8a) xi =x0+hs∑j=1αij˙xj, (8b) fis =−E(xi)˙xi, (8c) eis =z(xi), (8d) eid =−W(xi)fid, (8e) ui =u(xi), (8f)

we get a system that is equivalent to applying the collocation method and computing the discrete input and output , for . Let us define the collocation flows, efforts, input and output in as

 ~fs(t0+hτ) =s∑i=1fisℓi(τ), ~es(t0+hτ) =s∑i=1eisℓi(τ), ~fd(t0+hτ) =s∑i=1fidℓi(τ), ~ed(t0+hτ) =s∑i=1eidℓi(τ), ~y(t0+hτ) =s∑i=1yiℓi(τ), ~u(t0+hτ) =s∑i=1uiℓi(τ).

Note that, by construction, in all collocation points . Let us consider the evolution of the Hamiltonian along the collocation polynomial . In particular, let : we can then write

 H(t)−H(t0)=∫tt0˙H(s)\textupds.

In the collocation points, the PBE is satisfied:

 ˙H(ti)=∇H(xi)T˙xi=z(xi)TE(xi)˙xi==−⟨eis|fis⟩=⟨eid|fid⟩+⟨yi|ui⟩,

for . In particular, if we apply the quadrature rule associated to the collocation method, we get

 H(tf)−H(t0)=hs∑j=1βj˙H(tj)+O(hp+1)==−hs∑j=1βj⟨ejs|fjs⟩+O(hp+1)==hs∑j=1βj⟨ejd|fjd⟩+hs∑j=1βj⟨yj|uj⟩+O(hp+1),

where is the degree of exactness of the quadrature rule. We observe that, for the same reason,

 hs∑j=1βj⟨ejd|fjd⟩ =∫tft0⟨~ed|~fd⟩+O(hp+1), (9a) hs∑j=1βj⟨yj|uj⟩ =∫tft0⟨~y|~u⟩+O(hp+1), (9b)

so

 H(tf)−H(t0)=∫tft0(⟨~ed(s)|~fd(s)⟩+⟨~y(s)|~u(s)⟩)\textupds+O(hp+1).

We note that, if , then equations (9a) and (9b) are exact. Furthermore, if for (and this is the case for many collocation methods), we can deduce that , thus the dissipation component retains its qualitative behaviour.

### 3.3 The quadratic Hamiltonian case

Let us consider the case where the Hamiltonian is a polynomial of degree (at most) 2, i.e. it can be written as

 H(x)=12xTQx+vTx+c,

for some , and . Since , we also have and . It is known that the maximum degree of exactness for quadrature rules with nodes is , 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 will be exact, i.e.

 H(tf)−H(t0)=hs∑j=1βj⟨ejd|fjd⟩+hs∑j=1βj⟨yj|uj⟩==∫tft0(⟨~ed(s)|~fd(s)⟩+⟨~y(s)|~u(s)⟩)\textupds.

In particular, since we always have for Gauss-Legendre collocation, the dissipation term is always non-positive, and we can write the discrete dissipation inequality

 H(tf)−H(t0)≤hs∑j=1βj⟨yj|uj⟩=∫tft0⟨~y(s)|~u(s)⟩\textupds.

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 are resistances, is an inductor, are capacitors and is a controlled voltage source. This can be interpreted as a basic representation of a DC generator (,), connected to a load () with a transmission line (the model given by ).

By means of Kirchhoff’s circuit laws, this system can be naturally written as the following DAE:

 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)=12LI2+12C1V21+12C2V22. (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 , , , , and

 J=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0−1100100−10−1000−10100000100⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,R=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣RL00000000000000000RG00000RR⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

The power balance equation reads as

 ˙H=−RLI2−RGI2G−RRI2R+IGEG;

in particular, if we shut down the generator (), the Hamiltonian will decrease and converge to a solution such that , that is, . The only state compatible with (10) that satisfies that condition is , so the system will always converge to that asymptotically stable point.

### 4.2 Controlling the circuit

Suppose now that represents a consumer, that requires a fixed amount of power to be delivered to them. We would like then to control the voltage of the generator , so that the state of the system will converge to . If we assume that a solution with exists, then we also get , , and .

This can be interpreted in the following way, exploiting the port-Hamiltonian framework: let

 x∗=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣I