Structure-preserving discretization for port-Hamiltonian descriptor systems

TL;DR

This paper extends port-Hamiltonian descriptor systems to more general forms, explores structure-preserving discretization methods, and demonstrates their effectiveness through numerical examples, enhancing modeling and control capabilities.

Contribution

It introduces a generalized framework for port-Hamiltonian systems, including discretization techniques that preserve their energy structure and stability properties.

Findings

Structure-preserving discretization schemes maintain energy properties.

Port-Hamiltonian systems are invariant under nonlinear transformations.

Numerical examples confirm the effectiveness of the proposed methods.

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. Numerical examples are provided.

Equations

$$\begin{array}{rl}E(t,x)\dot{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){)}^{T}z(t,x)+(S(t,x)-N(t,x))u,\end{array}$$

$$\Gamma :=\left[\begin{array}{cc}J& B\\ -B^T& N\end{array}\right],W:=\left[\begin{array}{cc}R& P\\ P^T& S\end{array}\right]$$

$$\frac{\text{d}}{\text{d}t}\mathcal{H}(t,x(t))=-{\left[\begin{array}{c}z\\ u\end{array}\right]}^{T}W\left[\begin{array}{c}z\\ u\end{array}\right]+u^{T}y$$

$$\mathcal{H}(t_2,x(t_2))-\mathcal{H}(t_1,x(t_1))\le {\int }_{t_1}^{t_2}u(\tau {)}^{T}y(\tau )\text{d}\tau$$

$$\left[\begin{array}{c}E\dot{x}+r\\ 0\end{array}\right]=(\Gamma -W)\left[\begin{array}{c}z\\ u\end{array}\right]+\left[\begin{array}{c}0\\ y\end{array}\right],$$

$$\frac{\text{d}}{\text{d}t}\mathcal{H}(t,x(t))$$

$$={\left[\begin{array}{c}z\\ u\end{array}\right]}^T\left((\Gamma -W)\left[\begin{array}{c}z\\ u\end{array}\right]+\left[\begin{array}{c}0\\ y\end{array}\right]\right)=-{\left[\begin{array}{c}z\\ u\end{array}\right]}^{T}W\left[\begin{array}{c}z\\ u\end{array}\right]+u^{T}y,$$

$$\begin{array}{rl}\tilde{E}\dot{\tilde{x}}+\tilde{r}& =(\tilde{J}-\tilde{R})\tilde{z}+(\tilde{B}-\tilde{P})u,\\ y& =(\tilde{B}+\tilde{P}{)}^T\tilde{z}+(S-N)u,\end{array}$$

$$\begin{array}{rl}\tilde{W}& =\left[\begin{array}{cc}\tilde{R}& \tilde{P}\\ {\tilde{P}}^T& S\end{array}\right]=\left[\begin{array}{cc}{U}^T(R\circ \phi )U& U^T(P\circ \phi )\\ (P\circ \phi {)}^{T}U& S\circ \phi \end{array}\right]=\\ & ={\left[\begin{array}{cc}U& 0\\ 0& I\end{array}\right]}^T(W\circ \phi )\left[\begin{array}{cc}U& 0\\ 0& I\end{array}\right]\ge 0.\end{array}$$

$$\begin{array}{rl}\frac{\partial \tilde{\mathcal{H}}}{\partial \tilde{x}}(t,\tilde{x})& ={\frac{\partial \phi }{\partial \tilde{x}}}^T\left(\frac{\partial \mathcal{H}}{\partial x}\circ \phi \right)={\frac{\partial \phi }{\partial \tilde{x}}}^T\left(E^{T}z\circ \phi \right)=\\ & ={\frac{\partial \phi }{\partial \tilde{x}}}^T(E^T\circ \phi )U{U}^{-1}(z\circ \phi )={\tilde{E}}^T\tilde{z}\end{array}$$

$$\begin{array}{rl}\frac{\partial \tilde{\mathcal{H}}}{\partial t}(t,\tilde{x})& =\frac{\partial \mathcal{H}}{\partial t}\circ \phi +{\left(\frac{\partial \mathcal{H}}{\partial x}\circ \phi \right)}^T\frac{\partial \phi }{\partial t}=\\ & =z^{T}r\circ \phi +\left(z^{T}E\circ \phi \right)\frac{\partial \phi }{\partial t}=\\ & =(z\circ \phi {)}^T\left(r\circ \phi +(E\circ \phi )\frac{\partial \phi }{\partial t}\right)=\\ & ={\tilde{z}}^T{U}^T{U}^{-T}\tilde{r}={\tilde{z}}^T\tilde{r}.\end{array}$$

$$\begin{array}{rl}\left[\begin{array}{cc}E& r\\ 0& 1\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{t}\end{array}\right]& =\left[\begin{array}{cc}J-R& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}z\\ 0\end{array}\right]+\left[\begin{array}{cc}B-P& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}u\\ 1\end{array}\right],\\ \left[\begin{array}{c}y\\ 0\end{array}\right]& ={\left[\begin{array}{cc}B+P& 0\\ 0& 1\end{array}\right]}^T\left[\begin{array}{c}z\\ 0\end{array}\right]+\left[\begin{array}{cc}S-N& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}u\\ 1\end{array}\right].\end{array}$$

$$E_i{\dot{x}}_i$$

$$y_i$$

$$\left[\begin{array}{ccc}E& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{\hat{u}}\\ \dot{\hat{y}}\end{array}\right]$$

$$y$$

$$\Pi =\left[\begin{array}{cccc}{I}_{{\ell }_1}& 0& 0& 0\\ 0& 0& I_{{\ell }_2}& 0\\ 0& I_{n_1}& 0& 0\\ 0& 0& 0& I_{n_2}\end{array}\right].$$

$$\text{set@color}⟨\text{leavevmode}\text{set@color}⟨(f_1,e_1),(f_2,e_2)\text{set@color}⟩\text{leavevmode}\text{set@color}⟩:=⟨e_1|f_2⟩+⟨e_2|f_1⟩,$$

$${\mathcal{D}}_x=\left\{(f,e)\in {\mathcal{V}}_x\times {\mathcal{V}}_x^{\ast }:f+J(x)e=0\right\},$$

$$\text{set@color}⟨\text{leavevmode}\text{set@color}⟨(f,e),(f^{\prime },e^{\prime })\text{set@color}⟩\text{leavevmode}\text{set@color}⟩$$

$$=⟨e|f^{\prime }⟩-⟨e^{\prime }|Je⟩=⟨e|f^{\prime }+J{e}^{\prime }⟩.$$

$$\displaystyle\mathcal{D}_{x}=\Bigg{\{}(f,e)\in\mathcal{V}_{x}\times\mathcal{V}_{x}^{*}\;\Bigg{|}\;f+\begin{bmatrix}\Gamma(x)&I_{\ell+m}\\ -I_{\ell+m}&0\end{bmatrix}e=0\Bigg{\}}$$

$$f_s$$

$$f_p$$

$$e_d$$

$$\left[\begin{array}{c}E(x)\dot{x}\\ -y\end{array}\right]=(\Gamma (x)-W(x))\left[\begin{array}{c}z(x)\\ u\end{array}\right].$$

$$\left[\begin{array}{c}-f_s\\ -f_p\end{array}\right]=\Gamma (x)e_s+e_p,f_d=(e_s,e_p),$$

$$\begin{array}{rl}\left[\begin{array}{c}E(x)\dot{x}\\ -y\end{array}\right]& =(\Gamma (x)-W(x))\left[\begin{array}{c}z(x)\\ u\end{array}\right],\\ f& =\left(-E(x)\dot{x},y,\left[\begin{array}{c}z(x)\\ u\end{array}\right]\right),\\ e& =\left(z(x),u,-W(x)\left[\begin{array}{c}z(x)\\ u\end{array}\right]\right),\end{array}$$

$$\begin{array}{rl}0& =⟨z(x)|-E(x)\dot{x}⟩+⟨u|y⟩+⟨-W(x)\left[\begin{array}{c}z(x)\\ u\end{array}\right]|\left[\begin{array}{c}z(x)\\ u\end{array}\right]⟩=\\ & =\frac{\text{d}}{\text{d}t}\mathcal{H}(x)+y^{T}u-{\left[\begin{array}{c}z(x)\\ u\end{array}\right]}^{T}W(x)\left[\begin{array}{c}z(x)\\ u\end{array}\right],\end{array}$$

$${\ell }_i(\tau ):=\prod _{\begin{array}{c}j=1\\ j\ne i\end{array}}^s\frac{\tau -{\gamma }_j}{{\gamma }_i-{\gamma }_j}.$$

$$\begin{array}{rl}\dot{\tilde{x}}(t_0+\tau h)& =\sum _{i=1}^s{\dot{x}}_i{\ell }_i(\tau ),\\ \tilde{x}(t_0+\tau h)& =x_0+h\sum _{j=1}^s{\dot{x}}_j{\int }_0^{\tau }{\ell }_j(\sigma )\text{d}\sigma ,\end{array}$$

$$x_i$$

$$x_f$$

$${\mathcal{D}}_{x_i}$$

Full text

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 $\mathbb{I}\subseteq\mathbb{R}$, a state space $\mathcal{X}\subseteq\mathbb{R}^{n}$ and the extended state space $\mathcal{S}=\mathbb{I}\times\mathcal{X}$. A port-Hamiltonian descriptor system (or pHDAE) is a system of differential (-algebraic) equations of the form

[TABLE]

associated with an Hamiltonian function $\mathcal{H}\in\mathcal{C}^{1}(\mathcal{S},\mathbb{R})$, where $x(t)\in\mathcal{X}$ is the state, $u(t),y(t)\in\mathbb{R}^{m}$ are the input and output, $E\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{\ell,n})$ is the flow matrix, $r,z\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{\ell})$ are the time-flow and effort functions, $J,R\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{\ell,\ell})$ are the structure and dissipation matrices, $B,P\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{\ell,m})$ are the port matrices and $S,N\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{m,m})$ are the feed-through matrices. Furthermore, the following properties must hold:

The extended structure and dissipation matrices $\Gamma,W\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{\ell+m,\ell+m})$, defined as

[TABLE]

satisfy $\Gamma=-\Gamma^{T}$ and $W=W^{T}\geq 0$ pointwise. 2. 2.

The gradient of the Hamiltonian satisfies $\partial_{x}\mathcal{H}=E^{T}z$ and $\partial_{t}\mathcal{H}=z^{T}r$ 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)

[TABLE]

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

[TABLE]

holds.

Proof.

The pHDAE (1) can be written as

[TABLE]

in particular

[TABLE]

which is (3). By time integration, since $W=W^{T}\geq 0$, we immediately get (4). ∎

Note that, with no input ($u=0$), if the Hamiltonian $\mathcal{H}$ is locally positive-definite in an equilibrium point $x^{*}$ (up to shifting it by a constant), then $\mathcal{H}$ 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 $\tilde{\mathcal{X}}\subseteq\mathbb{R}^{\tilde{n}}$ be a second state space, let $\tilde{\mathcal{S}}:=\mathbb{I}\times\tilde{\mathcal{X}}$, let $x=\varphi(t,\tilde{x})\in\mathcal{C}^{1}(\tilde{\mathcal{S}},\mathcal{X})$ be a local diffeomorphism (with respect to $\tilde{x}$) and let $U\in\mathcal{C}(\tilde{\mathcal{S}},\mathbb{R}^{\ell,\ell})$ be pointwise invertible. Consider the input-output DAE

[TABLE]

with $\tilde{E}=U^{T}(E\circ\varphi)\partial_{\tilde{x}}\varphi$, $\tilde{J}=U^{T}(J\circ\varphi)U$, $\tilde{R}=U^{T}(R\circ\varphi)U$, $\tilde{B}=U^{T}(B\circ\varphi)$, $\tilde{P}=U^{T}(P\circ\varphi)$, $\tilde{z}=U^{-1}(z\circ\varphi)$ and $\tilde{r}=U^{T}(r\circ\varphi+(E\circ\varphi)\partial_{t}\varphi)$, where we denote $(F\circ\varphi)(t,\tilde{x})=F(t,\varphi(t,\tilde{x}))$ for any $F\in\mathcal{C}(\mathcal{S},\cdot)$, and let $\tilde{\mathcal{H}}(t,\tilde{x}):=(\mathcal{H}\circ\varphi)(t,\tilde{x})$. Then (5) is a pHDAE with Hamiltonian function $\tilde{\mathcal{H}}$, and to any solution $(\tilde{x},u,y)$ of (5) corresponds a solution $(x,u,y)$ of (1) with $x(t)=\varphi(t,\tilde{x}(t))$. Furthermore, if $\varphi(t,\cdot)$ is a global diffeomorphism for all $t\in\mathbb{I}$, then the two systems are equivalent.

Proof.

The transformed DAE system is obtained from (1) by setting $x=\varphi(t,\tilde{x})$, pre-multiplying with $U^{T}$ and inserting $UU^{-1}$ in front of $z$ in the first equation. It is then clear that if $(\tilde{x},u,y)$ is a solution of (5), then $(x,u,y)$ is a solution of the original system. If $\varphi(t,\cdot)$ is a global diffeomorphism, then we can apply the inverse tranformation $(\varphi^{-1},U^{-1})$, where $\varphi^{-1}$ is to intended as $\varphi(t,\varphi^{-1}(t,x))=x$, and get a solution $(\tilde{x},u,y)$ for any solution $(x,u,y)$ of the original system, making the two systems equivalent.

To show that (5) is still a pHDAE, we must check that conditions 1 and 2 in the definition are satisfied.

By substitution, we get

[TABLE] 2. 2.

By substitution, we get

[TABLE]

and

[TABLE]

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:

[TABLE]

Note that condition 2 becomes simply $\nabla_{\tilde{x}}\tilde{\mathcal{H}}=\tilde{E}^{T}\tilde{z}$, where the tilde denotes the quantities in the autonomous system.

2.2.4 Structure-preserving interconnection

Let us consider two autonomous pHDAEs of the form

[TABLE]

with Hamiltonian $\mathcal{H}_{i}$, for $i=1,2$, and assume that the aggregated input $u=(u_{1},u_{2})$ and output $y=(y_{1},y_{2})$ satisfy a linear interconnection relation $Mu+Ny=0$ for some $M,N\in\mathcal{C}^{*}(\mathcal{S},\mathbb{R}^{k,m})$. Then the aggregated system can be written as a pHDAE of the form

[TABLE]

with Hamiltonian $\mathcal{H}=\mathcal{H}_{1}+\mathcal{H}_{2}$, where we have introduced new state variables $\hat{u},\hat{y}\in\mathbb{R}^{n}$ that copy $u,y$, and we aggregated $E=\operatorname{diag}(E_{1},E_{2})$, $x=(x_{1},x_{2})$, $z=(z_{1},z_{2})$, $m=m_{1}+m_{2}$, $\Gamma=\Pi\operatorname{diag}(\Gamma_{1},\Gamma_{2})\Pi^{T}$ and $W=\Pi\operatorname{diag}(W_{1},W_{2})\Pi^{T}$, where $\Pi\in\mathbb{R}^{\ell+m,\ell+m}$ is the permutation matrix

[TABLE]

Note that, in general, we may not be able to reduce the number of inputs and outputs. If we also assume that the interconnection is energy-preserving (e.g. if $Mu+Ny=0$ defines a Dirac structure for $(y,u)$), then index reduction [5] and row operations can usually be applied to make the system smaller.

2.3 Dirac structure

Port-Hamiltonian systems are ofter described through differential geometric structures known as Dirac structures [11]. We present first the basic definitions.

Definition 2** (linear Dirac structure).**

Let $\mathcal{F}$ be a linear space and $\mathcal{E}:=\mathcal{F}^{*}$ its dual space. Let \mathopen{\hbox{\set@color{\langle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\langle}}}\cdot,\cdot\mathclose{\hbox{\set@color{\rangle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\rangle}}} be the bilinear form on $\mathcal{F}\times\mathcal{E}$ defined as

[TABLE]

where $\langle\cdot\,|\,\cdot\rangle$ denotes the duality pairing. A Dirac structure on $\mathcal{F}\times\mathcal{E}$ is then a linear subspace $\mathcal{D}\subseteq\mathcal{F}\times\mathcal{E}$, such that \mathcal{D}=\mathcal{D}^{\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}}.

In particular, in finite dimension, one only needs to prove $\dim\mathcal{D}=\dim\mathcal{F}$ and $\langle e\,|\,f\rangle=0$ for all $(f,e)\in\mathcal{D}$. If $(f,e)\in\mathcal{D}$, then $f$ and $e$ are called flow and effort, respectively. In [11], the more general definition of modulated Dirac structure over $\mathcal{X}$ is also presented, as a subbundle of $T\mathcal{X}\oplus T^{*}\mathcal{X}$, that denotes the Whitney sum between the tangent and cotangent bundles of $\mathcal{X}$. To introduce a Dirac structure for pHDAEs, we need to extend this further.

Definition 3**.**

Consider a state space $\mathcal{X}$ and a vector bundle $\mathcal{V}$ over $\mathcal{X}$ with fibers $\mathcal{V}_{x}$. A Dirac structure over $\mathcal{V}$ is a subbundle $\mathcal{D}\subseteq\mathcal{V}\oplus\mathcal{V}^{*}$ such that, for all $x\in\mathcal{X}$, $\mathcal{D}_{x}\subseteq\mathcal{V}_{x}\times\mathcal{V}_{x}^{*}$ is a linear Dirac structure.

Note that modulated Dirac structures are a special case of our definition, where $\mathcal{V}=T\mathcal{X}$. To associate a Dirac structure to our pHDAE system, we first prove the following lemma:

Lemma 1**.**

Let $\mathcal{D}\subseteq\mathcal{V}\oplus\mathcal{V}^{*}$ be a vector subbundle with fibers defined by

[TABLE]

where $J:\mathcal{X}\to\mathcal{L}(\mathcal{V}_{x}^{*},\mathcal{V}_{x})$ is a skew-symmetric operator. Then $\mathcal{D}$ is a Dirac structure.

Proof.

For generic $(f,e)\in\mathcal{D}_{x}$ and $(f^{\prime},e^{\prime})\in\mathcal{V}_{x}\times\mathcal{V}_{x}^{*}$,

[TABLE]

We show that $(f^{\prime},e^{\prime})\in\mathcal{D}_{x}$ if and only if $\langle e\,|\,f^{\prime}+Je^{\prime}\rangle=0$ for all $(f,e)\in\mathcal{D}_{x}$. On one hand, if $(f^{\prime},e^{\prime})\in\mathcal{D}_{x}$, then $\langle e\,|\,f^{\prime}+Je^{\prime}\rangle=0$ holds for any $e\in\mathcal{E}$. On the other hand, if $(f^{\prime},e^{\prime})\notin\mathcal{D}_{x}$, then $f^{\prime}+Je^{\prime}\neq 0$ and so $\exists e\in\mathcal{E}$ such that $\langle e\,|\,f^{\prime}+Je^{\prime}\rangle=1$, but $(f,e)\in\mathcal{D}_{x}$ with $f=-Je$. ∎

We can now associate a Dirac structure to our pHDAE system in the following way:

Theorem 3**.**

Given an autonomous pHDAE, let us define the flow fiber $\mathcal{V}_{x}=\mathcal{F}_{x}^{s}\times\mathcal{F}_{x}^{p}\times\mathcal{F}_{x}^{d}$ for all $x\in\mathcal{X}$, where $\mathcal{F}_{x}^{s}:=E(x)T_{x}\mathcal{X}\subseteq\mathbb{R}^{\ell}$ is the storage flow fiber, $\mathcal{F}_{x}^{p}:=\mathbb{R}^{m}$ is the port flow fiber and $\mathcal{F}_{x}^{d}:=\mathbb{R}^{\ell+m}$ is the dissipation flow fiber. Let us partition $f=(f_{s},f_{p},f_{d})\in\mathcal{V}$ and $e=(e_{s},e_{p},e_{d})\in\mathcal{V}^{*}$. Then the subbundle $\mathcal{D}\subseteq\mathcal{V}\oplus\mathcal{V}^{*}$ with

[TABLE]

is a Dirac structure over $\mathcal{V}$. Furthermore, the system of equations

[TABLE]

is equivalent to the original pHDAE, and $\langle e\,|\,f\rangle=0$ represents the power balance equation.

Proof.

$\mathcal{D}$ is a Dirac structure because of Lemma 1. Note that the pHDAE can be written in compact form as

[TABLE]

The condition $(f,e)\in\mathcal{D}_{x}$ can be written as

[TABLE]

that together with the conditions (7) is equivalent to

[TABLE]

which is exactly the compact form of the pHDAE, plus the definition for the flow and effort. Finally, note that the equation $\langle e\,|\,f\rangle=0$ can be written as

[TABLE]

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 $\mathcal{H}(x)$ are needed. These can also be lifted to a geometric interpretation, by the means of a Lagrangian submanifold and of a dissipative structure [12].

3 Time discretization

Let us consider a finite-dimensional pHDAE of the form (1). Under some regularity assumptions [5], many classes of Runge-Kutta methods can be applied to compute time-discretizations of differential-algebraic equations. An important class of Runge-Kutta methods is the class of collocation methods. We extend the results from [4] to pHDAEs, proceeding in a similar way.

3.1 Collocation methods for DAEs

Assume that an input function $u:\mathbb{I}\times\mathcal{X}\to\mathbb{R}^{m}$, depending on time and possibly space, is given. Given a time interval $\mathbb{I}=[t_{0},t_{f}]$ of length $h=t_{f}-t_{0}$, and a consistent initial condition $x_{0}$ at time $t_{0}$, we approximate the solution $x(t)$ of the pHDAE on $\mathbb{I}$ with a polynomial $\tilde{x}(t)\in\mathbb{R}[t]_{s}$, where $\mathbb{R}[t]_{s}$ denotes the vector space of polynomials of degree at most $s$. The polynomial is chosen such that $\tilde{x}(t_{0})=x_{0}$, and that it satisfies the differential-algebraic equation in $s$ collocation points $t_{i}=t_{0}+h\gamma_{i}$ with $\gamma_{i}\in[0,1]$ for $i=1,\ldots,s$.

Let $\ell_{i}$ denote the $i$-th Lagrange interpolation polynomial with respect to the nodes $\gamma_{1},\ldots,\gamma_{s}$, i.e.,

[TABLE]

Then we can write

[TABLE]

for certain $\dot{x}_{i}=\dot{\tilde{x}}(t_{i})$ that we have to compute, and also

[TABLE]

where $\alpha_{ij}:=\int_{0}^{\gamma_{i}}\ell_{j}(\sigma)\textup{d}\sigma$ and $\beta_{j}:=\int_{0}^{1}\ell_{j}(\sigma)\textup{d}\sigma$. In particular, the constants $\alpha_{ij},\beta_{j},\gamma_{i}\in\mathbb{R}$ for $i,j=1\ldots s$ are the coefficients of the Butcher diagram of the associated Runge-Kutta method [3].

3.2 Dirac structure associated to discretization

Consider for simplicity the autonomous case. The non-autonomous case is similar, with a more cumbersome notation. Let $\mathcal{D}_{x}$ be the Dirac structure associated to (6) (as in Theorem 3). We define the Dirac structure associated to the time-discretization as $\left\{\mathcal{D}_{x_{i}}:i=1,\ldots,s\right\}$, i.e.,

[TABLE]

with $f^{i}=(f_{s}^{i},y_{i},f_{d}^{i})$ and $e^{i}=(e_{s}^{i},u_{i},e_{d}^{i})$. Taking $(f^{i},e^{i})\in\mathcal{D}_{x_{i}}$, together with

[TABLE]

we get a system that is equivalent to applying the collocation method and computing the discrete input and output $u_{i},y_{i}$, for $i=1,\ldots,s$. Let us define the collocation flows, efforts, input and output in $\mathbb{R}[t]_{s-1}$ as

[TABLE]

Note that, by construction, $(\tilde{f}_{s},\tilde{y},\tilde{f}_{d};\tilde{e}_{s},\tilde{u},\tilde{e}_{d})\in\mathcal{D}_{\tilde{x}}$ in all collocation points $t_{i}$. Let us consider the evolution of the Hamiltonian $\mathcal{H}$ along the collocation polynomial $\tilde{x}(t)$. In particular, let $H(t):=\mathcal{H}(\tilde{x}(t))$: we can then write

[TABLE]

In the collocation points, the PBE is satisfied:

[TABLE]

for $i=1,\ldots,s$. In particular, if we apply the quadrature rule associated to the collocation method, we get

[TABLE]

where $p\in\mathbb{N}$ is the degree of exactness of the quadrature rule. We observe that, for the same reason,

[TABLE]

so

[TABLE]

We note that, if $p\geq 2s-2$, then equations (9a) and (9b) are exact. Furthermore, if $\beta_{j}\geq 0$ for $j=1,\ldots,s$ (and this is the case for many collocation methods), we can deduce that $h\sum_{j=1}^{s}\langle e_{d}^{j}\,|\,f_{d}^{j}\rangle\leq 0$, thus the dissipation component retains its qualitative behaviour.

3.3 The quadratic Hamiltonian case

Let us consider the case where the Hamiltonian $\mathcal{H}(x)$ is a polynomial of degree (at most) 2, i.e. it can be written as

[TABLE]

for some $Q=Q^{T}\in\mathbb{R}^{n,n}$, $v\in\mathbb{R}^{n}$ and $c\in\mathbb{R}$. Since $\tilde{x}\in\mathbb{R}[t]_{s}^{n}$, we also have $H=\mathcal{H}\circ\tilde{x}\in\mathbb{R}[t]_{2s}$ and $\dot{H}\in\mathbb{R}[t]_{2s-1}$. It is known that the maximum degree of exactness for quadrature rules with $s$ nodes is $2s-1$, and that it is attained only with Gaussian quadrature rules, that are associated with Gauss-Legendre collocation methods. In particular, if we apply such a method, the integration of $\dot{H}$ will be exact, i.e.

[TABLE]

In particular, since we always have $\beta_{j}\geq 0$ for Gauss-Legendre collocation, the dissipation term is always non-positive, and we can write the discrete dissipation inequality

[TABLE]

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 $R_{G},R_{L},R_{R}>0$ are resistances, $L>0$ is an inductor, $C_{1},C_{2}>0$ are capacitors and $E_{G}$ is a controlled voltage source. This can be interpreted as a basic representation of a DC generator ($E_{G}$,$R_{G}$), connected to a load ($R_{R}$) with a transmission line (the $\pi$ model given by $C_{1},C_{2},L,R_{L}$).

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

[TABLE]

The energy in the system can be stored in the inductor and in the two capacitors, giving the Hamiltonian

[TABLE]

The system can then be written equivalently as an autonomous pHDAE of the form

[TABLE]

with $x=(I,V_{1},V_{2},I_{G},I_{R})$, $E=\operatorname{diag}(L,C_{1},C_{2},0,0)$, $B=e_{4}$, $u=E_{G}$, $y=I_{G}$ and

[TABLE]

The power balance equation reads as

[TABLE]

in particular, if we shut down the generator ($E_{G}=0$), the Hamiltonian will decrease and converge to a solution such that $\dot{\mathcal{H}}=0$, that is, $I=I_{G}=I_{R}=0$. The only state compatible with (10) that satisfies that condition is $x=0$, so the system will always converge to that asymptotically stable point.

4.2 Controlling the circuit

Suppose now that $R_{R}$ represents a consumer, that requires a fixed amount of power $P=R_{R}I_{R}^{2}$ to be delivered to them. We would like then to control the voltage of the generator $E_{G}$, so that the state of the system will converge to $I_{R}=I_{R}^{*}:=-\sqrt{P/R_{R}}$. If we assume that a solution with $I_{R}\equiv I_{R}^{*}$ exists, then we also get $I=I_{G}\equiv-I_{R}^{*}$, $V_{1}\equiv(R_{R}+R_{L})I_{R}^{*}$, $V_{2}\equiv R_{R}I_{R}^{*}$ and $E_{G}\equiv-(R_{R}+R_{L}+R_{G})I_{R}^{*}$.

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

[TABLE]

denote the desired state. By applying the change of variables $\tilde{x}=x-x^{*}$ to (12), we get the equivalent pHDAE

[TABLE]

with the same Hamiltonian. Since our goal is having $\tilde{x}\equiv 0$ as an asymptotically stable solution, and $(J-R)x^{*}=-\sqrt{P/R_{R}}(R_{R}+R_{L}+R_{G})e_{4}=:-e_{4}u^{*}$, by construction, we write the equivalent system

[TABLE]

with $\tilde{u}=u-u^{*}$ and $\tilde{y}=y-e_{4}^{T}x^{*}=I_{G}-I_{G}^{*}$, which is again a pHDAE, but with Hamiltonian

[TABLE]

As before, this shows that choosing $\tilde{u}=0$ (i.e. $u=-u^{*}$) would make the system converge to the desired state, that will be asymptotically stable. On the other hand, this is not the only input that would satisfy this goal. If we want to speed up the convergence, we have to increase the dissipation. To do so, we can for example apply some linear feedback of the form $\tilde{u}=-\alpha\tilde{y}$ for some $\alpha>0$ (i.e. $u=u^{*}-\alpha(I_{G}-I_{G}^{*})$): in this way, the dissipation inequality can be strengthened, since

[TABLE]

From another point of view, if we apply a feedback of the form $u=u^{*}-\alpha(y-I_{G}^{*})+\hat{u}$, then (14) can be written as

[TABLE]

with $R_{\alpha}=R+\alpha e_{4}e_{4}^{T}>R$.

4.3 Numerical simulations

We present numerical simulations on the toy example (10). We choose as constants $L=2$, $C_{1}=0.01$, $C_{2}=0.02$, $R_{L}=0.1$, $R_{G}=6$ and $R_{R}=3$. This choice is not based on real world data, but it is made to reflect the usual assumption that, for transmission lines, $L\gg R_{L}\gg C$. We discretize the system with respect to time, using the midpoint rule, which is the Gauss-Legendre collocation method with $s=1$ stages and order $p=2$ for ODEs. Since (10) is a semi-explicit DAE with index 1, the chosen method will also give convergence order 2 (see [5, Theorem 5.16]).

First, we apply the time discretization to a system without control ($E_{G}=0$), starting from consistent non-zero initial values (see Fig. 2). The evolution of the state and of the Hamiltonian show the expected qualitative behaviour: after some time, they all converge to zero. Furthermore, while in the first half of the simulation the state oscillates, the Hamiltonian always decreases monotonically, following the discrete dissipation inequality.

As a second example, we apply control to the system, starting from almost-zero consistent initial values, with the goal of delivering $P=10$ power to the consumer. This corresponds to the desired state

[TABLE]

and to the control value $u^{*}=9.1\sqrt{10/3}\approx 16.614$. To simulate the fact that the change of input does not happen instantly, we actually choose as control

[TABLE]

so that the input increases smoothly and rapidly from 0 to $u^{*}$. In Fig. 3 one can see that, after some initial oscillations, the state quickly converges to the desired values, represented by the dotted lines. At the same time, the second Hamiltonian $\tilde{\mathcal{H}}$ decreases monotonically, converging to zero, while the first Hamiltonian $\mathcal{H}$, representing actual energy in the system, converges to a positive value $\mathcal{H}^{*}$.

5 Conclusions and future works

5.1 Conclusions

We have presented a new definition for port-Hamiltonian descriptor systems, generalizing the one from [1] to include a larger class of equations. Extension to weak problems and to partial differential equations is also possible. We verified that this definition satisfies several properties, in particular a dissipation inequality, invariance under a large class of variable transformations, and structure-preserving interconnection. We generalized the definition of Dirac structure and we associated one to every system included in our formulation. We extended the results from [4] to apply structure-preserving collocation schemes to port-Hamiltonian differential-algebraic equations. Finally, we illustrated a simple example from the domain of electrical circuits that can be written with our formulation, and we presented some numerical experiments.

5.2 Future Works

Ongoing work is on the analysis of a larger class of numerical methods applied to pHDAEs, including more general Runge-Kutta schemes and partitioned methods, in particular in the case of semi-explicit DAEs with index 1. Structure-preserving model reduction and space discretization for PDEs are also being considered. Future work will include further analysis of the properties of pHDAEs, in particular the extension of the results from [7] and possibly the characterization of controllability and observability of pH systems.

6 Acknowledgments

Riccardo Morandin is supported by the Deutsche Forschungsgemeinschaft (DFG) in the subproject B03, within the SFB Transregio 154. Volker Mehrmann is supported by the German Federal Ministry of Education and Research (BMBF), in the project EiFer.