Model reduction techniques for linear constant coefficient port-Hamiltonian differential-algebraic systems

TL;DR

This paper develops and adapts model reduction techniques specifically for large-scale port-Hamiltonian differential-algebraic systems, ensuring structure preservation and constraint integrity, with applications to flow and mechanical systems.

Contribution

It introduces adapted model reduction methods for port-Hamiltonian differential-algebraic systems, combining Dirac structure reduction and moment matching techniques.

Findings

Methods effectively reduce model size while preserving structure.

Validated on flow and multibody system benchmarks.

Maintains algebraic constraints during reduction.

Abstract

Port-based network modeling of multi-physics problems leads naturally to a formulation as port-Hamiltonian differential-algebraic system. In this way, the physical properties are directly encoded in the structure of the model. Since the state space dimension of such systems may be very large, in particular when the model is a space-discretized partial differential-algebraic system, in optimization and control there is a need for model reduction methods that preserve the port-Hamiltonian structure while keeping the (explicit and implicit) algebraic constraints unchanged. To combine model reduction for differential-algebraic equations with port-Hamiltonian structure preservation, we adapt two classes of techniques (reduction of the Dirac structure and moment matching) to handle port-Hamiltonian differential-algebraic equations. The performance of the methods is investigated for benchmark examples originating from semi-discretized flow problems and mechanical multibody systems.

Figures

Equations

$$E\dot{x}$$

$$y$$

$$Q^{T}E\frac{d}{dt}-Q^{T}JQ:\mathcal{X}\subset {\mathcal{C}}^1(\mathbb{I},{\mathbb{R}}^n)\to {\mathcal{C}}^0(\mathbb{I},{\mathbb{R}}^n)$$

$$W=\left[\begin{array}{cc}{Q}^{T}RQ& Q^{T}P\\ P^{T}Q& S\end{array}\right]\in {\mathbb{R}}^{(n+m)\times (n+m)}$$

$$\mathcal{H}(x)=\frac{1}{2}{x}^T{Q}^{T}Ex.$$

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

$$f_x=-E\dot{x},e_x=Qx,f_R=-\left[\begin{array}{cc}R& P\\ P^T& S\end{array}\right]e_R,e_R=\left[\begin{array}{c}Qx\\ u\end{array}\right],f_P=y,e_P=u,$$

$$((f_x,f_R,f_P),(e_x,e_R,e_P))\in \mathcal{D}\subset \mathcal{F}\times \mathcal{E}\text{ for all }t\in \mathbb{I}.$$

$$-F_x{f}_x=E_x{e}_x+F_R{f}_R+E_R{e}_R+F_P{f}_P+E_P{e}_P$$

$$\tilde{E}\dot{\tilde{x}}$$

$$y$$

$${\tilde{U}}^{T}EV=\left[\begin{array}{cc}{E}_{11}& 0\\ 0& 0\end{array}\right]$$

$$T=\left[\begin{array}{cc}I& 0\\ T_{21}& I\end{array}\right],T_{21}=-(L_{22}{)}^{-T}({\tilde{L}}_{12}{)}^T$$

$$\left[\begin{array}{cc}{E}_{11}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]$$

$$y$$

$$E_{11}{\dot{x}}_1$$

$$y$$

$${\hat{B}}_1$$

$${\hat{P}}_1$$

$$\hat{S}$$

$$\hat{N}$$

$$L_{22}{Q}_{22}{x}_2=-(L_{21}{Q}_{11}+L_{22}{Q}_{21})x_1-(B_2-P_2)u,$$

$$Q_{11}{E}_{11}^{-1}=\bar{U}\left[\begin{array}{cc}{\Sigma }_Q& 0\\ 0& 0\end{array}\right]{\bar{U}}^T,{\Sigma }_Q={\Sigma }_Q^T>0,{\bar{U}}^T={\bar{U}}^{-1}.$$

$$\left[\begin{array}{ccc}I& 0& 0\\ 0& I& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_{1,a}\\ {\dot{x}}_{1,b}\\ {\dot{x}}_2\end{array}\right]$$

$$y$$

$$\displaystyle V=\tilde{V}_{1}\left[\begin{array}[]{cc}I&0\\ 0&\tilde{V}_{2}\end{array}\right]\left[\begin{array}[]{ccc}I&0&0\\ 0&I&0\\ -A_{3}^{-1}A_{1}&0&I\end{array}\right],$$

$${\tilde{U}}_2^T({\tilde{U}}_1^{T}QV)=\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& Q_{13}\\ Q_{21}& Q_{22}& Q_{23}\\ 0& 0& Q_{33}\end{array}\right],$$

$$\left[\begin{array}{ccc}{E}_{11}& 0& 0\\ E_{21}& 0& 0\\ E_{31}& 0& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right]$$

$$y$$

$$\left[\begin{array}{cc}{E}_{11}& 0\\ E_{21}& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]$$

$$y$$

$$\left[\begin{array}{ccc}{I}_{n_1}& 0& 0\\ 0& I_{n_2}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right]$$

$$+\left[\begin{array}{c}{B}_1\\ B_2\\ B_3\end{array}\right]u,$$

$$y$$

$${\dot{x}}_i=-f_{x_i},{\nabla }_{x_i}H(x)=e_{x_i},i\in \{r,s\}$$

$$\left[\begin{array}{c}{f}_{x_{1,r}}\\ f_{x_{1,s}}\\ f_{x_2}\\ f_{x_3}\end{array}\right]=-\left[\begin{array}{cccc}I& 0& 0& 0\\ 0& I& 0& 0\\ 0& 0& I& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\dot{x}}_{1,r}\\ {\dot{x}}_{1,s}\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right],\left[\begin{array}{c}{e}_{x_{1,r}}\\ e_{x_{1,s}}\\ e_{x_2}\\ e_{x_3}\end{array}\right]=\left[\begin{array}{cccc}{Q}_{11}^{rr}& Q_{11}^{rs}& 0& 0\\ (Q_{11}^{rs}{)}^T& Q_{11}^{ss}& 0& 0\\ 0& 0& 0& 0\\ Q_{31}^r& Q_{31}^s& Q_{32}& Q_{33}\end{array}\right]\left[\begin{array}{c}{x}_{1,r}\\ x_{1,s}\\ x_2\\ x_3\end{array}\right].$$

$$V^{-1}R{V}^{-T}=\left[\begin{array}{cc}C& \hat{C}\end{array}\right]\left[\begin{array}{cc}\hat{R}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{C}^T\\ {\hat{C}}^T\end{array}\right]=C\hat{R}{C}^T,$$

Full text

Keywords. structure-preserving model reduction; index reduction; port-Hamiltonian differential-algebraic system; moment matching; effort constraint method; flow constraint method

AMS-Classification. 34H05, 41A22, 65L80, 93A15, 65F22

1. Introduction

Port-Hamiltonian differential-algebraic systems (pHDAEs) arise from port-based network modeling of multi-physics problems. For this, a physical system is decomposed into smaller subsystems that are interconnected through energy exchange. The subsystems may belong to various different physical domains, such as electrical, mechanical, or hydraulic ones. The energy-based formulation is advantageous as it brings different scales on a single level, the port-Hamiltonian character is inherited by the coupling and the physical properties (e.g., stability, passivity, energy and momentum conservation) are encoded directly in the structure of the pH model equations [4, 35]. Algebraic constraints naturally come from the interconnections in form of network conditions, such as Kirchhoff’s laws, or from constraints that are directly modeled, like position or velocity constraints in mechanical systems, or from mass balances in chemical engineering problems, see e.g., [9, 23, 29]. The state space dimension of pHDAEs can be very large, e.g., for constraint finite element models in structural mechanics [17], semi-discretized problems arising in fluid dynamics [13, 14, 21, 33] or multibody problems [25]. In this case, for optimization and control, model order reduction techniques are needed that preserve the port-Hamiltonian structure and keep the explicit and hidden algebraic constraints unchanged. To present such methods and to study their properties is the main topic of this paper, which brings together model reduction for differential-algebraic equations with structure preservation.

The properties of pHDAEs have recently been studied in [4, 35]. For systems of port-Hamiltonian ordinary differential equations (pHODEs), structure-preserving reduction methods have been developed based on ideas of tangential interpolation [18, 19], moment matching [26, 27, 39] as well as effort and flow constraint reduction methods [28]. Structure-preserving model reduction for nonlinear systems has been studied in [12] and for linear damped wave equations in [13], where particular Galerkin projections have been constructed for the pHDAEs arising in gas transport networks. Surveys of model reduction techniques for general DAEs are given in [5, 7]. A crucial step for model reduction is that the dynamic and algebraic equations are exactly identified and only the dynamic equations are reduced, otherwise the system may loose important properties, such as stability or passivity.

In this paper we generalize structure-preserving techniques that have been developed for pHODEs to pHDAEs. We focus on the effort and flow constraint model reduction methods as well as on moment matching. To do this, we follow the regularization concept of [4] to identify and decouple the algebraic constraints and the dynamical equations in a structure-preserving manner to develop the corresponding reduction methods. To illustrate the performance of the reduction methods, we apply them to benchmark problems originating from semi-discretized flow calculations and multibody systems. We discuss the advantages but also the limitations of these methods for pHDAEs.

This paper is structured as follows. The model framework of port-Hamiltonian differential-algebraic systems and a structure-preserving regularization concept is presented in Section 2. Following this, structure-preserving reduction techniques for the differential-algebraic pH systems are generalized from their ordinary differential equation counterparts in Section 3. The performance of the methods is numerically investigated on the basis of various benchmark examples in Section 4. The paper closes with a summary in Section 5.

2. Model framework of port-Hamiltonian differential-algebraic systems

In this section we review the structural properties and simplified representations of pHDAEs according to [4] and [37]. In particular, we study the decoupling of the dynamic and algebraic equations and variables, that will be used in the next section to derive structure-preserving model reduction techniques.

2.1. Port-Hamiltonian systems

Port-Hamiltonian systems can be derived in two different ways, via a formulation as descriptor systems with special structured coefficient matrices or via an energy-based formulation on top of a Dirac structure. Since each formulation is a basis for a model reduction technique, we briefly discuss their relation.

Definition 1** (pHDAE).**

A linear constant coefficient DAE system of the form

[TABLE]

with $E$, $Q$, $J$, $R\in\mathbb{R}^{n\times n}$, $B$, $P\in\mathbb{R}^{n\times m}$, $S=S^{T}$, $N=-N^{T}\in\mathbb{R}^{m\times m}$, on a compact interval $\mathbb{I}\subset\mathbb{R}$ is called a port-Hamiltonian differential-algebraic system (pHDAE) if the following properties are satisfied.

1. (1)

The differential-algebraic operator

[TABLE]

is skew-adjoint, i.e., we have that $Q^{T}J^{T}Q=-Q^{T}JQ$ and $Q^{T}E=E^{T}Q$, 2. (2)

the product $Q^{T}E=E^{T}Q$ is positive semidefinite, i.e., $Q^{T}E=E^{T}Q\geq 0$, and 3. (3)

the passivity matrix

[TABLE]

is symmetric positive semi-definite, i.e., $W=W^{T}\geq 0$.

The quadratic Hamiltonian function ${\mathcal{H}}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ of the system is given by

[TABLE]

Theorem 1**.**

Consider a pHDAE of the form (2.1). If for given input function $u$ the system has a (classical) solution $x\in\mathcal{C}^{1}(\mathbb{I},\mathbb{R}^{n})$ in $\mathbb{I}$, then

[TABLE]

Furthermore, if $W=0$, then $\frac{d}{dt}{\mathcal{H}}(x)=u^{T}y$.

Theorem 1 implies some important properties of a pHDAE. First of all, its Hamiltonian is an energy storage function, and the system is passive. A pHDAE satisfies a dissipation inequality. Furthermore, it is implicitly Lyapunov stable as $\mathcal{H}$ defines a Lyapunov function. The physical properties are encoded in the algebraic structure of the coefficient matrices and the geometric structures associated with the flow of the system. In this sense, $E^{T}Q$ is the energy matrix, $Q^{T}RQ$ is the dissipation matrix, $Q^{T}JQ$ the structure matrix describing the energy flux among the energy storage elements, $B\pm P$ are the port matrices for energy in- and output, and $S$, $N$ are the matrices associated to a direct feed-through from input $u$ to output $y$. In the case that $E=I$ is the identity matrix, the pHDAE reduces to a standard pHODE as studied in [36].

In the alternative energy-based formulation a port-Hamiltonian system is characterized by the fact that the flow and effort variables of its energy-storing port, its energy-dissipating port and its external port are linked together in a power-conserving manner by a Dirac structure. Given a finite-dimensional linear space $\mathcal{F}$ with its dual space $\mathcal{E}=\mathcal{F}^{*}$, a Dirac structure is a subset $\mathcal{D}\subset\mathcal{F}\times\mathcal{E}$ satisfying $e^{T}f=0$ for all $(f,e)\in\mathcal{D}$ and $\dim\mathcal{D}=\dim\mathcal{F}$. For a pHDAE in the form (2.1), the flow and effort variables are defined on $\mathcal{F}=\mathcal{F}_{x}\times\mathcal{F}_{R}\times\mathcal{F}_{P}$ and $\mathcal{E}=\mathcal{E}_{x}\times\mathcal{E}_{R}\times\mathcal{E}_{P}\subset\mathbb{R}^{n}\times\mathbb{R}^{n+m}\times\mathbb{R}^{m}$, respectively. They are given by

[TABLE]

The variables $(f_{x},e_{x})\in\mathcal{F}_{x}\times\mathcal{E}_{x}$ of the energy-storing port are related to the evolution of the state and the Hamiltonian $\mathcal{H}$. If $E=I$, then the constitutive relations read as $(f_{x},e_{x})=(-\dot{x},\nabla_{x}{\mathcal{H}}(x))$. The port variables of the energy-dissipating elements satisfy a resistive relation, $e_{R}^{T}f_{R}\leq 0$ for all $(f_{R},e_{R})\in\mathcal{R}\subset\mathcal{F}_{R}\times\mathcal{E}_{R}$, which is encoded in the stated positive semidefinite matrix $R$. The external port variables $(f_{P},e_{P})\in\mathcal{F}_{P}\times\mathcal{E}_{P}$ correspond to the out- and inputs of the system. The energy-conservation property follows directly from Theorem 1. Based on the notion of the Dirac structure, the pH system possesses a DAE representation, see [36, 37], i.e.,

[TABLE]

with matrices $F_{x},E_{x}\in\mathbb{R}^{q\times n},$ $F_{R},E_{R}\in\mathbb{R}^{q\times(n+m)}$ and $F_{P},E_{P}\in\mathbb{R}^{q\times m}$ where $q=n+(n+m)+m$ and $\sum_{i=x,R,P}E_{i}F_{i}^{T}+F_{i}E_{i}^{T}=0$.

2.2. Structure-preserving regularization

A pHDAE system typically contains explicit as well as implicit (hidden) constraints. Since in model reduction all constraints need to be kept unchanged in order not to destroy crucial properties, we need to identify all constraints. If the differentiation-index is larger than one, then an index reduction, e.g., via derivative arrays or minimal extension, should be performed, see [23] for general DAEs. For pHDAEs this index reduction has to be performed in a structure-preserving way, see [4]. It has been shown in [24] for the linear constant coefficient case and in [32] for the linear time-varying case that the differentiation-index will be at most two, i.e., in simple terms at most the second derivative of the input function $u$ is required to transform the system into a pHODE. In contrast to the numerical solution of pHDAEs/pHODEs via time-integration methods, which is still partially an open problem, see, e.g. [22], in the context of model reduction also a structure-preserving decoupling of the dynamic and algebraic variables should be performed. This may be a very critical step for linear time-varying or nonlinear systems, since it may require time-varying changes of variables, with all its difficulties, in particular of having to provide derivatives of the transformation functions [23]. But even in the case of constant coefficients, changes of variables with ill-conditioned transformation matrices may have to be handled.

In the following we discuss the structure-preserving regularization for linear constant coefficient pHDAEs of differentiation-index one or two. We refer to [4] for a detailed analysis of the regularization concept for linear time-varying systems. The concept is particularly based on the fact that the port-Hamiltonian structure and the associated Hamiltonian are preserved under basis change and scaling with invertible matrices (cf. Lemma 1).

Lemma 1**.**

Consider a pHDAE of the form (2.1) with Hamiltonian $\mathcal{H}$ (2.2). Let $U$, $V\in\mathbb{R}^{n\times n}$ be invertible. Then the transformed system

[TABLE]

with $\tilde{E}=U^{T}EV$, $\tilde{J}=U^{T}JU$, $\tilde{R}=U^{T}RU$, $\tilde{B}=U^{T}B$, $\tilde{P}=U^{T}P$, as well as $\tilde{Q}=U^{-1}QV$ and $\tilde{x}=V^{-1}x$ is still a pHDAE with the same Hamiltonian $\tilde{\mathcal{H}}(\tilde{x})=\frac{1}{2}\tilde{x}^{T}\tilde{Q}^{T}\tilde{E}\tilde{x}={\mathcal{H}}(x)$.

Decoupling of pHDAE of index at most one

For the decoupling of a pHDAE (2.1) of index at most one, two orthogonal matrices $\tilde{U}$, $V\in\mathbb{R}^{n\times n}$ are determined (e.g., via a singular decomposition) such that

[TABLE]

with $E_{11}$ invertible. We set $L=J-R$ and apply the transformation induced by $\tilde{U}$ and $V$ to (2.1) (cf. Lemma 1). Partitioning as in (2.4) yields block-structured system matrices whose blocks we denote by $\tilde{\quad}$ in case that they change in the decoupling procedure. As $Q^{T}E$ is real symmetric, we have $Q^{T}_{11}E_{11}=E^{T}_{11}Q_{11}$ and $Q_{12}=0$. Furthermore, as the system is of differentiation-index at most one, the block matrix $L_{22}Q_{22}$ is either not present – in case of an implicitly defined pHODE – or it is invertible, i.e., $L_{22}$ and $Q_{22}$ both are invertible. Setting $U=\tilde{U}T$ with

[TABLE]

and transforming (2.1) with $U$ and $V$ as in Lemma 1 yields the block-structured pHDAE

[TABLE]

where $L_{ij}=J_{ij}-R_{ij}$.

Theorem 2** (Decoupled pHDAE).**

Suppose that the pHDAE (2.1) is of differentiation-index at most one. Let the system be transformed to the form (2.5) via $U$ and $V$ and define $V^{-1}x=[x_{1}^{T}\,\,x_{2}^{T}]^{T}$. Then for any input $u$ and initial condition $x_{1}(t_{0})=x_{1,0}$, the output $y$ and the state $x_{1}$ of (2.5) are given by the implicit pHODE

[TABLE]

with Hamiltonian $\hat{\mathcal{H}}(x_{1})=\frac{1}{2}x_{1}^{T}Q^{T}_{11}E_{11}x_{1}={\mathcal{H}}(x)$, and coefficients

[TABLE]

The state $x_{2}$ is uniquely determined by the explicit algebraic constraint

[TABLE]

which implies a consistency constraint for the respective initial condition.

Typically the original pencil $(sE-LQ)$, $s\in\mathbb{C}$, is regular, i.e., its determinant is not identically zero, which means that a unique solution exist for every sufficiently smooth input function $u$ and every consistent initial condition. If this is not the case, then a complex regularization procedure can be performed, which consists of transformations, feedbacks and renaming of variables, see [11]. Since this procedure is not yet available for pHDAEs, in the following we assume that $(sE-LQ)$ is regular. Then the pencil $(sE-Q)$ is regular as well as shown in [24]. In this case (2.5) can be decoupled even further by identifying the zero eigenvalues of the system. For this, a change of basis is applied to the dynamic state $x_{1}$. From $E^{T}Q=Q^{T}E\geq 0$ it follows that the block matrix $Q_{11}E_{11}^{-1}$ is symmetric positive semidefinite and allows for an ordered Schur decomposition that can be obtained from the generalized singular value decomposition [16]. Even though we will not carry out this transformation explicitly, it follows that the system can be transformed as

[TABLE]

Introducing $\bar{U}^{T}E_{11}x_{1}=[x_{1,a}^{T}\,\,x_{1,b}^{T}]^{T}$ yields then a pHDAE of the form

[TABLE]

Remark 1**.**

In the decoupled form (2.6), some further transformations can be applied to achieve $Q_{21}^{a}=0$. However, since the inverse of $Q_{22}$ is involved in this transformation, we stay with the form (2.5) ((2.6), respectively).**

Decoupling for pHDAE of index two

To decouple the algebraic and differential variables for a pHDAE of differentiation-index two we make use of the index-reduction procedure developed in [4] for the linear time-varying case. Assume that the state equation with $u=0$ forms a DAE of index two. It has been shown in [10] that the extra constraint equations (hidden constraints) that arise from derivatives are uncontrollable, because otherwise the index reduction could have been done via feedback. This means that these hidden constraint equations are of the form $\hat{A}x=0$. We just add these constraint equations to our original pHDAE (2.1) and obtain an overdetermined DAE system, see also [23]. Then we perform a singular value decomposition of $E$ by means of orthogonal matrices $\tilde{U}_{1}$, $\tilde{V}_{1}\in\mathbb{R}^{n\times n}$ as in (2.4) and partition the matrix associated to the extra constraints accordingly ${\hat{A}}\tilde{V}_{1}=[A_{1}\,\,\tilde{A}_{2}]$. Note that we denote blocks by $\tilde{\quad}$ in case that they change in the reduction procedure. The equations $\hat{A}x=0$ include all the constraints that are needed for index reduction. Since $\tilde{E}_{11}$ is invertible, these extra equations must arise from the full row-rank part of $\tilde{A}_{2}$. In the following we assume w.l.o.g. that $\tilde{A}_{2}$ has full row-rank. This can be always achieved by transforming and omitting hidden constraint equations that do not contribute to the index reduction, see [4]. Then there exists an orthogonal matrix $\tilde{V}_{2}$ such that $\tilde{A}_{2}\tilde{V}_{2}=[0\,\,A_{3}]$ with $A_{3}$ invertible. Introducing $[x_{1}^{T}\,\,x_{2}^{T}\,\,x_{3}^{T}]^{T}=V^{-1}x$ with

[TABLE]

the hidden constraint equations become $A_{3}x_{3}=0$, consequently it follows that $x_{3}=0$. In addition we use an orthogonal matrix $\tilde{U}_{2}$ such that

[TABLE]

transforming (2.1) with $U=\tilde{U}_{1}\tilde{U}_{2}$ and $V$ as in Lemma 1 yields then the following block-structured pHDAE

[TABLE]

together with the constraint $x_{3}=0$, where $L_{ij}=J_{ij}-R_{ij}$. Since the constraint does not change the solution, the subsystem given by the first two block rows is a pHDAE of differentiation-index at most one, i.e.,

[TABLE]

This system (2.7) can then be decoupled as described in the previous paragraph.

3. Model reduction techniques

In this section we present two different classes of model order reduction methods for pHDAEs. The first class is based on the reduction of the underlying Dirac structure and the associated power conservation, whereas the second one, the well-known moment matching, aims at the approximation of the transfer function. To derive the reduction methods, we assume that the pHDAE is in analogous form as (2.6), but for convenience we neglect the feed-through terms, i.e., we assume that $S=N=0$ and then as consequence $P=0$, since $W\geq 0$. However, generalizations to systems with feed-through term are straightforward.

In the following we use an adapted notation, considering the pHDAE

[TABLE]

for the states $x_{i}\in\mathbb{R}^{n_{i}}$, $\sum_{i=1}^{3}n_{i}=n$, with the accordingly block-structured coefficient matrices $E$, $Q$, $Q^{T}JQ=-Q^{T}J^{T}Q$, $Q^{T}RQ=Q^{T}R^{T}Q\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times p}$, where $(J_{j3}-R_{j3})=0$ for $j=1,2$, $Q_{11}>0$ and $Q_{33}$ is nonsingular.

3.1. Power conservation based model order reduction

The power conservation based methods were originally developed for standard pHODEs in [28], i.e., systems of the form (2.1) with $E=I$ being the identity. If for such a system a splitting of the dynamic state as $x=[x_{r}^{T}\,\,x_{s}^{T}]^{T}$ exists, with $x_{r}\in\mathbb{R}^{r}$ and $x_{s}\in\mathbb{R}^{n-r}$, where $x_{s}$ does not contribute much to the input-output behavior of the system, then the general idea is to cut the interconnection between the part of the energy storage port belonging to $x_{s}$ and the Dirac structure, such that no power is transferred. Then the power is exclusively exchanged via the energy storage of $x_{r}$, which will act as reduced state variable, whereas $x_{s}$ will be skipped. The constitutive relations become

[TABLE]

with $e_{x}=[e_{x_{r}}^{T}\,\,e_{x_{s}}^{T}]^{T}$ and $f_{x}=[f_{x_{r}}^{T}\,\,f_{x_{s}}^{T}]^{T}$. To cut the interconnection, one forces one of the power products $(\nabla_{x_{s}}H)^{T}\dot{x}_{s}$ or $e_{x_{s}}^{T}f_{x_{s}}$ to be zero. Choosing $\nabla_{x_{s}}H(x)=e_{x_{s}}=0$ yields the effort constraint reduction method (ECRM), while setting $\dot{x}_{s}=-f_{x_{s}}=0$ results in the flow constraint reduction method (FCRM).

We now adopt these reduction procedures of [28] to reduce pHDAEs in the form (2.3). Proceeding from (3), let $\hat{V}^{-1}x_{1}=[x_{1,r}^{T}\,\,x_{1,s}^{T}]^{T}$, $x_{1,r}\in\mathbb{R}^{r}$, $x_{1,s}\in\mathbb{R}^{n_{1}-r}$, be an appropriate splitting of the dynamic part of the state variable with respect to its relevance for the input-output behavior. In the resulting system transformed by means of $V^{-T}$ and $V$ as in Lemma 1 with a block diagonal matrix $V=\mathrm{diag}(\hat{V},I_{n_{2}},I_{n_{3}})\in\mathbb{R}^{n\times n}$, we denote the state by $x=[{x^{T}_{1,r}}\,\,{x^{T}_{1,s}}\,\,x_{2}^{T}\,\,x_{3}^{T}]^{T}$. By definition, the flow and effort variables of the energy-storing and energy-dissipating ports inherit the partitioning, and the coefficient matrices are structured accordingly. The constitutive relations then read $f_{x}=-(V^{-1}EV)\dot{x}$ and $e_{x}=(V^{T}QV)x$, with

[TABLE]

For the model reduction we have to open the resistive port. The transformed symmetric positive semi-definite dissipation matrix $V^{-1}RV^{-T}\in\mathbb{R}^{n\times n}$ admits an ordered Schur decomposition

[TABLE]

with $0<\hat{R}=\hat{R}^{T}\in\mathbb{R}^{\ell\times\ell}$ and $C\in\mathbb{R}^{n\times\ell}$, with $\ell$ being the number of energy-dissipating elements. Plugging (3.3) into the transformed system and introducing the associated flow and effort variables accordingly, i.e., $f_{R}=-\hat{R}e_{R}$ and $e_{R}=C^{T}(V^{T}QV)x=C^{T}e_{x}$, yields a pHDAE with opened resistive port. Inserting the constitutive relations (3.2) and introducing the external port variables $(f_{P},e_{P})=(y,u)$, where $y=(V^{-1}B)^{T}(V^{T}QV)x=(V^{-1}B)^{T}e_{x}$, we obtain a new representation of (3) as

[TABLE]

In the Effort Constraint Reduction Method (ECRM) the energy transfer between the energy-storing elements and the Dirac structure is cut by setting $e_{x_{1,s}}=0$. Here, the Hamiltonian is considered only weakly influenced by $x_{1,s}$. The relation $x_{1,s}=-(Q_{11}^{ss})^{-1}(Q_{11}^{rs})^{T}x_{1,r}$ follows in a straightforward way from (3.2) as $Q_{11}^{ss}=(Q_{11}^{ss})^{T}>0$. The reduced Dirac structure is obtained by multiplying (3.4) from the left with any matrix $D^{\mathrm{ec}}$ of maximal rank satisfying $D^{\mathrm{ec}}F_{x_{1,s}}=0$, such as,

[TABLE]

Rewriting the resulting system again as DAE and closing the resistive port yields the reduced model (3).

Theorem 3** (Reduced model by ECRM).**

Consider a pHDAE of the form (3) with its representation (3.4). Then the reduced model obtained by ECRM with reduced state $x^{\mathrm{ec}}=[x^{T}_{1,r}\,\,x^{T}_{2}\,\,x^{T}_{3}]^{T}\in\mathbb{R}^{(r+{n_{2}}+{n_{3}})}$, $r\ll n_{1}$ is a pHDAE given by

[TABLE]

with $J^{\mathrm{ec}}=-(J^{\mathrm{ec}})^{T}$, $R^{\mathrm{ec}}=(R^{\mathrm{ec}})^{T}\geq 0$, where ${\hat{Q}}_{11}=Q_{11}^{rr}-Q_{11}^{rs}(Q_{11}^{ss})^{-1}(Q_{11}^{rs})^{T}$ and ${\hat{Q}}_{31}=Q_{31}^{r}-Q_{31}^{s}(Q_{11}^{ss})^{-1}(Q_{11}^{rs})^{T}$.

Proof.

In the described reduction the port-Hamiltonian properties are inherited. The skew-symmetry of $J^{\mathrm{ec}}$ follows trivially, since

[TABLE]

The symmetry and positive semi-definiteness

[TABLE]

follows from ${\hat{Q}}_{11}={\hat{Q}}_{11}^{T}>0$ as it is constructed from the Schur complement of a positive definite matrix. Finally,

[TABLE]

where the component $x_{1,s}$ is projected out by means of $D_{b}$. ∎

In the Flow Constraint Reduction Method (FCRM), the energy transfer between the energy-storing elements and the Dirac structure is cut by setting $f_{x_{1,s}}=0$ and $\dot{x}_{1,s}=0$. Thus, $x_{1,s}$ is constant and can particularly be chosen as $x_{1,s}=0$. The reduced Dirac structure is obtained by multiplying (3.4) from the left with any matrix $D^{\mathrm{fc}}$ of maximal rank satisfying $D^{\mathrm{fc}}E_{x_{1,s}}=0$, e.g., with

[TABLE]

in the case that $(J_{11}^{ss})^{-1}$ exists. Rewriting the system as a DAE, analogously to ECRM, yields the reduced model (3.6).

Theorem 4** (Reduced model by FCRM).**

Consider a pHDAE of the form (3) with its general representation (3.4) and suppose that $J_{11}^{ss}$ in (3.4) is invertible. Then the reduced order model obtained by FCRM is port-Hamiltonian and, for state $x^{\mathrm{fc}}=[x^{T}_{1,r}\,\,x^{T}_{2}\,\,x^{T}_{3}]^{T}\in\mathbb{R}^{(r+{n_{2}}+{n_{3}})}$, $r\ll n_{1}$, given by

[TABLE]

with

[TABLE]

as well as $Z=\hat{R}(I-\mathcal{D}\hat{R})^{-1}$ with its symmetric and skew-symmetric parts, $Z_{R}=(Z+Z^{T})/2$ and $Z_{J}=(Z-Z^{T})/2$.

Proof.

In the described reduction the port-Hamiltonian properties are inherited. It follows trivially that $J^{\mathrm{fc}}$, $N^{\mathrm{fc}}$ are skew-symmetric and $R^{\mathrm{fc}}$, $S^{\mathrm{fc}}$ are symmetric by construction. From $Q_{11}^{rr}=(Q_{11}^{rr})^{T}>0$ it results that $(E^{\mathrm{fc}})^{T}Q^{\mathrm{fc}}=(Q^{\mathrm{fc}})^{T}E^{\mathrm{fc}}\geq 0$. Finally,

[TABLE]

holds, since $Z_{R}$ is positive (semi-)definite, see [28]. ∎

The reduced models obtained by ECRM and FCRM have similarities but also show crucial differences. Obviously, $E^{\mathrm{ec}}=E^{\mathrm{fc}}$, implying the same size of the reduced states. The energy matrices $Q^{\mathrm{ec}}$, $Q^{\mathrm{fc}}$ only differ in the first column. In particular, ECRM generates an additive term in the matrix $\hat{Q}_{11}$, and also in the matrix $\hat{Q}_{31}$ associated to the algebraic constraints due to the elimination of $x_{1,s}$, cf. (3). Clearly, the models differ in the feed-through term which is only present in FCRM (see $S^{\mathrm{fc}}$, $N^{\mathrm{fc}}$ and $P^{\mathrm{fc}}$ in (3.6)), even though the original system (3) did not have such a term. It is further important to note that the construction of ECRM is always applicable, whereas in the presented form FCRM requires the skew-symmetric matrix $J^{ss}_{11}$ to be invertible, which is, e.g., impossible if the size $(n_{1}-r)$ is odd. In the case that $J^{ss}_{11}$ is singular, the procedure has to be modified, but we do not present this modification here because it gets rather technical.

Remark 2**.**

In the presented power conservation based methods the general representation (3.4) for the pHDAE differs from the one for a standard pHODE [28] by the two additional block rows and columns for the equations associated with the kernel of the energy matrix and with the algebraic constraints. Hence, ECRM yields a reduced model (3) with equivalent block matrices in the dynamic part. However, although the reduction is only applied to the dynamic state, additional terms in $\hat{Q}_{31}$ associated to the algebraic constraints are generated. Also the reduced model (3.6) by FCRM contains additional blocks, i.e., $\mathcal{J}_{1j}$, $\mathcal{J}_{2j}$, $\mathcal{J}_{33}$, ${\mathcal{B}}_{j}$ and ${\mathcal{C}}_{j}$, $j\in\{2,3\}$, which arise from the additional equations. The coefficients $\mathcal{J}_{11}$, ${\mathcal{B}}_{1}$, ${\mathcal{C}}_{1}$, $\mathcal{G}$, $\mathcal{D}$ and $\mathcal{N}$ are analogous to their counterparts for a pHODE.

3.2. Moment Matching

The model reduction procedure of moment matching (MM) derives a reduced order model by means of a Galerkin projection in such a way that the leading coefficients of the series expansion of its transfer function (its moments) match those of the full order system, i.e., $G(s)=B^{T}(sE-(J-R)Q)^{-1}B=\sum_{j=0}^{\infty}m_{j}(s_{0}-s)^{j}$ with moments $m_{j}$ associated with a given shift parameter $s_{0}$. For details of MM for DAEs we refer to [15] for $s_{0}\in\mathbb{C}$ and to [6] for $s_{0}=\infty$. To apply these techniques in a structure-preserving way to the pHDAE (3), the symmetric positive definite energy matrix block $Q_{11}$ associated to the dynamic state is first transformed to become an identity, as done in the works on MM for pHODEs [26, 27]. Performing a Cholesky factorization $Q_{11}=KK^{T}$ and transforming appropriately, by Lemma 1 with $U=\mathrm{diag}(K,I_{n_{2}},I_{n_{3}})$, $V=U^{-T}\in\mathbb{R}^{n\times n}$, the resulting system is still port-Hamiltonian. Then a Galerkin projection matrix for the dynamic part $V_{r}\in\mathbb{R}^{n_{1}\times r}$, $r\ll n_{1}$, can be computed, e.g., by an Arnoldi method [15], such that $V_{r}^{T}V_{r}=I_{r}$ and its columns span a Krylov space of associated to the system shifted by $s_{0}$ [30]. Applying finally the Galerkin projection with $V^{\mathrm{m}}_{r}=\mathrm{diag}(V_{r},I_{n_{2}},I_{n_{3}})\in\mathbb{R}^{n\times(r+n_{2}+n_{3})}$ yields the reduced model.

Theorem 5** (Reduced model by MM).**

Consider a pHDAE of the form (3) with energy matrix block $Q_{11}=KK^{T}$. Let the projection matrix $V_{r}\in\mathbb{R}^{n_{1}\times r}$ be computed by the Arnoldi method. Then the reduced system is port-Hamiltonian, matches the first $r$ moments of the full order system and is for the state $x^{\mathrm{m}}=[x_{1,r}^{T}\,\,x_{2}^{T}\,\,x_{3}^{T}]^{T}\in\mathbb{R}^{(r+n_{2}+n_{3})}$, $r\ll n_{1}$, with $x_{1,r}=V_{r}^{T}K^{T}x_{1}$ given by

[TABLE]

with $L_{ij}=J_{ij}-R_{ij}$, $i,j=1,2,3$.

Proof.

The port-Hamiltonian structure is trivially preserved by the $Q_{11}$-associated transformation and the subsequent Galerkin projection.. The matching of the first $r$ moments is proved in [15] for $s_{0}\in\mathbb{C}$ and in [6] for $s_{0}\in\infty$. ∎

4. Numerical results

In this section we investigate the performance of the model reduction methods, using benchmark examples from the literature, see, e.g., [3, 8, 13, 20, 21, 25] or [33]. In order to perform model reduction, it is essential to identify all constraints arising from the physics of the problem as discussed in Section 2.2. In many applications this can be done directly by exploiting the structure of the equations coming from the physical properties.

Considering the transfer function $G$ of the original full order system, we study the approximation quality of the various reduced order systems by comparing the relative errors $(G-G_{r})/G$ with $G_{r}\in\{G^{\mathrm{ec}},G^{\mathrm{fc}},G^{\mathrm{m}}\}$ being the transfer function of the reduced system. Usually this is done in the $L_{2}$-norm in the state-space formulation or in the $H_{\infty}$-norm in the frequency domain. Since the latter norm is only defined for pHDAEs of index at most one, for pHDAEs of higher index we present the errors for the index reduced system. The transfer function of a pHDAE of the form (2.1) is given by

[TABLE]

For pHDAEs of index at most one it is either a proper rational function or the sum of a proper rational function with a term that is constant in $s$.

The numerical results have been computed with Matlab 2017a on a Linux 64-Bit machine with an Intel® Core™ i7-6700 processor. In the context of computing the error norms, see [1, 34], it is necessary to solve Lyapunov equations, for which we have used the M.M.E.S.S. Toolbox [31].

4.1. Flow problems

Consider an instationary incompressible fluid flow, prescribed in terms of velocity $v:\Omega\times[0,T]\rightarrow\mathbb{R}^{2}$ and pressure $p:\Omega\times[0,T]\rightarrow\mathbb{R}$ on the spatial domain $\Omega=(0,1)^{2}$ with boundary $\partial\Omega$ for the time period $[0,T]$, that is driven by external forces $f:\Omega\times[0,T]\rightarrow\mathbb{R}^{2}$ and dynamic viscosity $\nu>0$. The subsequently described models of partial differential equations (Stokes as well as Oseen equations) are closed by non-slip boundary conditions and appropriate initial conditions $v^{0}$. Spatial discretization by a finite difference method on a uniform staggered grid yield pHDAEs for the state $x=[v_{h}^{T}\,\,p_{h}^{T}]^{T}$ with the semi-discretized values of the velocities $v_{h}(t)\in\mathbb{R}^{n_{v}}$ and pressures $p_{h}(t)\in\mathbb{R}^{n_{p}}$, $t\in[0,T]$ ($n_{v},n_{p}\in\mathbb{N}$, cf. Remark 3).

Stokes equations

A laminar flow can be modeled by the linear Stokes equations,

[TABLE]

A spatially semi-discretized differential-algebraic system for $x=[v_{h}^{T}\,\,p_{h}^{T}]^{T}$ completed with an appropriate output equation is given by

[TABLE]

with the symmetric negative definite discrete Laplace operator $\mathcal{L}\in\mathbb{R}^{n_{v}\times n_{v}}$, as well as the discrete divergence $D$ and gradient operators $D^{T}\in\mathbb{R}^{n_{v}\times n_{p}}$, $n_{p}<n_{v}$. The operator $D$ usually has full row rank if the freedom in the pressure (which only occurs in differentiated form in the system) is removed, see Remark 3 on page 3. The initial conditions are $v_{h}(0)=v_{h}^{0}$ and consistently $p_{h}(0)=p_{h}^{0}$. The input $u$ with input matrix $F\in\mathbb{R}^{n_{v}\times m}$ results from the external forces. The output equation is supplemented accordingly regarding the port-Hamiltonian form (2.1). System (4.1) is obviously port-Hamiltonian with $Q=I$, $P=0$ and $S=N=0$, as $E^{T}Q=Q^{T}E\geq 0$, $J=-J^{T}$ and $R=R^{T}\geq 0$ holds.

For the index reduction of the pHDAE (4.1) (which is of differentiation-index two) we do not need the whole derivative array; instead we can easily identify the equations that have to be differentiated from the special structure of the system by performing, e.g., a singular value decomposition

[TABLE]

with orthogonal matrices $U\in\mathbb{R}^{n_{v}\times n_{v}}$, $V\in\mathbb{R}^{n_{p}\times n_{p}}$ and singular values $\sigma_{i}>0$, $i=1,\dots,n_{p}$. Setting ${Z}=V\Sigma_{D}$, performing an equivalence transformation with $U$ and splitting the state variable accordingly into three parts, we get the system

[TABLE]

Obviously $x_{1}=0$, as the last equation is $0={Z}x_{1}$ with $Z$ invertible. This is the equation that has to be differentiated and inserted into the first equation to derive the second (hidden) algebraic constraint

[TABLE]

as well as a consistency condition for the initial value which relates the initial condition for $u$ and $x_{2}$ to that for $x_{3}$. The second equation yields the underlying ODE of the system for the variable $x_{2}:[0,T]\rightarrow\mathbb{R}^{(n_{v}-n_{p})}$ to be reduced and the output equation

[TABLE]

Note that this equation can be interpreted as the discretized heat equation in the set of divergence-free velocities [14].

In system (4.2) the skew-symmetric interconnection matrix is zero. Thus, FCRM cannot be applied directly as discussed in the previous section. Concerning ECRM, the splitting of the dynamic state is provided by Lyapunov balancing. The transformation matrices are computed by the Square-Root Algorithm [2]. Note that in this case the reduced model by ECRM is equivalent to the one obtained by using Balanced Truncation due to the symmetry of the system matrix and the relation of the input and output matrices.

Oseen equations

A flow model that is closer to the nonlinear Navier-Stokes equations is given by the Oseen equations

[TABLE]

The Oseen equations differ from the Stokes equations by the additional convective term with driving velocity $a:\Omega\rightarrow\mathbb{R}^{2}$. The associated spatially semi-discretized pHDAE for $x=[v_{h}^{T}\,\,p_{h}^{T}]^{T}$ is given by

[TABLE]

where the appropriately discretized convective term is decomposed into its skew-symmetric part $A$ and its symmetric part. The last forms, together with the discrete Laplacian, the symmetric negative definite operator $\tilde{\mathcal{L}}$.

Analogously to the index reduction performed for the Stokes equations, we obtain $x_{1}=0$ as well as $Z^{T}x_{3}=(A_{12}+\mathcal{\tilde{L}}_{12})x_{2}+B_{1}u$. The underlying ODE and the output equation are given by

[TABLE]

In system (4.3) the skew-symmetric interconnection matrix is prescribed by $A_{22}$ which, depending on the discretization scheme, may or may not be invertible. If it is invertible, then in contrast to the Stokes problem, FCRM can be applied for model order reduction. For ECRM and FCRM the balancing transformations may be computed via the Balancing Free Square-Root Algorithm [38] to preserve the structure and to avoid creating an additional energy matrix.

Remark 3**.**

In the numerical examples the flow domain $\Omega=(0,1)^{2}$ is partitioned into uniform quadratic cells of edge length $h=1/M$, $M\in\mathbb{N}$. We use finite differences on a staggered grid where the velocity components $v=[v^{\xi}\,\,v^{\eta}]^{T}$ are evaluated at the center of the cell faces to which they are normal and the pressure is taken at the cell centers. This procedure provides small discretization stencils and ensures numerical stability. The unknowns are

[TABLE]

for $t\in[0,T]$. Since the pressure is non-unique in the flow equations, we fix without loss of generality the value $p_{M-0.5,M-0.5}(t)=0$ for the numerical solution and discard this quantity from the unknowns. The unknowns are ordered row-wise (and for the velocity component-wise) in the vectors $v_{h}(t)\in\mathbb{R}^{n_{v}}$ and $p_{h}(t)\in\mathbb{R}^{n_{p}}$, where $n_{v}=2M(M-1)$ and $n_{p}=M^{2}-1$, yielding the state $x=[v_{h}^{T}\,\,p_{h}^{T}]^{T}:[0,T]\rightarrow\mathbb{R}^{n}$, $n=3M^{2}-2M-1$, of the pHDAE systems. **

We illustrate the model reduction techniques by comparing the approximation quality of their reduced transfer functions for both flow problems. The presented results are given for the example flow setup, where we use a random normally distributed input matrix $B\in(\mathcal{N}(0,10^{2}))^{n_{v}\times 1}$, the dynamic viscosity $\nu=1$ and, in case of the Oseen equations, the constant convective velocity $a\equiv[1\,\,1]^{T}$. Applying a spatial resolution of $M=23$, the states have a size of $n=1540$ for the full order pHDAE models and of $n_{2}=484$ for the underlying ODEs to be reduced. Figure 4.1 shows the relative errors in the spectral norm for the reduced models of size $r=16$ and in the $H_{\infty}$-norm for $r\in[2,20]$. The results of MM are similar for the Stokes and Oseen equations. As expected, MM at $s_{0}=\infty$ and at $s_{0}=0$ yields only negligibly small errors for high or low frequencies, respectively. In the $H_{\infty}$-norm MM at $s_{0}=0$ performs better than MM at $s_{0}=\infty$.

For the Stokes problem the error of ECRM is small for all frequencies, oscillating around $\mathcal{O}(10^{-13})$, whereas for the Oseen problem it decreases from $\mathcal{O}(10^{-8})$ to $\mathcal{O}(10^{-13})$ for increasing frequencies. In the $H_{\infty}$-norm ECRM outperforms MM for both flow problems. The same error trends can be observed in the $H_{2}$-norm. Consequently, ECRM yields better reduced models than the moment matching methods globally, as even for low and high frequencies the relative errors only differ slightly from the respective errors of MM. FCRM is not applicable to the Stokes problem and to the Oseen problem only if the invertibility requirement on the skew-symmetric interconnection submatrix $J^{ss}_{11}$ is satisfied. This implies in this example setup that the size $r$ of the reduced model has to be even. Here, FCRM yields smaller errors for low frequencies than ECRM, but for higher frequencies the error increases monotonically because of the additional feed-through term in the reduced system, see also Figure 4.2. However, in the $H_{\infty}$-norm FCRM even outperforms ECRM, whereas the $H_{2}$-norm is unbounded for systems with nonzero feed-through terms.

4.2. Damped mass-spring system

The holonomically constrained damped mass-spring system is a multibody problem that describes the one-dimensional dynamics of $g$ connected mass points in terms of their positions $p:[0,T]\rightarrow\mathbb{R}^{g}$, velocities $v:[0,T]\rightarrow\mathbb{R}^{g}$ and a Lagrange multiplier $\lambda:[0,T]\rightarrow\mathbb{R}$, see Figure 4.3. In the chain of mass points the $i$th mass of weight $m_{i}$ is connected to the $(i+1)st$ mass by a spring and a damper with constants $k_{i}$ and $d_{i}$ and also to the ground by a spring and a damper with the constants $\kappa_{i}$ and $\delta_{i}$, respectively, where $m_{i}$, $k_{i}$, $d_{i}$, $\kappa_{i}$, $\delta_{i}>0$. Furthermore, the first and the last mass points are connected by a rigid bar. The vibrations are driven by an external force $u:[0,T]\rightarrow\mathbb{R}$ (control input) acting on the first mass point. The resulting DAE system of size $n=2g+1$ is not port-Hamiltonian and has differentiation-index three. In first order formulation it is given by

[TABLE]

with mass matrix $M=\mathrm{diag}(m_{1},\dots,m_{g})$, tridiagonal stiffness and damping matrices $K$, $D\in\mathbb{R}^{g\times g}$, constraint matrix $G=[1\,\,0\dots 0\,\,-1]\in\mathbb{R}^{1\times g}$ and input matrix $F=e_{1}\in\mathbb{R}^{g\times 1}$. Assuming $K$ and $D$ to be symmetric negative semi-definite, the multibody problem can be formulated as a pHDAE of differentiation-index two by replacing the algebraic constraint $Gp=0$ by its first derivative $Gv=0$ and adding an appropriate output equation,

[TABLE]

Obviously, $J=-J^{T}$, $R=R^{T}\geq 0$, $E^{T}Q=Q^{T}E\geq 0$ and $Q^{T}RQ\geq 0$ hold. But note that $(J-R)Q$ is singular such that for $s=0$ the matrix $(sE-(J-R)Q)$ is not invertible.

The structure of the equations simplify a further index reduction for (4.4) analogous to that for the flow problems. The singular value decomposition of $G^{T}$, i.e., $GV=[Z\,\,0]$ with orthogonal matrix $V=[V_{1}\,\,V_{2}]$, $V_{1}\in\mathbb{R}^{g\times 1}$, $V_{2}\in\mathbb{R}^{g\times(g-1)}$ as well as $Z\in\mathbb{R}^{1\times 1}$ invertible, yields

[TABLE]

The last equation $Zv_{1}=0$ implies that $v_{1}=0$. Differentiating this equation and inserting it into the second equation yields the hidden constraint for the Lagrange multiplier

[TABLE]

which also imposes a consistency condition for the initial value. The underlying ODE of size $n_{1}=2g-1$ together with the output equation are given by

[TABLE]

We present numerical results for an example setup, where the parameters of the mass-spring system are set to be $m_{i}=100$ for $i=1,...,g$ and $k_{i}=\kappa_{j}=2$, $d_{i}=\delta_{j}=5$ for $i=1,...,g-1$, $j=2,...,g-1$ as well as $\kappa_{1}=\kappa_{g}=4$, and $\delta_{1}=\delta_{g}=10$. Choosing $g=6000$ yields a state of size $n=12001$ for the original DAE and of $n_{1}=11999$ for the underlying ODE to be reduced. Since the dimension is odd, the skew-symmetric interconnection matrix is singular and also the relevant submatrix for FCRM is singular as well. But MM at almost all $s_{0}\in\mathbb{C}\backslash\{0\}$ and ECRM can be applied. Figure 4.4 shows the respective relative errors of the transfer functions in spectral norm for the reduced size $r=10$ and in the $H_{\infty}$-norm for $r\in[2,20]$. As expected, MM at $s_{0}=\infty$ and $s_{0}=10^{-10}$ yields outstanding approximations (errors of order $\mathcal{O}(10^{-15})$) for high and low frequencies, respectively. ECRM, in contrast, provides a uniformly good approximation quality of order $\mathcal{O}(10^{-5})$, independently of the chosen frequency. In the $H_{\infty}$- and the $H_{2}$-norms it even outperforms the moment matching versions by one up to two orders.

Remark 4**.**

Alternatively to (4.4), the damped mass-spring system can be also formulated by keeping the original constraint $Gp=0$ and adding the additional constraint $Gv=0$ to match the symmetry structure. This is called index reduction by minimal extension, see [23], and reduces the differentiation-index to two in the system

[TABLE]

The resulting system (4.6) is of size $n=2g+2$, has a regular system matrix and a regular matrix pencil, but is not in port-Hamiltonian form. However, using again the singular value decomposition of $G^{T}$, solving the last equation and inserting the derivative yields a pHODE for $[p_{2}^{T}\,\,v_{2}^{T}]:[0,T]\rightarrow\mathbb{R}^{n_{1}}$, $n_{1}=2(g-1)$ of the form

[TABLE]

whose interconnection matrix is invertible. The other variables satisfy $v_{1}=p_{1}=\lambda_{2}=0$ and $Z^{T}\lambda_{1}=-M_{12}\dot{v}_{2}+K_{12}p_{2}+D_{12}v_{2}+B_{1}u$. To this formulation, also FCRM and MM at $s_{0}=0$ are applicable. While ECRM and MM at $s_{0}=\infty$ yield analogous results independent of the problem formulation, MM at $s_{0}=0$ for (4.7) provides a slightly better $H_{\infty}$-approximation than MM at $s_{0}=10^{-10}$ for (4.5). FCRM shows in general a similar approximation behavior as ECRM but suffers from an error drift off for high frequencies caused by its additional feed-through terms (cf. Figure 4.4). **

5. Conclusion

The power conservation methods (ECRM and FCRM) as well as moment matching via Galerkin projections are established structure-preserving model reduction techniques for standard port-Hamiltonian systems of ordinary differential equations. In this paper we have adapted them to handle also port-Hamiltonian differential-algebraic systems of differentiation-index one or two. Making use of an appropriate decoupling of differential and algebraic variables, the dynamic state is reduced, while the properties and all explicit and hidden constraints of the pHDAE are preserved. The performance of the techniques has been illustrated for benchmark problems arising from spatially discretized flow problems and multibody systems. ECRM shows similarities to Balanced Truncation, if a Lyapunov balancing is performed. Therefore, as expected, ECRM outperforms moment matching when studying the reduction errors in $H_{\infty}$- and/or $H_{2}$-norms, whereas moment matching yields better local approximations in the spectral norm. The performance of FCRM is comparable to ECRM, but it may suffer from an error increase for high frequencies caused by the feed-through terms generated in the reduced model. Moreover, its applicability is limited.

Acknowledgements

The authors acknowledge the support by the German BMBF, Project EiFer, and the German BMWi, Project MathEnergy.