Differential-algebraic systems with dissipative Hamiltonian structure

Different representations of dissipative Hamiltonian and port-Hamiltonian differential-algebraic equations (DAE) systems are presented and compared. Using global geometric and algebraic points of view, translations between the different representations are presented. Characterizations are also derived when a general DAE system can be transformed into one of these structured representations. Approaches for computing the structural information and the described transformations are derived that can be directly implemented as numerical methods. The results are demonstrated with a large number of examples.


Introduction
Since the original introduction of the energy based modeling concept of port-Hamiltonian (pH) systems in [28,35], see also [6,22,24,33,38,34,40], various novel definitions and formulations have been made to incorporate on the one hand systems defined on manifolds, see [14,41] and on the other hand systems with algebraic constraints defined in the form of Differential-Algebraic Equations (DAEs), see [4,31].Constraints in pH systems typically arise in electrical or transport networks, where Kirchhoff's laws constrain the models at network nodes, as balance equations in chemical engineering, or as holonomic or non-holonomic constraints in mechanical multibody systems.Furthermore, they arise in the interconnection of pH systems when the interface conditions are explicitly formulated and enforced via Lagrange multipliers, see [4,31,39,40,41,42] for a variety of examples, or the recent survey [32].
In this paper we recall different model representations of pHDAEs, respectively dissipative Hamiltonian DAEs (dHDAEs), and provide a systematic geometric and algebraic theory that also provides a 'translation' between the different representations.We also discuss approaches for explicitly computing the structural information that can be directly implemented as numerical methods.We mainly restrict ourselves to finite-dimensional linear time-invariant systems without inputs and outputs, but indicate in several places where extensions to general systems are possible.We also discuss only real systems although most of the results can also be formulated for complex problems in a similar way.
For general linear time-invariant homogeneous differential-algebraic equations where E, A ∈ R n,n (without additional geometric and algebraic structures), as well as for their extensions to linear time-varying and nonlinear systems, the theory is well understood [25].On the other hand, in contrast to general DAEs (1), the DAEs that arise in energy based modeling have extra structure and symmetries, and it is natural to identify and exploit this extra structure for purposes of analysis, simulation and control.
In this paper we will restrict our attention to regular systems, i.e., systems with det(λE − A) not identically zero, although many of the results are expected to extend to over-and underdetermined systems; cf. the Conclusions.The structural properties for regular systems are characterized via the (real) Weierstraß canonical form of the matrix pair (E, A), see e.g.[20].To determine this canonical form one computes nonsingular matrices U, W such that where J is in real Jordan canonical form and N is nilpotent in Jordan canonical form which is associated to the eigenvalue ∞ An important quantity that will be used throughout the paper is the size of the largest Jordan block in N (the index of nilpotency), which is called the (differentiation) index of the pair (E, A) as well as the associated DAE, where, by convention, ν = 0 if E is invertible.
The goal of this paper is to study in detail the relationships and differences between different classes of (extended) Hamiltonian DAEs; without making any a priori invertibility assumptions.Furthermore, we exploit at the same time a geometric point of view on Hamiltonian DAEs using Lagrange and Dirac structures (as well maximally monotone subspaces; cf.Section 4.2), and an algebraic point of view using (computationally feasible) condensed forms.From a DAE analysis point of view, similar to the result known already for Hamiltonian DAEs of the form (2) in [29], we also show that the index of an (extended) dissipative Hamiltonian DAEs can be at most two.Moreover, we will show that index two algebraic constraints may only arise from singularity of P , i.e., from the Lagrange structure, and thus are strictly linked to singular Hamiltonians.
Another important question that we will answer is the characterization when a general DAE is equivalent to a structured DAE the form ((2) or ( (10)).
The paper is organized as follows.
In Section 2 we discuss dissipative Hamiltonian DAEs and present several examples.The concept of extended Hamiltonian DAEs is discussed in Section 3. In Section 4 we will present the geometric theory of dissipative Hamiltonian DAEs, we introduce Dirac and Lagrange structures as well maximally monotone subspaces to incorporate dissipation.Section 5 presents different coordinate representations and their relation.To derive algebraic characterizations, in Section 6 we present different types of equivalence transformations and corresponding condensed forms for the different classes of (extended dissipative) Hamiltonian DAEs.In Section 7 it is analyzed when general linear DAEs can be represented as (extended dissipative) Hamiltonian systems.In all sections we present examples that illustrate the properties of the different representations.In several appendices we present proofs that can be implemented as numerically stable algorithms.

Dissipative Hamiltonian DAEs
A well established representation of DAE systems with symmetry structure is that of dissipative Hamiltonian DAEs (dHDAEs) which has been studied in detail in [4,29,30].These systems have the form where J = −J , R = R ∈ R , , and E, Q ∈ R ,n with E Q = Q E. In this representation the associated Hamilton function (Hamiltonian) is given by a quadratic form which (as an energy) is typically nonnegative for all z, but more general Hamiltonians also arise in practice.Note that in the ODE case, i.e. if E = I, then this reduces to the standard quadratic Hamiltonian H z (z) = 1 2 z Qz.

Remark 1
The equivalent formulations d dt (Ez) = E ż are both used in the literature and obviously lead to the same results in the linear constant coefficient case, but not anymore in the linear time-varying or nonlinear case.Actually, for the structured DAEs considered in this paper E is the product E = KP of two matrices K and P , and for generalizations it is appropriate to consider K d dt (P z).
Furthermore, the dissipation term R is typically positive semidefinite, denoted as R ≥ 0, but we will first discuss the case that R = 0 which we call Hamiltonian DAE (HDAE).

Remark 2
The definition of (dHDAEs) in (2) can be easily extended to systems with ports in which case it takes the form with see [4].We will call these systems dissipative pHDAE (dpHDAE) systems.
Remark 3 Note that by the symmetry of E Q the Hamiltonian H z (z) = 1 2 z E Qz can be also written as a function of Qz or Ez.Remark 4 In standard port-Hamiltonian modeling the matrix R in ( 2) is always assumed to be symmetric (as well as positive semi-definite).Alternatively, one could start from a general matrix R only satisfying R + R ≥ 0.
Then by splitting R into its symmetric and skew-symmetric part, one could add the skew-symmetric part to the skew-symmetric structure matrix J and continue as in (2) with the symmetric part 1 2 (R + R ) of R. Structured DAEs of the form (2) arise naturally in all physical domains.

Example 5
The subclass of RLC networks as considered in e.g.[4,12,19] has the form where the real positive definite diagonal matrices L, C, G describe inductances, capacitances, and conductances, respectively.The matrices D C , D L , D R , D S are the parts of the incidence matrix corresponding, respectively, to the capacitors, inductors, resistors (conductors), and current sources of the circuit graph, where D S is of full column rank.Furthermore V are the node potentials and I L , I S denote the currents through the inductors and sources, respectively.This system has the form (2), with E = E ≥ 0, Q equal to the identity matrix, and where J and −R are defined to be the skew-symmetric and symmetric part, respectively, of the matrix on the right hand side of (5).The Hamiltonian is given by H(V, L LI L and it does not involve the variables I S , which are in the kernel of E.
Example 6 Space discretization of the Stokes equation in fluid dynamics, see, e.g., [18], leads to a dissipative Hamiltonian system where A = A is a positive semidefinite discretization of the negative Laplace operator, B is a discretized gradient, and M = M is a positive definite mass matrix.The homogeneous system has the form (2), with The Hamiltonian is given by H = 1 2 v h M v h ; it does not involve the variables p h .
Example 7 Space discretization of the Euler equation describing the acoustic wave propagation in a gas pipeline network [15,16] where p h is the discretized pressure, q h is a discretized flux, and λ is a Lagrange multiplier that penalizes the violation of the conservation of mass and momentum at the pipeline nodes.
The homogeneous system has the form (2) with Q = I, and Hamiltonian ).The Hamiltonian does not involve the Lagrange multiplier λ.
Example 8 Consider a linear mechanical system (with q denoting the vector of position coordinates) M q + D q + W q = f , together with kinematic constraints G q = 0, see e.g.[17].The constraints give rise to constraint forces G λ, with λ a vector of Lagrange multipliers.This yields the dynamics M q + D q + W q = f + G λ, and the resulting system can be written in first order form as Here, E = E and Q = Q are commuting matrices, E, W , and R = R are positive semidefinite, so the homogeneous system is of the form (2). The Hamiltonian is given by 1 2 ( q M q + q W q) (kinetic plus potential energy) and does not involve the Lagrange multiplier λ.
As indicated in the presented examples, singularity of E, and thus the presence of algebraic constraints, implies that the Hamiltonian (the total stored energy) H z (z) = 1  2 z E Qz does not involve all of the variables contained in the vector z.In fact in some of the examples the variables that do not show up in the Hamiltonian are Lagrange multipliers.Conversely, singularity of E may arise as a limiting situation of otherwise regular Hamiltonians, as shown by the following simple example.
Example 9 Consider the model of standard mass-spring-damper system with model equation and Hamiltonian H(q, p) To compute the limit m → 0, we rewrite the system in coordinates q and v : with Hamiltonian H(q, v) = 1 2 kq 2 + 1 2 mv 2 .For m → 0 this converges to the which is of (differentiation) index one if d = 0, and of index two if d = 0.The limiting Hamiltonian H(q, v) = 1 2 kq 2 is not a function of v anymore.Alternatively, we can compute the limit for k → ∞.For this we rewrite the system in coordinates F := kq and p as , which is index two for any d, with Hamiltonian H(F, p) = 1 2m p 2 , which is not involving F .
We may also take the limits m → 0 and k → ∞ simultaneously.By rewriting the system in the variables F, v, with Hamiltonian H(F, p) = 1 2k F 2 + 1 2 mv 2 .This leads to the purely algebraic system with zero Hamiltonian, and having a single solution v = 0, F = 0, irrespective of the damper.

Extended Hamiltonian DAEs
Classically, Hamiltonian systems and also their extension to systems with inputs and outputs, called port-Hamiltonian systems, see e.g.[13,22,28,35,36,38,40], are defined by a Dirac structure, representing the powerconserving interconnection structure, an energy-dissipation relation, and a Hamilton function (Hamiltonian), capturing the total energy storage in the system.The Dirac structure formalizes the generalized junction structure known from port-based modeling theory.Importantly, the Dirac structure may entail linear constraints on the so-called effort variables, called effort constraints.Through the gradient vector of the Hamilton function, these effort constraints induce algebraic constraints on the state variables, see especially [31,39] for a discussion of general nonlinear port-Hamiltonian DAEs.In the special case of a linear autonomous system without energy dissipation and without inputs and outputs, and after choosing a basis such that X = R n , the equations of motion take the DAE form where the pair of matrices K, L ∈ R n,n , satisfying KL + LK = 0 and rank K L = n, describes the Dirac structure, and H x (x) = 1 2 x Qx defines the Hamilton function for some Q = Q .Denoting by ker and im the kernel and image of a linear map or its matrix representation, geometrically the Dirac structure is defined by the subspace D ⊂ X × X * given as ker K L , where X is the linear state space, and X * its dual space.
Remark 10 More precisely D ⊂ X × X * , where X is the tangent space to X at x ∈ X .However, since X is linear, tangent spaces at any x ∈ X can be identified with each other and with X .

Remark 11
In this paper we frequently switch between coordinate free representations with state space X and coordinate representations that are obtained for the case X = R n after choosing a basis.Furthermore, in this case it will be tacitly assumed that the dual basis is chosen for the dual space X * .Although this is not the most general setup (from an abstract linear algebraic or functional analytic point of view), it makes the presentation of results much more convenient.
Algebraic constraints occur if the matrix K is singular, and are represented via Qx ∈ im K .
Singularity of K often results from the network structure as the following example demonstrates.
Example 12 Consider a general linear LC-electrical circuit.Let the circuit graph be determined by an incidence matrix D, defining Kirchhoff's current laws I ∈ ker D and voltage laws V ∈ im D .Split the currents I into currents I C through the capacitors and currents I L through the inductors, and the voltages V into voltages V C across the capacitors and voltages V L across the inductors.Furthermore, let I C = − qC and V L = − φ, with q the vector of charges at the capacitors and ϕ the flux linkages of the inductors.Split the incidence matrix D accordingly as D = D C D L .Then Kirchhoff's current laws take the form Furthermore, let F be a maximal annihilator of D , i.e., ker F = im D .
Then Kirchhoff's voltage laws are given as F V C V L = 0, and after splitting F = F C F L accordingly, we have Writing the linear constitutive equations for the capacitors as q = CV C for some positive definite diagonal capacitance matrix C and those for the inductors as ϕ = LI L for some positive definite diagonal inductance matrix L, we finally obtain the system of equations which is in the form (6) with Hamilton function sponds to parallel interconnection of capacitors or series interconnection of inductors.(Note that since ker F = im D the linear space im F = ker D is spanned by the cycles of the circuit graph.)See e.g.[23] for pHDAE modeling of electrical circuits.
Motivated by [4] the port-Hamiltonian point of view on DAE systems was extended in [41] by replacing the gradient vector Qx of the Hamiltonian function H(x) = 1 2 x Qx by a general Lagrangian subspace, called Lagrange structure in the present paper, since we want to emphasize the similarity with Dirac structures.As we will see, this extension allows to bridge the gap with the (dissipative) Hamiltonian formulation (4) of the previous Section 2.
In the linear homogeneous case without dissipation, one starts with a Dirac structure D ⊂ X × X * and a Lagrange structure L ⊂ X × X * .The composition of the Dirac structure D and the Lagrange structure L, over the shared variables e ∈ X * , is then defined as In order to obtain a coordinate representation of the dynamics (8), the simplest option (later on in Section 5.1 we will discuss another one) is to start from the image representation of the Lagrange structure L, defined by a pair of matrices P, S ∈ R n,n satisfying P S = S P and rank P S = n.Taking coordinates x for X = R n and dual coordinates e for its dual space X * = R n , the image representation of L is given by for some parameterizing vector z in a space Z (of the same dimension as X ).
We will call this class extended Hamiltonian differential-algebraic systems (extended HDAEs).If dissipation is incorporated, see Section 4.2, then it is called extended dHDAEs.
Importantly, the presence of algebraic constraints in (10) may arise both by singularity of K (as was already the case for ( 6)) as well as by singularity of P (as was the case for ( 2)).This motivated the introduction of the notions of Dirac algebraic constraints (corresponding to singularity of K) and of Lagrange algebraic constraints (corresponding to singularity of P ) in [41].
Similar to dHDAE systems (2), the Hamiltonian of the extended HDAE system ( 10) is specified by the Lagrange structure, and is given by Indeed, one immediately has the energy conservation property and thus e f = 0.While singularity of K in physical systems modeling typically arises from interconnection due to the network structure, singularity of P often arises as a limiting situation.An elaborate example will be provided later as Example 17.
Remark 13 Note that if K is invertible, then by multiplying (10) with K −1 from the left, we obtain the lossless version of the dHDAE system (2) in Section 2 with x Qx by a Lagrange subspace (9) constitutes a first step towards an overarching formulation of (dissipative) Hamiltonian DAE systems.

Geometric theory of (dissipative) Hamiltonian DAEs
In this section we take a systematic geometric view on Hamiltonian DAE systems, extending the existing geometric treatment of extended HDAE systems, as already discussed in Section 3. We also incorporate the discussed classes of dHDAE systems from Section 2. Consider an n-dimensional linear state space X with elements denoted by x.Let X denote the tangent space to X at x ∈ X , with elements denoted by f and called flow vectors.As mentioned before, since X is linear, tangent spaces at different x ∈ X can be identified with each other and with X ; implying that f ∈ X as well.Furthermore, let X * be the dual space of X , with elements denoted by e and called effort vectors.

Dirac and Lagrange structures
The product space X × X * is endowed with the two canonical bilinear forms represented by the two matrices where we recognize Π − as the standard symplectic form on X .
Definition 14 A subspace D ⊂ X × X * is called a Dirac structure if the bilinear form •, • + is zero on D and moreover D is maximal with respect to this property.A subspace L ⊂ X × X * is a Lagrange structure if the bilinear form •, • − is zero on L and moreover L is maximal with respect to this property.A Lagrange structure L ⊂ X × X * is called nonnegative if the quadratic form defined by Π + is nonnegative on L.

Remark 15
In this paper we have chosen the terminology 'Lagrange structure', instead of the more common terminology 'Lagrangian subspace', in order to emphasize the similarity to Dirac structures.Also note that the definition of Dirac structures can be extended to manifolds instead of linear state spaces X ; in this context Dirac structures on linear spaces are often referred to as constant Dirac structures.
We have the following characterizations of Lagrange and Dirac structures, see e.g.[13,40].
Proposition 16 Consider an n-dimensional linear state space X and its dual space X * .i) A subspace D ⊂ X × X * is a Dirac structure if and only if D = D ⊥ ⊥ + , where ⊥ ⊥ + denotes the orthogonal complement with respect to the bilinear form •, • + .Furthermore, D ⊂ X × X * is a Dirac structure if and only e f = 0 for all (f, e) ∈ D and dim D = n.
ii) A subspace L ⊂ X × X * is a Lagrange structure if and only if L = L ⊥ ⊥ − , where ⊥ ⊥ − denotes the orthogonal complement with respect to the bilinear form •, • − .Any Lagrange structure satisfies dim L = n.
Dirac and Lagrange structures admit structured coordinate representations, see e.g.[40,41].For this paper the following representations are most relevant.Using matrices K, L ∈ R n,n , any Dirac structure D ⊂ X × X * admits the kernel/image representation with K, L satisfying rank K L = n and the generalized skew-symmetry condition Conversely any such pair K, L defines a Dirac structure.Analogously, any Lagrange structure L ⊂ X × X * can be represented as for certain matrices S, P ∈ R n,n satisfying rank P S = n as well as the generalized symmetry condition A Lagrange structure is, furthermore, nonnegative if and only if S P ≥ 0.
As already described in Section 3, by using the image representation x = P z, e = Sz of the Lagrange structure L, and the kernel representation Kf + Le = 0 of the Dirac structure D one is led to the representation (10) of the extended HDAE system defined by D and L.
The following is a physical example where both K and P turn out to be singular.The singularity of K is due to the presence of kinematic constraints, while the singularity of P is caused by a limiting argument in the energy expression.
Example 17 Consider two masses m 1 and m 2 connected by a spring with spring constant k, where the right mass m 2 is subject to the kinematic constraint v 2 = 0 (velocity is zero).With positions q 1 , q 2 and momenta p 1 , p 2 , the Hamiltonian is given by Denoting by e q1 , e q2 the spring forces at both ends of the spring, and by e p1 , e p2 the velocities of the two masses, we obtain the relation To consider the limit k → ∞, meaning that the spring is replaced by rigid connection, we first express the system in different coordinates.
This yields the transformed representation and the limiting Hamiltonian is just the kinetic energy The system has a Lagrange algebraic constraint due to the linear dependency in the rows of P .The Dirac structure D is given as After elimination of the Lagrange multiplier λ this yields and hence Finally subtracting the second equation from the first equation, we obtain the HDAE system Here the first equation is the Lagrange algebraic constraint z3 = 0 (and eventually z1 = 0) obtained by letting k → ∞ (corresponding to singularity of P ), and the last equation is the Dirac algebraic constraint − z3 m 2 + z4 m 1 +m 2 = 0, i.e., z4 = 0 resulting from the kinematic constraint, leading to singularity of K and resulting in the trivial dynamics ż2 = 0.
Remark 18 Instead of using a parametrization of the Lagrange structure L one can also use a parameterization of the Dirac structure D, with v ∈ V, where V is an n-dimensional parameter space.This yields an extended dHDAE system (but now in the parameter vector v) given by which is the adjoint system of (10).See [26] for a detailed discussion of adjoint systems of DAEs.

Incorporation of dissipation
As noted in Section 4, extended HDAE systems (10), geometrically defined by a Dirac and Lagrange structure, already include HDAE systems (2) without dissipation.Conversely, any HDAE system with K invertible can be rewritten into the form (2) with R = 0.In order to complete the geometric viewpoint towards the inclusion of dissipation (and thus to (2)), we recall the geometric definition of a port-Hamiltonian system [35,37,40].By replacing the Hamiltonian function by a Lagrange structure as in [41], and specializing to the case without external variables (inputs and outputs), such systems will be called extended dHDAE systems.
Definition 19 Consider a state space X with linear coordinates x and a linear space of resistive flows F R .Furthermore, consider a Dirac structure D on X ×F R , a Lagrange structure L ⊂ X ×X * , and a nonnegative Lagrange structure R ⊂ F R × F * R .Then an extended dissipative Hamiltonian DAE (extended dHDAE) system is defined as the tuple (X , F R , D, L, R) with If L is represented as in (15), i.e.
then it immediately follows from the properties of the Dirac structure D and the nonnegative Lagrange structure R that the dynamics of the extended dHDAE satisfies More generally we will now introduce the notion a maximally monotone subspace, which is overarching the notions of a Dirac structure D and a nonnegative Lagrange structure R.
Definition 20 Consider a linear space for all (f, e) ∈ M, and it is maximally monotone if additionally M is maximal with respect to this property (i.e., there does not exist a monotone subspace M ⊂ X × X * with M M ).

Remark 21
The definition of a monotone subspace is a special case of the notion of a monotone relation M , which is defined as a subset of X × X * satisfying (e 1 − e 2 ) (f for all (f 1 , e 1 ), (f 2 , e 2 ) ∈ M .Clearly if M is a subspace then (22) reduces to (21).(Maximally) monotone subspaces with a sign change were employed before in [21], using the terminology of '(maximally) linear dissipative relations'.Nonlinear port-Hamiltonian systems with respect to a general (maximally) monotone relation were coined as incremental port-Hamiltonian systems in [10]; see [9] for further developments.
Obviously, a subspace M is monotone if and only if the quadratic form defined by Π + is nonnegative on M, since (f, e), (f, e) + = 2e f .This yields Proposition 22 Consider a state space X with dim X = n.Then any monotone subspace of X ×X * has dimension less than or equal to n, and any maximally monotone subspace of X × X * has dimension n.Any maximally monotone subspace M ⊂ X × X * can be represented as for M, N ∈ R n,n satisfying rank N M = n and Conversely, any subspace defined by M, N satisfying (24) is a maximally monotone subspace.
Proof.The proof follows, since Π + has n positive and n negative eigenvalues.
Obviously any Dirac structure D given by a pair of matrices K, L is maximally monotone (take M = K and N = L).In a similar way any nonnegative Lagrange structure R given by a pair of matrices P, S with S P ≥ 0 is maximally monotone by taking N = P , M = S.
Importantly, also the composition of two maximally monotone subspaces is again maximally monotone.In order to prove this we first state the following lemma.
Lemma 23 Let A : F → G be a linear map between two linear spaces F, G. Let M G ⊂ G × G * be a maximally monotone subspace.Then the pull-back of M G via A, defined as is maximally monotone.Furthermore, let M F ⊂ F × F * be a maximally monotone subspace.Then the push-forward of M F via A, defined as and thus b A (M G ) is maximally monotone.The proof to show that f A (M F ) is maximally monotone is analogous.
Using Lemma 23 we can show that maximally monotone subspaces satisfy the following composition property.This same property was recently derived for maximally monotone relations in [9], assuming additional regularity conditions.
Proposition 24 Consider an extended dHDAE system as in (19) and let M a and M b be maximally monotone subspaces Proof.Let V a := F and V b := F. Define the linear maps Then for the maximally monotone subspace it can be readily checked that where M a × M I × M b is clearly maximally monotone.Then the proof finishes by applying Lemma 23.
We immediately have the following corollary.
Then the composition of D and R defined via Remark 26 We conjecture that conversely any maximally monotone subspace M can be generated this way, i.e., as the composition of a certain Dirac structure D and a certain nonnegative Lagrangian subspace R.
The presented analysis of maximally monotone subspaces leads to the following geometric definition of an extended dHDAE system, covering both dHDAE systems (2) and extended HDAE systems (10).See [21] for related results (using the terminology of (maximally) dissipative linear relations).
Definition 27 Consider a linear state space X with coordinates x, a maximally monotone subspace M ⊂ X × X * , and a Lagrange structure L ⊂ X × X * .Then an extended dHDAE system is a system (X , M, L) satisfying {( ẋ, x) | there exists e ∈ X * such that (− ẋ, e) ∈ M, (x, e) ∈ L}.
A coordinate representation of an extended dHDAE system is obtained as follows.Consider a coordinate expression (23) of the maximally monotone subspace M (with M, N satisfying ( 24)).This means that any element (f, e) ∈ M can be represented as for some v ∈ R n .Furthermore, any (x, e) ∈ L can be represented as in (20).
Substituting −f = ẋ = P ż this yields Then pre-multiplication by such a maximal annihilator C D eliminates the auxiliary variables v, and one obtains the coordinate representation Remark 28 The geometric construction of extended dissipative Hamiltonian system can be immediately generalized to extended dissipative port-Hamiltonian DAE (dpHDAE) systems with external port variables (inputs and outputs), by extending the maximally monotone subspace M ⊂ X × X * to a maximally monotone subspace M e ⊂ X ×X * ×F P ×F * P , where F P ×F * P is the space of external port variables.
Two particular cases of Definition 27 are of special interest.The first one is where the maximally monotone subspace M is actually a Dirac structure as in (6) with K, L satisfying (14).In this case one can take C D = K L , and thus the extended dHDAE system reduces to the extended HDAE system KP ż = LSz.
The other special case is where the maximally monotone subspace M in (23) is such that M is invertible.In this case, without loss of generality M can be taken to be the identity matrix, and the maximal annihilator C D can be taken to be of the form I D .Hence D = −N , and thus (27) reduces to P ż = −N Sz.
Furthermore, M N + N M ≥ 0 reduces to N + N ≥ 0, and hence Thus in this case the extended dpHDAE system takes the familiar form ( 2) with E = P , S = Q expressed as Remark 29 Similar to the theory exposed in [41] for HDAE systems (10), the algebraic constraints of the dpHDAE system (27) can be split into two classes: one corresponding to singularity of P (Lagrange algebraic constraints in [41]), and one corresponding to singularity of C. In case of (10) the second class of algebraic constraints are called Dirac algebraic constraints in [41], but now they correspond to the maximally monotone subspace.Furthermore, mimicking the developments in [41], one can transform algebraic constraints associated with index one belonging to one class into algebraic constraints in the other, by the use of additional state variables (serving as Lagrange multipliers).

Representation of DAE systems generated by
Dirac and Lagrange structures in the state variables x The representation (10) of an extended HDAE system as discussed in the previous sections does not use the state variables x of the state space X , but instead an equally dimensioned vector z ∈ Z parameterizing the Lagrange structure, cf. ( 9) and (15).In this section we show how a different DAE representation involving the original state vector x ∈ X can be obtained.Furthermore we discuss in what sense this representation in x is equivalent with the representation (10) involving z.

A coordinate representation in the original state variables x
Consider a Dirac structure D ⊂ X × X * , a Lagrange structure L ⊂ X × X * , and the resulting dynamics specified (in coordinate-free form) as D • L ⊂ X ×X .Let x be coordinates for the state space X and let the Dirac structure represented by a pair of matrices K, L and the Lagrange structure by a pair of matrices P, S. To derive a coordinate representation employing directly the state vector x, we first consider the combined representations of D and L, both in kernel representation, i.e., where e are dual coordinates for X * .In order to obtain a DAE system only involving x we need to eliminate the variables e.This can be done by considering a maximal annihilator (left null-space) M N of L −P , i.e., and thus, in particular, M L = N P . Since premultiplication of the equations ( 29) by M N thus yields Hence the resulting DAE system is given by Remark 30 Also for extended dHDAE systems (including dissipation) we can consider, instead of the coordinate representation (27) involving the parametrizing vector z, a representation that is using the original state x.
In fact, let as before, cf. ( 26 The analysis performed in the current subsection for (32) can be performed, mutatis mutandis, for (30) as well.
Recall that in the coordinate representation (10) we have the expression H z (z) = 1  2 z S P z for the Hamiltonian.In the representation (32) we do not yet have a Hamiltonian associated with the extended HDAE system.To define such a Hamiltonian H x in (32) we would need that P is invertible, in which case it is given by If P is invertible then there is a direct relation between the Hamiltonians (11) and (33).In fact, substituting x = P z, we immediately obtain Alternatively, if S is invertible then one can use the co-energy (Legendre transform) of H x given by for which H e (e) = H z (z) with e = Sz.Note that if P is invertible, then also M is invertible.This follows, since then M L = N P implies that the columns of N are in im M , and since rank M N = n this means M is invertible.The converse that M invertible implies P invertible follows analogously.In a similar fashion, it follows that L is invertible if and only if N is invertible.These observations imply the following simplifications of the representation (32) under additional assumptions.1. a) If P is invertible then N = M LP − and by multiplying (32) from the left by M −1 we obtain the system where the last equality follows from S P = P S. This is exactly the form of a Hamiltonian DAE system in case of a general Dirac structure and a Lagrange structure that is given as the graph of a symmetric matrix Q := SP −1 , see [35,38,39,40].Indeed, the Lagrange structure simplifies to the gradient of the Hamiltonian function H x (x) = 1 2 x SP −1 x. b) If in this case additionally K is invertible, then we obtain the Poisson formulation of Hamiltonian systems, see e.g.[2], with J = −J = K −1 L, and Q = Q .

a)
If L and thus also N is invertible, then by (31) we have M = N P L −1 and multiplying with N −1 from the left we get the DAE with b) If additionally P is invertible, then with Q = SP −1 = P − S this may be rewritten as which is the standard symplectic formulation of a Hamiltonian system in case additionally J is invertible, see e.g.[2].

Relation between the representations (10) and (32)
An immediate question that arises is how the representations (10) and (32) are related.We have already seen that if P is invertible then the relationship is obvious, since in this case x = P z defines an ordinary state space transformation.However, if P is not invertible then the representations are not state space equivalent, as the following simple example demonstrates.However, representations (10) (in the parameterizing z variables) and (32) (in the original state variables x) can be shown to be equivalent in the following generalized sense.First note that for any representation P, S of a Lagrange structure there exist nonsingular matrices V, W such that This is a direct consequence of Lemma 37 that will be presented in the next section.Setting and z 2 = e 2 .After such a transformation the system KP ż = LSz takes the form If we add to the vector x 1 e 2 the subvector x 2 , and if we consider the equations ( 40) together with the original Lagrange algebraic constraint x 2 = 0, then the so extended system can be rewritten as On the other hand, as shown in Subsection 5.1, the extended dHDAE system defined by the Dirac structure D and the Lagrangian structure L in the state space variables x can be expressed as with e 1 , e 2 serving as auxiliary variables.Instead of eliminating e 1 , e 2 from these equations, as discussed in Subsection 4.1, we can only eliminate e 1 by premultiplication of ( 42) by the full row rank matrix which directly leads to the system (41).This extended equivalence between ( 10) and ( 32) is summarized in the following proposition.
Proposition 32 Consider the pHDAE representations (10) and (32) defined by the same Lagrange structure L represented by matrices P, S, and by the same Dirac structure D represented by K, L. Consider a transformation such that P, S and K, L are transformed into the form (39) with corresponding partitioning where x 1 = z 1 .Adding to (10) the Lagrange algebraic constraint x 2 = 0 corresponding to x = P z, the resulting dHDAE system is given by (41).This system is equivalent to the representation (29) of (32) after elimination of the variables e 1 .
Note that the subvector e 2 = z 2 can be regarded as the Lagrange multiplier vector corresponding to the constraint x 2 = 0.As such, e 2 = z 2 does not contribute to the expression of the Hamiltonian H z (z).
Let us illustrate the previous discussion with some further examples.
Example 33 Consider a system in the form (10)

Equivalence transformations and condensed forms
To characterize the properties of extended dHDAEs we use transformations to condensed forms from which the properties can be read off.
For general DAEs (1) given by matrix pairs (E, A), E, AR ,n ( or the representation via matrix pencils λE − A) we can perform equivalence transformations of the coefficients of the form with U ∈ R ,ell , W ∈ R n,n nonsingular.This corresponds to a scaling of the equation with U and a change of variables z = W z. Under such transformations there is a one-to-one relationship between the solution spaces, see [25] and the canonical form is the Weierstraß canonical form.
For structured systems of the form (2), the associated equivalence transformation that preserves the structure is of the form n nonsingular.A condensed form for this case has been presented in [29].
Finally for systems of the form (10), the equivalence transformations have the form where The geometric interpretation of the set of transformations in (44) is clear: V defines a coordinate transformation on the state space X while V − is the corresponding dual transformation on the dual state space X * .Also note that the combination of V and V − on the product space X × X * leaves the canonical bilinear forms defined by the matrices Π − and Π + invariant.(In fact, it can be shown that any transformation on X × X * that leaves both canonical bilinear forms invariant is necessarily of this form for some invertible V .)Finally, U is an invertible transformation on the equation space for the kernel representation of the Dirac structure D, while W is an invertible transformation on the parametrization space Z for the Lagrange structure L.
In all three cases, in view of an implementation of the transformations as numerically stable procedures, we are also interested in the case that U, V, W are real orthogonal matrices.We then have that V −1 = V and for both pairs (K, L) and (P, S) this is a classical orthogonal equivalence transformation.
Using the described equivalence transformations we can derive condensed forms for pencils λP − S with P S = S P associated with Lagrange subspaces (or isotropic subspaces if the dimension is not n, see e.g.[29]).Here we slightly modify the representation and also give a constructive proof that can be implemented as numerically stable algorithm in Appendix A.
Lemma 37 Let P, S ∈ R n,m be such that P S = S P .Then there exist invertible matrices V ∈ R n,n , W ∈ R m,m such that with S 51 S 5,2 S 5,3 of full row rank n 5 .(Note that block sizes may be zero).Moreover, if the pencil λP − S is regular then the condensed form is unique, except for the order of blocks, and just contains the first four block rows and columns.

Proof. See Appendix A.
Note that the condensed form is in general not unique in the fifth block row, but the block sizes m 1 , m 2 , m 3 , m 4 and the row dimension n 5 are.
Corollary 38 Let P, S ∈ R n,m be such that P S = S P .Then there exist real orthogonal matrices V ∈ R n,n , W ∈ R m,m such that Proof.The proof follows by performing Steps 1. and 2. of the proof of Lemma 37, see Appendix B, which yields followed by a singular value decomposition V 3 Ŝ11 Ŵ3 = Š11 0 0 0 with Š11 nonsingular diagonal and a full rank decomposition V 4 Ŝ31 Ŵ3 = S 41 S 42 .
Corollary 38 shows that the characteristic quantities m 1 + m 2 , m 3 and m 4 , as well as n 5 can be obtained by purely real orthogonal transformations.The quantities m 1 , m 2 can then be determined from the real orthogonal staircase form of the symmetric pencil λP 11 S 11 − S 11 P 11 which has been presented in [8] and implemented as production software in [7].
There is an analogous condensed form for pencils of the form λK − L satisfying KL = −LK .For the case of regular pairs this directly follows from the canonical form presented in [11], but again we present the construction so that it can be directly implemented as a numerical method, see Appendix B.
Lemma 39 Let K, L ∈ R ,n be such that KL = −LK .Then there exist invertible matrices U ∈ R , , V ∈ R n,n such that  of full column rank n 5 .(Note that block sizes may be zero).
Moreover, if the pencil λK − L is regular then the condensed form is unique except for the order of blocks and just contains the first four block rows and columns.
Proof.See Appendix B. Note again that the form (47) is not unique in general but the block sizes Corollary 40 Let K, L ∈ R ,n be such that LK = −KL .Then there exist real orthogonal matrices U ∈ R , , V ∈ R n,n such that Proof.The proof follows by performing Steps 1. and 2. of the proof of Lemma 39, which yields followed by a singular value decomposition Û 3 L11 V4 = Ľ11 0 0 0 with Ľ11 nonsingular diagonal and a full rank decomposition Û 3 L13 V3 = L 14 L 24 .
Corollary 40 shows that the characteristic quantities 1 , 3 , 4 , as well as n 5 can be obtained by purely real orthogonal transformations.
The presented condensed forms can now be used in generating a condensed form for systems of the form (10).
Lemma 41 Consider a system of the form (10) with K, P, L, S ∈ R n,n and regular pencil λKP − LS.Then there exists invertible matrices U, V, W as in (44) such that where S 11 = S 11 and K 11 L 11 = −L 11 K 11 .
Proof.Since the pencil λKP − LS is square and regular, it is square, and also the pencil λP − S is regular, otherwise by Lemma 37 there would be common kernel of P and S which would imply the pencil λKP − LS to be singular.Thus, by Lemma 37 there exist nonsingular matrices with n 1 = m 1 + m 2 + m 3 , n 2 = m 4 and S11 symmetric.The regularity of the pencil λKP − LS implies that L12 L22 has full column rank and hence there exist invertible matrices U 2 ∈ R n,n , Ṽ2 ∈ R n 2 ,n 2 , and With we then get that Transforming the system as in (49) and setting partitioned accordingly, from the first block row of the coefficient matrices we obtain a reduced system given by with P = I n 1 , S = S = S 11 , K = K 11 and L = L 11 , together with an equation z 2 = K 21 ż1 , where z 2 does not contribute to the Hamiltonian H z (z) = 1 2 z P Sz.Note that the second equation is an index two constraint, because it uses the derivative of z 1 , [25].It arises from the Lagrange structure due to the singularity of P .
An analogous representation can be constructed from the condensed form of Lemma 39.
Lemma 42 Consider a system of the form (10) with K, P, L, S ∈ R n,n and regular pencil λKP − LS.Then there exist invertible matrices U, V, W as in (44) such that where L 11 = −L 11 and P 11 S 11 = S 11 P 11 .
Proof.Since the pencil λKP − LS is square and regular, also the pencil λK − L is regular, otherwise by Lemma 39 there would be a common left nullspace of K and L which would imply the pencil λKP −LS to be singular.Thus by Lemma 39, there exist nonsingular matrices with n 1 = 2 1 + 3 , n 2 = 4 and L11 skew-symmetric.The regularity of the pencil λKP − LS implies that S21 S22 has full row rank and hence there exist invertible matrices W 2 ∈ R n,n , Ṽ2 ∈ R n 2 ,n 2 , and With we then get that has the desired form with Transforming KP ż = LSz as in (51) and setting partitioned accordingly, from the first block row of the coefficient matrices we obtain a reduced system given by with P = P 11 , S = S 11 , K = I n 1 and L = L 11 = −L 11 , together with a differential algebraic equation 0 = z 2 , so that z 2 does not contribute to the Hamiltonian.
Remark 43 System (50) is a dHDAE of the form (35) in which the Lagrange structure is spanned by the columns of See also [30,32] for similar constructions in the context of removing the factor Q in systems of the form (2). Similarly, System (52) is a pHDAE of the form (2) (with R = 0) and the Dirac structure is spanned by the columns of .
We also perform a similar construction for systems of the form (27). Since C and D are chosen to be a maximal annihilator such that CN + DM = 0 in (26) and M N + N M ≥ 0 with rank N M = n we can use the same construction as in the proof of Lemma 39 to first transform N and M in such a way that This implies that we may choose C and D such that If λCP − DS is regular, then it follows that the last n 2 rows of CP have full row rank and hence altogether we have the following condensed form.
Lemma 44 Consider a system of the form (27) with C, P, L, S ∈ R n,n and regular pencil λCP − DS and C D a maximal annihilator as in (26).
Then there exist invertible matrices U, V, W as in (44) such that As a consequence, by transforming CP ż = DSz as in (53) and setting partitioned accordingly, from the second block row of the coefficient matrices we obtain ż2 = 0, i.e. z 2 is a constant function and the first block row gives an inhomogeneous reduced system with P = P 11 , S = S 11 , P S = S P , C = −M T 11 and D = I n 1 .It will then depend on the initial condition for z 2 whether z 2 = 0 in which case it does not contribute to the Hamiltonian, otherwise the Hamiltonian is still a quadratic function in z 1 plus some linear and constant terms.
Remark 45 The condensed forms in this section require rank decisions.Even if they are done in a numerically stable way using singular value decompositions, they can give wrong decisions in finite precision arithmetic.It is a common strategy to use in the case of doubt the worst case scenario.In the case of condensed forms this would be to assume that the problem is a DAE of index two.
In this section we have derived structured condensed forms and shown that these can also be used to identify a subsystem which is of one of the well-established forms plus an algebraic constraint whose solution does not contribute to the Hamiltonian.In the next section we analyze, when general DAEs can be transformed to the forms (2) or (10).For general DAE systems E ẋ = Ax it has been characterized in [30] when they are equivalent to a dHDAE system of the form (2). We present here a simplified result for the regular case.
Theorem 46 i) A regular pencil L(λ) = λ Ê − Â is equivalent to a pencil of the form λE − (J − R)Q as in (2) with λE − Q being regular if and only if the following conditions are satisfied: 1.The spectrum of L(λ) is contained in the closed left half plane.
2. The finite nonzero eigenvalues on the imaginary axis are semisimple and the partial multiplicities of the eigenvalue zero are at most two.
3. The index of L(λ) is at most two.
ii) A regular pencil L(λ) = λ Ê − Â is equivalent to a pencil of the form λE − (J − R) as in (2) (i.e., with Q = I) if and only if the following conditions are satisfied: 1.The spectrum of L(λ) is contained in the closed left half plane.
2. The finite eigenvalues on the imaginary axis (including zero) are semisimple.
3. The index of L(λ) is at most two.
As a Corollary for the case without dissipation we have the following result.
Corollary 47 A regular pencil L(λ) = λ Ê − Â is equivalent to a pencil of the form λE − J as in (2) (with Q = I, R = 0) if and only if the following conditions are satisfied: 1.All finite eigenvalues are on the imaginary axis and semisimple.
2. The index of L(λ) is at most two.
To study when general regular DAEs of the form (1) can be expressed as extended dHDAEs of the form (10) we first consider a condensed form under orthogonal equivalence.
Theorem 48 Consider a regular pencil λE −A with E, A ∈ R n,n of index at most two.Then there exist real orthogonal matrices U ∈ R n,n and V ∈ R n,n such that Proof.The proof is presented in Appendix C.
Transforming the DAE (1) as U EV V ẋ = U AV V x and setting V x = [x 1 , . . ., x 4 ] , it follows that x 1 = 0, x 2 is determined form the implicit ordinary differential equation (note that E 22 is invertible) x 3 = −A −1 33 A 32 x 2 , and x 4 is uniquely determined in terms of x 2 , ẋ2 , x 3 .Initial conditions can be prescribed freely for x 2 only.
Corollary 49 Consider a general regular pencil λE − A with E, A ∈ R n,n that is of index at most two and for which all finite eigenvalues are in the closed left half plane and those on the imaginary axis are semi-simple.Then there exist invertible matrices U ∈ R n,n and V ∈ R n,n such that ) where Proof.The proof follows by considering the condensed form (55), and using block elimination with the invertible matrices A 33 , A 14 , A 41 , E 22 to transform pencil in (55) to the form Let U 1 Ẽ11 V 1 = I n1 0 0 0 be the echelon form of Ẽ11 .We scale the first block row of with U 1 , the fourth block row with V −1 1 , the first block column by V 1 , and the fourth block column by U − 1 and obtain a form I n1 0 0 0 0 0 0 0 0 0 0 0 0 0 I n 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 one can multiply the second block row by X and obtain that E 33 = X and A 33 = X Â33 has the desired form, see e.g.[1,3].Note that the transformation to a system of the form (2) can also be achieved in a similar way for singular pencils with zero minimal indices.
If there is no dissipation, i.e. if R 33 = 0, then A 33 is skew-symmetric.In Corollary 49 we have shown that general systems E ẋ = Ax can be transformed to a very special canonical form and the following remark shows that for the case R 33 = 0 each of the blocks in the canonical form can be expressed as a pencil of the form λKP − LS.
Remark 50 Consider a regular pencil λE − A in the form (57) then after a permutation one gets four blocks which all can be written in the form λKP − LS as in (10).

1) We have
In this case we can insert the derivative of the second equation into the first (index reduction) and obtain 0 0 0 0 without changing the Hamiltonian.
2) We have with different possibilities of representation, e.g.a) Here the Hamiltonian is H z 4 = 1 2 z 4 z 4 .Note that the presented representations are in no way unique, but if the condensed form is available or computable, and the properties of Theorem 48 hold then we can express the general DAE in the representation (10) or (2).This discussion yields the following useful corollary.
Corollary 51 Consider a regular pencil of the form λKP − LS associated with the dHDAE (10).Then it has index at most two, and index two can only occur if the system has a singular Lagrange structure.
Proof.Consider the representations in Remark 50.Then the index two structure occurs only in the first case where K 1 , L 1 are invertible, but the product P T 1 S 1 is singular.Thus index two arises only from a singular Lagrange structure.

Conclusion and Outlook
Different definitions of (extended, dissipative) Hamiltonian or port Hamiltonian differential-algebraic systems lead to different representations.We have collected all the known representations as well as a few new ones and analyzed them from a geometric as well as an algebraic point of view.The latter leads to condensed forms that can be directly implemented in numerical algorithms to compute the structural properties of the systems.We have also studied the effect that the different representations have on the index of the differential-algebraic system as well as on the associated Hamilton function.In general it can be seen that certain algebraic constraints do not contribute to the Hamiltonian and therefore can be separated from the system in an appropriate coordinate system.We have also characterized when a general differential-algebraic system can be transformed to the different representations.Several important tasks remain open.These include extensions to the case of non-regular systems.These can be based on the results and methods in Appendices A and B that are already proved for the non-square case.For systems with inputs and outputs, linear time-varying and nonlinear systems the extensions are currently under consideration.
, W 2 = I m1 0 0 W2 and form and form where by the structure (16) now Š11 is symmetric.Note that although we are working with nonorthogonal transformation matrices in this step the numerical errors can be controlled, since we are inverting diagonal matrices.
where K11 ∈ R ˜ 1 , ˜ 1 is diagonal with positive diagonal elements.This transformation can be constructed via a singular value decomposition of K and a numerical rank decision.Then, using the structure ( 14), it follows that K 11 L11 = L 11 K11 and L21 = 0.
Step 2. Let Ũ2 , Ṽ2 be real orthogonal matrices such that where L22 ∈ R 4 , 4 is diagonal with positive diagonal elements.This transformation can be constructed again via a singular value decomposition and a numerical rank decision.Set where by the structure (14) now Ľ11 is skew-symmetric.Note that although we are working with nonorthogonal transformation matrices in this step, the numerical errors can be controlled since we are inverting diagonal matrices.
Step 2. Let Ũ2 , Ṽ2 be real orthogonal matrices such that Ũ 2 Ã22 Ṽ2 = Â22 0 0 0 , where Â22 ∈ R n 3 ,n 3 is diagonal with positive diagonal elements.This transformation can be constructed again via a singular value decomposition and a numerical rank decision.Set and form with Ê11 = Ẽ11 .The regularity of the pencil implies that A 13 has full column rank n 1 = ñ1 − n 3 and that A 31 has full row rank n 1 = ñ1 − n 3 , because otherwise there would be common right or left nullspace, respectively.
Step 3. Let Ũ31 , Ṽ31 , Ũ13 , Ṽ13 be real orthogonal matrices of appropriate dimensions such that Ũ 13 Ã13 Ṽ13 = Â14 0 , Ũ 31 Ã31 Ṽ31 = Â41 0 where Â41 ∈ R n 1 ,n 1 and Â14 ∈ R n 1 ,n 1 are diagonal with positive diagonal elements.This transformation can be constructed again via a singular value decomposition and numerical rank decisions.Set 3 are as claimed in (55).The invertibility of E 22 then follows from the assumption that the pencil has index at most two, see [25].
Analogously, we consider a kernel representation of the Dirac structure D given by matrices K, L satisfying KL = −LK and rank K L = n, such that D = {(f, e) ∈ X × X * | Kf + Le = 0}.(More details will be given in Section 4.) Substituting f = − ẋ = P ż and e = Sz this leads to the coordinate representation KP ż = LSz.
R ) ∈ D, and (f R , e R ) ∈ R} is maximally monotone.In particular, for any (f, e) ∈ D • R, one has e f = e R f R ≥ 0.
), C D denote a maximal annihilator of N M , i.e., ker C D = im N M .Then consider, similarly to (29), the stacked matrix annihilator V W to D −P , that is ker V W = im D −P .Then premultiplication by V W yields the representation V C ẋ = W S x

V 11
form of the skew-symmetric matrix Š11 under congruence which can be obtained by first computing the spectral decomposition and then scaling the nonsingular diagonal parts by congruence to be ±I, see e.g.[27].Furthermore let Š31 = V4 S 51 S 5,2 S 5,3 0 0 0 be a full rank decomposition partitioned accordingly, with V3 real orthogonal.Then set form of the skew-symmetric matrix Ľ11 under congruence which can be obtained by first computing the spectral decomposition and then scaling the nonsingular diagonal parts by congruence to be ±I.This procedure is implemented in a numerically robust way in [5rank decomposition partitioned accordingly, with V3 real orthogonal.Then set and where L 11 = −L 11 and P 11 S 11 = S 11 P 11 .