Structurepreserving discretization for portHamiltonian descriptor systems
Abstract
We extend the modeling framework of portHamiltonian descriptor systems to include under and overdetermined systems and arbitrary differentiable Hamiltonian functions. This structure is associated with a Dirac structure that encloses its energy balance properties. In particular, portHamiltonian 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 timediscretization 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 portHamiltonian (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 GaussLegendre 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 differentialalgebraic 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 timevarying 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 infinitedimensional states and time and spacevarying 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 structurepreserving time discretization to portHamiltonian 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 GaussLegendre collocation.
2 PortHamiltonian descriptor systems
2.1 Formulation
Definition 1 (pHDAE).
Let us consider a time interval , a state space and the extended state space . A portHamiltonian descriptor system (or pHDAE) is a system of differential (algebraic) equations of the form
(1) 
associated with an Hamiltonian function , where is the state, are the input and output, is the flow matrix, are the timeflow and effort functions, are the structure and dissipation matrices, are the port matrices and are the feedthrough matrices. Furthermore, the following properties must hold:

The extended structure and dissipation matrices , defined as
(2) satisfy and pointwise.

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 partialdifferentialalgebraic equations (PDAEs). For simplicity, the results and examples presented in this paper will be limited to the finitedimensional 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 nonnegative [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
PortHamiltonian 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 portHamiltonian 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)
(3) 
holds along any solution , for any input . In particular, the dissipation inequality
(4) 
holds.
Proof.
Note that, with no input (), if the Hamiltonian is locally positivedefinite 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 inputoutput DAE
(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 , premultiplying 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.

By substitution, we get

By substitution, we get
and
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:
(6) 
Note that condition 2 becomes simply , where the tilde denotes the quantities in the autonomous system.
2.2.4 Structurepreserving interconnection
Let us consider two autonomous pHDAEs of the form
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
with Hamiltonian , where we have introduced new state variables that copy , and we aggregated , , , , and , where is the permutation matrix
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 energypreserving (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
PortHamiltonian 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
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
where is a skewsymmetric operator. Then is a Dirac structure.
Proof.
For generic and ,
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
is a Dirac structure over . Furthermore, the system of equations
(7)  
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
The condition can be written as
that together with the conditions (7) is equivalent to
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. ∎
3 Time discretization
Let us consider a finitedimensional pHDAE of the form (1). Under some regularity assumptions [5], many classes of RungeKutta methods can be applied to compute timediscretizations of differentialalgebraic equations. An important class of RungeKutta 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 differentialalgebraic equation in collocation points with for .
Let denote the th Lagrange interpolation polynomial with respect to the nodes , i.e.,
Then we can write
for certain that we have to compute, and also
where and . In particular, the constants for are the coefficients of the Butcher diagram of the associated RungeKutta method [3].
3.2 Dirac structure associated to discretization
Consider for simplicity the autonomous case. The nonautonomous 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 timediscretization as , i.e.,
with and . Taking , together with
(8a)  
(8b)  
(8c)  
(8d)  
(8e)  
(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
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
In the collocation points, the PBE is satisfied:
for . In particular, if we apply the quadrature rule associated to the collocation method, we get
where is the degree of exactness of the quadrature rule. We observe that, for the same reason,
(9a)  
(9b) 
so
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
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 GaussLegendre collocation methods. In particular, if we apply such a method, the integration of will be exact, i.e.
In particular, since we always have for GaussLegendre collocation, the dissipation term is always nonpositive, and we can write the discrete dissipation inequality
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:
(10) 
The energy in the system can be stored in the inductor and in the two capacitors, giving the Hamiltonian
(11) 
The system can then be written equivalently as an autonomous pHDAE of the form
(12a)  
(12b) 
with , , , , and
The power balance equation reads as
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 portHamiltonian framework: let