On structure-preserving model reduction for damped wave propagation in transport networks

TL;DR

This paper develops structure-preserving model reduction techniques for damped wave equations in pipeline networks, ensuring key physical properties are maintained during discretization and reduction, with practical algorithms based on Krylov subspaces.

Contribution

It introduces a novel structure-preserving Galerkin projection combined with Krylov subspace methods for model reduction of PDE-DAE systems in pipeline networks.

Findings

Reduced models preserve energy dissipation and passivity.

Algorithms ensure well-posedness and stability of reduced systems.

Numerical tests confirm the effectiveness of the proposed methods.

Abstract

We consider the discretization and subsequent model reduction of a system of partial differential-algebraic equations describing the propagation of pressure waves in a pipeline network. Important properties like conservation of mass, dissipation of energy, passivity, existence of steady states, and exponential stability can be preserved by an appropriate semi-discretization in space via a mixed finite element method and also during the further dimension reduction by structure preserving Galerkin projection which is the main focus of this paper. Krylov subspace methods are employed for the construciton of the reduced models and we discuss modifications needed to satisfy certain algebraic compatibility conditions; these are required to ensure the well-posedness of the reduced models and the preservation of the key properties. Our analysis is based on the underlying infinite dimensional problem and its Galerkin approximations. The proposed algorithms therefore have a direct interpretation in function spaces; in principle, they are even applicable directly to the original system of partial differential-algebraic equations while the intermediate discretization by finite elements is only required for the actual computations. The performance of the proposed methods is illustrated with numerical tests and the necessity for the compatibility conditions is demonstrated by examples.

Figures

Tables

Table 1. Table 6.1. Initial values of m h ​ ( 0 ) subscript 𝑚 ℎ 0 m_{h}(0) and E h ​ ( 0 ) subscript 𝐸 ℎ 0 E_{h}(0) for the mass and energy for full and the reduced order models obtained by projection in the energy norm with and without additional mass constraint. Only the left input u 1 subscript 𝑢 1 u_{1} at vertex v 1 subscript 𝑣 1 v_{1} was used in the Krylov iteration and the shift parameter was set to s 0 = 0 subscript 𝑠 0 0 s_{0}=0 .

		exact	${𝕍}_i={𝕎}_i^L$			${𝕍}_i={𝕎}_i^L+{\mathbb{Z}}_i$
	$L$		$1$	$3$	$10$	$1$	$3$	$10$
projection	$m_h$	$1.000$	$0.750$	$0.902$	$0.949$	$1.000$	$1.000$	$1.000$
	$E_h$	$0.500$	$0.375$	$0.451$	$0.475$	$0.500$	$0.500$	$0.500$
mass constraint	$m_h$	$1.000$	$1.000$	$1.000$	$1.000$	$1.000$	$1.000$	$1.000$
	$E_h$	$0.500$	$0.667$	$0.554$	$0.527$	$0.500$	$0.500$	$0.500$

Table 2. Table 6.2. Energy decay E h ​ ( t ) subscript 𝐸 ℎ 𝑡 E_{h}(t) for the full order model, the standared reduced model 𝕎 i = 𝕍 i subscript 𝕎 𝑖 subscript 𝕍 𝑖 \mathbb{W}_{i}=\mathbb{V}_{i} , and the improved reduced model 𝕎 i = 𝕍 i + ℤ i subscript 𝕎 𝑖 subscript 𝕍 𝑖 subscript ℤ 𝑖 \mathbb{W}_{i}=\mathbb{V}_{i}+\mathbb{Z}_{i} . Only the left input u 1 subscript 𝑢 1 u_{1} and a shift parameter s 0 = 0 subscript 𝑠 0 0 s_{0}=0 is used for the Krylov iteration.

	exact	${𝕍}_i={𝕎}_i^L$			${𝕍}_i={𝕎}_i^L+{\mathbb{Z}}_i$
$t\setminus L$		$1$	$3$	$10$	$1$	$3$	$10$
$0$	$0.5000$	$0.3750$	$0.4512$	$0.4745$	$0.5000$	$0.5000$	$0.5000$
$1$	$0.1528$	$0.3750$	$0.1665$	$0.1514$	$0.1876$	$0.1552$	$0.1527$
$2$	$0.0512$	$0.3750$	$0.0708$	$0.0509$	$0.0690$	$0.0511$	$0.0511$
$3$	$0.0174$	$0.3750$	$0.0384$	$0.0173$	$0.0245$	$0.0173$	$0.0174$
$4$	$0.0059$	$0.3750$	$0.0273$	$0.0059$	$0.0084$	$0.0059$	$0.0059$

Equations

$$a^e{\partial }_t{p}^e+{\partial }_x{q}^e$$

$$b^e{\partial }_t{q}^e+{\partial }_x{p}^e+d^e{q}^e$$

$${\sum }_{e\in \mathcal{E}(v)}{n}^e(v)q^e(v)$$

$$p^e(v)$$

$$y_v$$

$$p(0)=p_0,q(0)=q_0\text{on }\mathcal{E}.$$

$$\frac{d}{dt}m$$

$$=-{\sum }_{e\in \mathcal{E}}{\int }_e{\partial }_x{q}^{e}dx=-{\sum }_{v\in {\mathcal{V}}_{\partial }}{q}^e(v)n^e(v)={\sum }_{v\in {\mathcal{V}}_{\partial }}{y}_v.$$

$$\frac{d}{dt}E$$

$$=-{\sum }_{e\in \mathcal{E}}{\int }_e{d}^e|q^e{|}^{2}dx+{\sum }_{v\in {\mathcal{V}}_{\partial }}{y}_v{u}_v.$$

$$E(t)\le C{e}^{-\gamma (t-s)}E(s),t\ge s,$$

$$\displaystyle M_{1}\dot{x}_{1}\$$

$$\displaystyle M_{2}\dot{x}_{2}\$$

$$\displaystyle\ Nx_{2}\$$

$$y=B_2^{\top }{x}_2.$$

$$E\dot{x}+Ax$$

$$m_h=o_1^{\top }{M}_1{x}_1\text{and}{E}_h=\frac{1}{2}\left(x_1^{\top }{M}_1{x}_1+x_2^{\top }{M}_2{x}_2\right),$$

$$\displaystyle\widehat{M}_{1}\dot{z}_{1}\$$

$$\displaystyle\widehat{M}_{2}\dot{z}_{2}\$$

$$\displaystyle\ \widehat{N}z_{2}\$$

$$\hat{E}\dot{z}+\hat{A}z=\hat{B}u,\hat{y}={\hat{B}}^{\top }\hat{z}.$$

$${\hat{m}}_h={\hat{o}}_1^{\top }{\hat{M}}_1{z}_1\text{and}{\hat{E}}_h=\frac{1}{2}\left(z_1^{\top }{\hat{M}}_1{z}_1+z_2^{\top }{\hat{M}}_2{z}_2\right)$$

$${\mathbb{V}}_1={\mathbb{W}}_1^L+{\mathbb{Z}}_1,{\mathbb{V}}_2={\mathbb{W}}_2^L+{\mathbb{Z}}_2,\text{and}{\mathbb{V}}_3=\mathcal{R}(I_3),$$

$$E\dot{x}+Ax=Bu,y=B^{\top }x,$$

$$H(s)$$

$$(s_{0}E+A)r_0$$

$$(s_{0}E+A)r_l$$

$$\hat{E}\dot{z}+\hat{A}z$$

$$\hat{H}(s)={\hat{B}}^{\top }(s\hat{E}+\hat{A}{)}^{-1}\hat{B}={\sum }_{l=0}^{\infty }{\hat{m}}_l(s_0-s{)}^l,$$

$$E=\left[\begin{array}{ccc}{M}_1& 0& 0\\ 0& M_2& 0\\ 0& 0& 0\end{array}\right],A=\left[\begin{array}{ccc}0& G& 0\\ -G^{\top }& D& -N^{\top }\\ 0& N& 0\end{array}\right],B=\left[\begin{array}{c}0\\ B_2\\ 0\end{array}\right],$$

$$0=x^{\top }(sE+A)x=s(x_1^{\top }{M}_1{x}_1+x_2^{\top }{M}_2{x}_2)+x_2^{\top }D{x}_2.$$

$$0=y_2=s{M}_2{x}_2-G^{\top }{x}_1+D{x}_2-N^{\top }{x}_3=-N^{\top }{x}_3.$$

$$\tilde{A}=\left[\begin{array}{ccc}{A}_{22}& A_{21}& A_{23}\\ A_{12}& A_{11}& A_{13}\\ A_{32}& A_{31}& A_{33}\end{array}\right]=\left[\begin{array}{ccc}D& -G^{\top }& -N^{\top }\\ G& 0& 0\\ N& 0& 0\end{array}\right].$$

$$V=\left[\begin{array}{ccc}{V}_1& 0& 0\\ 0& V_2& 0\\ 0& 0& V_3\end{array}\right]\text{and}{V}_3=I_3.$$

Full text

**Keywords: partial differential-algebraic equations, port-Hamiltonian systems, Galerkin projection, structure-preserving model reduction, passivity, exponential stability **

**AMS-classification (2000): 35L05, 35L50, 65L20, 65L80, 65F25, 65M60 **

1. Introduction

We study a system of partial differential-algebraic equations modeling the propagation of pressure waves in a pipeline network. The basic features of this problem are conservation of mass and dissipation of energy by friction which in turn yields passivity and exponential stability of the system and the convergence to unique steady states. All these properties can be preserved for an appropriate semi-discretization in space by mixed finite elements resulting in a finite dimensional differential-algebraic system with a port-Hamiltonian structure [17]. In this paper, we consider a further dimension reduction of these high dimensional models by structure-preserving Galerkin projection with the aim to obtain reduced models of smaller dimension which can be used for online simulation and control. These models should also yield a good approximation of the overall behavior and preserve the port-Hamiltonian structure and further relevant properties.

The model reduction of structured linear time-invariant systems has attracted significant interest in the literature, see e.g. [4, 20, 33, 35, 36, 41] and the references given there. Related results for second order systems have been obtained in [5, 6, 13, 38, 34], and the reduction of differential-algebraic equations has, for instance, been addressed in [3, 32]. Let us refer to [1, 9] for a general introduction to reduced order modeling and further references.

It is well-known that the port-Hamiltonian structure and thus passivity of the underlying system are inherited automatically by reduced models obtained via structure-preserving Galerkin projection [24, 37, 39]. The preservation of further properties, like conservation of mass or uniform exponential stability, however, requires the bases of the reduced models to satisfy additional compatibility conditions which have to be guaranteed explicitly.

The reduction of infinite dimensional systems described by partial differential or partial differential-algebraic equations has been considered, e.g., in [14, 26, 27], and in [23] the reduction of models arising in gas transport networks has been discussed. For such problems, or discretizations thereof, the bases for the reduced models have to be generated by some iterative process. Krylov subspace methods [2, 19, 22, 37] and proper orthogonal decomposition [14, 26, 27] are frequently employed for this purpose, and their analysis in a function space setting allows to obtain mesh independent results.

In this paper we consider a structure-preserving model reduction for large scale differential-algebraic systems obtained by discretization of a partial differential-algebraic model. We utilize Krylov subspace methods for the basis construction together with a structure-preserving space splitting and discuss appropriate modifications in order to satisfy some compatibility conditions required for the proof of mass conservation, uniform exponential stability, and the existence of steady states. While our algorithms are formulated in an algebraic setting, they also have an interpretation in function spaces. This is used already for the formulation of our algorithms and allows a complete analysis of the reduced models. Our methods therefore turn out to be almost independent of the intermediate finite element approximation used in computations and they are applicable, in principle, even directly to the underlying partial differential-algebraic system.

The outline of the paper is as follows: In the following section, we introduce the model problem under consideration and discuss the basic steps and arguments of our approach. The remainder of the manuscript is then split into three major parts: Part I is concerned with an outline and a partial analysis of the model reduction approach and Part II provides numerical illustration of these results. Part III contains the full analysis of the reduced order models obtained with our approach which requires us to consider the infinite dimensional problem and its approximation by mixed finite elements. The corresponding results mostly follow from those in [17] and they are therefore presented in the appendix for completeness and convenience of the reader.

2. Model problem and outline of the approach

The purpose of this section is to introduce in detail the problem under consideration and to give a rough idea of our approach and of the mutual relations between the underlying infinite dimensional system, the large scale finite dimensional systems arising after discretization in space, and the reduced models we are looking for.

2.1. Model problem

We consider the propagation of pressure waves in a one-dimensional network of pipes whose geometry shall be given as finite directed and connected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertices $v\in\mathcal{V}$ and edges $e\in\mathcal{E}$. On every pipe $e$, the conservation of mass and the balance of momentum are described by

[TABLE]

Here $p^{e}$, $q^{e}$ denote the pressure and mass flux which are functions of space and time, the coefficients $a^{e}$, $b^{e}$ encode properties of the fluid and the pipe, and $d^{e}$ models the damping due to friction at the pipe walls. The coefficients are assumed to be positive and, for ease of presentation, constant on every pipe $e$. At every inner vertex $v\in{\mathcal{V}_{0}}$ of the graph, corresponding to a junction of several pipes $e\in\mathcal{E}(v)$, we require that

[TABLE]

Here $n^{e}(v)=\mp 1$, depending on whether the pipe $e$ starts or ends at the vertex $v$; see Figure 2.1. Furthermore, $m^{e}(v)$, $p^{e}(v)$ denote the respective functions evaluated at the vertex $v$ but still depending on time.

These coupling conditions model the conservation of mass and momentum at the junctions. At the boundary vertices $v\in{\mathcal{V}_{\partial}}=\mathcal{V}\setminus{\mathcal{V}_{0}}$, which correspond to the the ports of the network, we set

[TABLE]

with values $u_{v}$ denoting the given input at the port $v\in{\mathcal{V}_{\partial}}$. As corresponding output of the system, we consider the mass flux via the ports, given by

[TABLE]

Other input and output configurations could be considered without difficulty as well. The specification of the model is completed by assuming knowledge of the initial conditions

[TABLE]

The system (2.1)–(2.7) models the propagation of pressure waves in a gas network on the acoustic time scale [12]. For sufficiently smooth initial data $p_{0}$, $q_{0}$, and appropriate compatible input functions $(u_{v})_{v\in{\mathcal{V}_{\partial}}}$, existence of a unique classical solution can be established [17].

2.2. Basic properties

The partial differential-algebraic system (2.1)–(2.5) encodes several interesting properties which are directly related to the underlying physical principles:

(P1) global conservation of mass, which here can be expressed as

[TABLE]

The total mass of the gas contained in the system can thus only be altered by flow of gas into or out of the system via the ports of the network.

(P2) a port-Hamiltonian structure, leading to energy dissipation and passivity, i.e.,

[TABLE]

The total energy in the system only changes by dissipation through damping and injection or extraction via the system ports. Apart from these basic properties, the system further admits

(P3) exponential stability and convergence to equilibrium for input $u\equiv 0$; more precisely,

[TABLE]

with constants $C$ and $\gamma$ that are independent of the particular solution. For constant input $u(t)\equiv const$, one can thus observe exponential convergence to

(P4) unique steady states for the corresponding stationary problem.

For a proof of these properties, let us refer to the appendix and to [17]. Systems of similar structure and with similar properties also model networks of electric transmission lines, vibrations of elastic multi-structures, or more general wave phenomena on multiply connected domains. Our arguments therefore may be useful in a wider context; see [9, 23, 30, 40] for further applications.

2.3. Full order model

An appropriate discretization of the partial differential-algebraic system in space by mixed finite elements leads to a differential-algebraic system

[TABLE]

which we will call the full order model in the sequel. In the context of reduced basis methods, the notion truth approximation is sometimes used instead. The vectors $x_{1}$ and $x_{2}$ are the algebraic representations of the states $p$ and $q$ after discretization, and $x_{3}$ resembles the Lagrange multiplier for the constraint (2.3). The output of the system is then given by

[TABLE]

If an appropriate discretization is used, the system matrices can be shown to have some basic structural properties. For ease of presentation, we formulate them here as assumptions:

• (A0)

$M_{1}$, $M_{2}$, $D$ are symmetric and positive definite and $[G^{\top},N^{\top}]$ has trivial null-space.

The latter condition is equivalent to requiring that $G$ and the restriction of $N$ to the nullspace of $G$ define surjective linear operators.

Remark 2.1** (Notation).**

Throughout the paper, we identify matrices with corresponding linear operators. We call a matrix $A$ injective or surjective, if the operator has the respective property, and we write $\mathcal{R}(A)$ and $\mathcal{N}(A)$ for the range and the kernel of the corresponding operator. Furthermore $A\mathbb{V}$ denotes the image of the space $\mathbb{V}$ under the map induced by $A$.

Remark 2.2**.**

The differential-algebraic system (2.8)–(2.10) can be shown to formally have differentiation-index two [10, 29]. The condition that $N^{\top}$ is injective, and hence that $N$ is surjective, however, allows to eliminate the Lagrange multiplier by purely algebraic manipulations and hence to reduce the system to an ordinary differential equation; see Section 6.2 and also refer to [8, 18, 25] for more general situations. Let us emphasize that the number of constraints amounts to the number of junctions in the network and thus is finite here.

Remark 2.3**.**

The system (2.8)–(2.10) can be written as linear time-invariant descriptor system

[TABLE]

From the particular form of the matrices $E$ and $A$ one can directly deduce the port-Hamiltonian structure, i.e., $E$ is symmetric and positive semi-definite and $A=J+R$ can be decomposed into an skew-symmetric part $J$ and a symmetric positive semi-definite part $R$. This immediately guarantees the passivity of the system and further useful properties [42, 43].

The mass and the energy of the semi-discrete system (2.8)–(2.10) can be expressed as

[TABLE]

where $o_{1}$ is the vector representing the constant one function on the network. The basic properties (P1)–(P4) can then be shown to hold almost verbatim also for the semi-discrete problem which may therefore serve as a replacement for the infinite dimensional partial differential-algebraic problem under investigation.

2.4. Structure preserving model reduction

The main focus of the current paper is a further dimension reduction of the differential-algebraic model (2.8)–(2.10) by structure-preserving Galerkin projection of the following form: Given projection matrices $V_{1},V_{2}$ of appropriate size and full rank, we set $\widehat{M}_{i}=V_{i}^{\top}M_{i}V_{i}$, $\widehat{D}=V_{2}^{\top}DV_{2}$, $\widehat{B}_{2}=V_{2}^{\top}B_{2}$, $\widehat{G}=V_{2}^{\top}GV_{1}$, and $\widehat{N}=NV_{2}$. The reduced model is then defined as

[TABLE]

The tuple $(V_{1}z_{1},V_{2}z_{2},z_{3})$ is the approximation for the exact solution $(x_{1},x_{2},x_{3})$ of the full order model and $\widehat{y}=\widehat{B}_{2}^{\top}z_{2}$ serves as approximation for the output $y=B_{2}^{\top}x_{2}$ of the full system.

Remark 2.4**.**

Note that the dimension of the space for the Lagrange multiplier $x_{3}$ has not been reduced in the above construction; the network topology is thus completely maintained. The reduced model can again be written in the form of a descriptor system

[TABLE]

It is well-known [24, 37, 39] and easy to see for the system considered here that the port-Hamiltonian structure and thus passivity are inherited automatically by this kind of Galerkin projection. Additional conditions will, however, be required to establish the well-posedness of the resulting reduced differential-algebraic system and to characterize its index; see Section 3.3.

Similar as before, we will denote by

[TABLE]

the mass and energy of the reduced problem (2.14)–(2.16). An appropriate vector $\widehat{o}_{1}$ representing the constant one function on the network will be needed and additional compatibility conditions will be required to ensure well-posedness of the reduced system and the validity of (P1)–(P4).

2.5. Algebraic compatibility conditions

As we will demonstrate by explicit examples below, the validity of some of the properties (P1)–(P4) and even the well-posedness of the reduced models can in general not be guaranteed, unless additional assumptions on the projection matrices $V_{i}$ are satisfied. We will therefore require that

• (A1)

$o_{1}\in\mathcal{R}(V_{1})$;

• (A2)

$\mathcal{R}(M_{1}V_{1})=\mathcal{R}(GV_{2})$;

• (A3)

$\mathcal{N}(G)\subset\mathcal{R}(V_{2})$ and $N\mathcal{N}(G)=\mathcal{R}(I_{3})$.

Recall that $\mathcal{R}(A)$ and $\mathcal{N}(A)$ denote the range and the nullspace of the linear operator induced by a matrix $A$ and $A\mathbb{V}$ is the the image of the space $\mathbb{V}$ under mapping induced by $A$. Further, $I_{3}$ here denotes the identity matrix for the third component and $o_{1}$ is the vector used to describe the total mass $m_{h}=o_{1}^{\top}M_{1}x_{1}$ of the full order system (2.8)–(2.10).

Remark 2.5**.**

Assumption (A1) will allow us to prove the conservation of mass also for the reduced models. The conditions (A2)–(A3), on the other hand, allow us to show that

• ($\widehat{\text{A0}}$)

$\widehat{M}_{1}$, $\widehat{M}_{2}$, and $\widehat{D}$ are symmetric and positive definite and $[\widehat{G}^{\top},\widehat{N}^{\top}]$ has trivial nullspace;

see Lemma 3.4 for details. The reduced system thus has the same algebraic properties as the full order model (2.8)–(2.10). The well-posedness of the reduced system (2.14)–(2.16) can therefore be obtained with similar arguments as that of the full order model. The differentiation-index of the reduced model is again two and, by elimination of the Lagrange multiplier, we can obtain a regular system of ordinary differential equations; see Section 6.2 for details.

2.6. Basis construction

For the actual construction of the projection matrices $V_{i}$, we consider an extension of the approach proposed in [20] together with some modifications in order to satisfy the compatibility conditions (A1)–(A3). The main steps can be sketched as follows:

• •

Krylov iteration: construct finite dimensional subspaces $\mathbb{W}^{L}$ with good approximation properties by a Krylov iteration applied to the full order model (2.8)–(2.10).

• •

Splitting: Decompose $\mathbb{W}^{L}$ as $\mathbb{W}^{L}=(\mathbb{W}_{1}^{L},0,0)+(0,\mathbb{W}^{L}_{2},0)+(0,0,\mathbb{W}^{L}_{3})$ according to the components of the state $x=(x_{1},x_{2},x_{3})$.

• •

Modification: choose appropriate subspaces $\mathbb{Z}_{1}$ and $\mathbb{Z}_{2}$ and define

[TABLE]

such that the properties (A1)–(A3) can be verified for any choice of $V_{i}$, $i=1,2,3$, whose columns form bases for the corresponding subspaces.

With similar arguments as in [20], the reduced models (2.14)–(2.16) can be shown to match certain moments of the transfer function and thus to have good approximation properties. By construction, the projection matrices also satisfy the compatibility conditions (A1)–(A3). This will allow us to show that the reduced models are well-posed and that they inherit the structural properties (P1)–(P4) from the full order model.

2.7. Overview

The derivation of the properties (A0) for the system matrices of the full order model and of the algebraic compatibility conditions (A1)–(A3), as well as the complete analysis of the resulting reduced order models require us to consider in detail the connection between

• •

the underlying partial differential-algebraic equations;

• •

their discretization by Galerkin approximations in a function space setting; and

• •

the corresponding linear time-invariant systems in algebraic form.

A sketch of these different viewpoints is depicted in Figure 2.2. The close relation of the differential-algebraic systems to the problem on the continuous level will allow us to establish properties of the reduced order models that are uniform and almost independent of the intermediate full order model which is only required for the actual computations.

Apart from these analytical considerations, we also investigate in detail the algorithms for the actual subspace construction on the algebraic level and we address the following issues:

• •

The splitting step in the subspace construction turns out to be sensitive to numerical errors. To overcome this, we utilize a cosine-sine decomposition in the final algorithm.

• •

Round-off errors affect the validity of (A2) after the modification step outlined above. We therefore take special care in the basis construction to satisfy (A2) explicitly.

In order to reflect the functional analytic setting of the underlying infinite dimensional problem, we will utilize appropriate scalar products in the formulation of our algorithms on the algebraic level; see [27, 28] for similar approaches. As a consequence, the vectors obtained in the basis construction process can be interpreted as functions on the continuous level which allows for a further evaluation and interpretation of the numerical results.

Part I: Model reduction

In the following two sections, we present our model reduction approach on the algebraic level. We discuss in detail the construction of the reduced models and investigate their approximation properties. Futhermore, we address some algorithmic details.

3. Structure preserving model reduction

Let us first recall some basic facts about model order reduction and then informally discuss the algebraic compatibility conditions which are at the core of our model reduction approach. The basis construction algorithms on the algebraic level will then be presented in the next section.

3.1. Model reduction basics

Consider a general linear time-invariant descriptor system

[TABLE]

where $E$, $A$, $B$ are given matrices, $E$ symmetric and positive semi-definite, and $sE+A$ defines a regular matrix pencil, i.e., $sE+A$ is regular for almost every $s\in\mathbb{R}$; see [10, 15, 29, 31] for details. A formal expansion of the transfer function of the system [1, 9] leads to

[TABLE]

here $s_{0}\in\mathbb{C}$ is some given shift parameter. It is not difficult to see that the generalized moments $m_{l}$ can be written as $m_{l}=B^{\top}r_{l}$ with vectors $r_{l}$ that can be computed recursively by

[TABLE]

Let us denote by $\mathbb{W}^{L}=\text{span}\{r_{0},\ldots,r_{L-1}\}$ the $L$th Krylov subspace generated by this iteration.

For any given projection matrix $V$ of appropriate dimension and maximal rank, we define $\widehat{E}=V^{\top}EV$, $\widehat{A}=V^{\top}AV$, and $\widehat{B}=V^{\top}B$, and consider the reduced system

[TABLE]

resulting from Galerkin projection of (3.1) onto the range of the matrix $V$. The transfer function of this reduced model may again be expanded as

[TABLE]

and the relation between the full and the reduced order model can be characterized as follows.

Lemma 3.1** (Moment matching).**

Let $\mathbb{W}^{L}\subset\mathcal{R}(V)$. Then $\widehat{m}_{l}=m_{l}$ for $l=0,\ldots,2L-1$.

A proof of this assertion and further results can be found in [1, 9, 22]. The lemma can be interpreted as an abstract approximation result and may therefore serve as a quality indicator.

3.2. Basic properties of the full order model

The system (2.8)–(2.10) can be written in the compact form (3.1) with state vector $x=(x_{1},x_{2},x_{3})$ and with system matrices defined by

[TABLE]

and under basic structural assumptions, the system (3.1) can be shown to be well-posed.

Lemma 3.2**.**

Let (A0) hold and let $E$ and $A$ be defined as above. Then $sE+A$ defines a regular matrix pencil. Moreover, $sE+A$ is regular for any $s\geq 0$ and, in particular, $A$ is regular.

Proof.

First consider the case $s>0$ and set $x=(x_{1},x_{2},x_{3})$. We show that $0=y=(sE+A)x$ implies $x=0$. Multiplication with $x^{\top}$ from the left yields

[TABLE]

Since $M_{i}$ and $D$ are positive definite, this implies that $x_{1}=0$ and $x_{2}=0$. But then

[TABLE]

From the assumptions on $G$ and $N$, we can deduce the injectivity of $N^{\top}$, and hence $x_{3}=0$.

Now consider the case $s=0$: By simple rearrangement of the blocks of the matrix $A$, we obtain

[TABLE]

Such a system is regular, if, and only if, $D$ is regular and $[G^{\top},N^{\top}]$ has trivial nullspace; see [11] for the corresponding result in infinite dimensions. These properties are guaranteed by assumption (A0), which yields the invertibility of $\tilde{A}$ and hence also of $A$. ∎

Remark 3.3**.**

Only the injectivity of $N^{\top}$, or equivalently, the surjectivity of $N$ is required to obtain a regular matrix pencil $sE+A$ for $s>0$ and thus to establish well-posedness for the time-dependent problem. The stronger condition that $[G^{\top},N^{\top}]$ is injective is, however, necessary to obtain regularity of the matrix $A$ and thus to ensure existence of unique steady states.

3.3. Properties of the reduced problem

The reduced descriptor system (3.5), representing the reduced order model (2.14)–(2.16), can be obtained by Galerkin projection of the descriptor system (3.1), which represents (2.8)–(2.10), with a projection matrix $V$ of the form

[TABLE]

Recall that $I_{3}$ is the identity matrix for the space of Lagrange multipliers and note again that we did not reduce the number of constraints here. For a projection matrix $V$ of this form, the particular algebraic structure of the full order model is directly passed on to the reduced order model. The compatibility conditions (A2)–(A3) further allow us to establish the structural properties corresponding to (A0) also for the system matrices of the reduced model.

Lemma 3.4**.**

Let (A0) and (A2)–(A3) hold and let $V_{1}$, $V_{2}$ be injective. Then ($\widehat{\text{A0}}$) holds.

Proof.

The conditions on $\widehat{M}_{i}$ and $\widehat{D}$ are clearly satisfied, if (A0) is valid and the columns of the projection matrices $V_{i}$ are linearly independent. Due to assumption (A2), we can find for any vector $z_{1}$ a corresponding vector $z_{2}$ such that $GV_{2}z_{2}=M_{1}V_{1}z_{1}$. With this choice, we obtain

[TABLE]

Since $\widehat{M}_{i}$ is symmetric positive definite, we obtain $z_{1}^{\top}\widehat{M}_{1}z_{1}>0$ whenever $z_{1}\neq 0$, and consequently $\widehat{G}^{\top}$ is injective, or equivalently, $\widehat{G}$ is surjective. Using assumptions (A0) and (A3), we can further find for any $x_{3}=z_{3}$ a vector $x_{2}\in N(G)$ and a vector $z_{2}$ such that

[TABLE]

This shows that the restriction of $\widehat{N}$ to the nullspace of $\widehat{G}$ is surjective. Together with $\widehat{G}$ being surjective this is equivalent to $[\widehat{G}^{\top},\widehat{N}^{\top}]$ being injective. ∎

With the same reasoning as in Lemma 3.2, we now obtain

Lemma 3.5**.**

Let (A0) and (A2)–(A3) hold and $V_{1}$, $V_{2}$ be injective. Then the matrix pencil of the reduced problem is regular. In particular, $s\widehat{E}+\widehat{A}$ is regular for $s\geq 0$ and thus $\widehat{A}$ is regular.

Remark 3.6**.**

As a consequence of Lemma 3.5, we see that the system (2.14)–(2.16) representing the reduced problem is well-posed and possesses unique steady states. By construction, the system also inherits the port-Hamiltonian structure and passivity. Furthermore, the reduced differential-algebraic system again has differentiation-index two, and due to injectivity of $N^{\top}$, the constraints can be eliminated algebraically and the reduced system can thus again be reduced to an ordinary differential equation.

4. Subspace and basis construction

It is well understood [1, 9] that the model reduction approach outlined above amounts to a Galerkin projection of the full order system onto subspaces $\mathbb{V}_{i}=\mathcal{R}(V_{i})$ generated by the columns of the projection matrices $V_{i}$, $i=1,2,3$. In the following, we will change between the algebraic viewpoint and that of function spaces viewpoint as convenient.

We consider a construction of the subspaces $\mathbb{V}_{i}$ in the form

[TABLE]

A Krylov iteration [20] together with an appropriate splitting is used for generation of the spaces $\mathbb{W}_{1}$, $\mathbb{W}_{2}$, which by Lemma 3.1 automatically ensures good approximation properties of the resulting reduced model. The spaces $\mathbb{Z}_{1}$, $\mathbb{Z}_{2}$, on the other hand, will be chosen in order to guarantee the compatibility conditions (A1)–(A3). After definition of the subspaces $\mathbb{V}_{1}$ and $\mathbb{V}_{2}$, we also present the algorithms for the actual computation of the projection matrices $V_{1}$ and $V_{2}$.

4.1. Construction of the spaces $\mathbb{W}_{1}$ and $\mathbb{W}_{2}$

We start with applying the Krylov subspace iteration (3.3)–(3.4) to the particular system (2.8)–(2.10). The first step now reads

[TABLE]

and for $l\geq 1$ the further iterations are defined accordingly by

[TABLE]

As a direct consequence of Lemma 3.2, we obtain

Lemma 4.1** (Krylov subspace iteration).**

Let the assumption (A0) be valid. Then the iterative scheme (4.2)–(4.7) is well-defined for any $s_{0}\geq 0$ and all $l\geq 0$.

Let us denote the subspaces spanned by the iterates $x_{i}^{l}$ generated by (4.2)–(4.7) with

[TABLE]

The special structure of the iteration allows us to show the following result.

Lemma 4.2** (Properties of subspaces).**

Let $s_{0}>0$, then $M_{1}\mathbb{W}_{1}^{L}=G\mathbb{W}_{2}^{L}$, while for $s=0$ we have $M_{1}\mathbb{W}_{1}^{L}=G\mathbb{W}_{2}^{L+1}$. Moreover, $Nx_{2}=0$ for any $x_{2}\in\mathbb{W}_{2}^{L}$ and any choice of $s_{0}\geq 0$.

Proof.

First consider $s_{0}>0$: Then from (4.2), we see that $M_{1}x_{1}^{0}=-Gx_{2}^{0}$, and thus the first assertion holds true for $L=0$. By (4.5) and induction, we obtain $M_{1}\mathbb{W}_{1}^{L}=G\mathbb{W}_{2}^{L}$ for all $L\geq 1$. Now consider the case $s_{0}=0$: Then from (4.2), one can deduce that $x_{2}^{0}\in\mathcal{N}(G)$. Using (4.5) with $L=1$, we further see that $M_{1}x_{1}^{0}=Gx_{2}^{1}$, and hence $M_{1}\mathbb{W}_{1}^{1}=G\mathbb{W}_{2}^{2}$. By (4.5) and induction, we finally obtain the result for all $L\geq 1$ again. ∎

4.2. Construction of the spaces $\mathbb{Z}_{1}$ and $\mathbb{Z}_{2}$

The structural assumption (A0) particularly implies that $G$ is surjective which allows us to ensure the following property.

Lemma 4.3**.**

Let assumption (A0) hold. Then for any choice of $g_{1}$ and $g_{3}$, the system

[TABLE]

has at least one solution $x_{2}$. We denote by $x_{2}^{\dagger}$ the unique solution that also minimizes $x_{2}^{\top}M_{2}x_{2}$. This minimum-norm solution $x_{2}^{\dagger}$ can be expressed as

[TABLE]

denoting the pseudo-inverse with respect to the scalar product induced by $M_{2}$.

We can now define $o_{2}$ as the minimum-norm solution of $Go_{2}=M_{1}o_{1}$ and $No_{2}=0$ and set

[TABLE]

If $s_{0}=0$, then we choose $\mathbb{Z}_{2}=\text{span}\{o_{2}\}+\mathcal{N}(G)+\text{span}\{x_{2}^{L}\}$ instead. As a consequence of this construction and the observations made in Lemma 3.1 and 4.2, we obtain the following result.

Lemma 4.4** (Compatibility).**

Define $\mathbb{V}_{1}=\mathbb{W}_{1}^{L}+\mathbb{Z}_{1}$ and $\mathbb{V}_{2}=\mathbb{W}_{2}^{L}+\mathbb{Z}_{2}$, and let $V_{i}$ denote matrices whose columns form bases for $\mathbb{V}_{i}$, $i=1,2$, orthogonal with respect to the scalar products induced by the matrices $M_{i}$, respectively. Then the assumptions (A1)–(A3) are satisfied and the reduced model (2.14)–(2.16) satisfies the moment matching conditions of Lemma 3.1.

The assertions follow directly from the construction. Let us finally mention an equivalent representation for the space $\mathbb{V}_{2}$ which will be used in our basis construction algorithm in Section 5.

Lemma 4.5**.**

Let $\mathbb{V}_{i}$ and $V_{i}$, $i=1,2$ be defined as in Lemma 4.4. Then

[TABLE]

Proof.

The assertion follows by construction of the space $\mathbb{V}_{2}$ and Lemmas 4.2 and 4.3. ∎

5. Algorithms for the basis construction

Following the above considerations, we can now formulate algorithms for the explicit construction of the projection matrices $V_{i}$. We will use Matlab notation throughout this section. Let us start with the construction of the Krylov subspaces, for which we use an Arnoldi method.

Algorithm 5.1** (Krylov iteration).**

% function W=krylov(E,A,B,s0,L,tol) r = (s0E + A)\B; r = ortho(r,[],E,tol); W = r; for l=1:L-1 r = (s0E + A)(E*r); r = ortho(r,W,E,tol); W = [W,r]; end

The orthogonalization can be realized efficiently via a modified Gram-Schmidt process with one re-orthogonalization step [21]. Note that orthogonality is understood with respect to the bilinear form induced by the matrix $E$ here, which also appears in the definition of the energy of the system. The corresponding algorithm reads

Algorithm 5.2** (Orthogonalization).**

% function V=ortho(V,W,E,tol) for k = 1:size(V,2) % orthogonalize to Wj for r=1:2 % use reorthonormalization for j = 1:size(W,2) hk1j = W(:, j)’ * E * V(:, k); V(:, k) = V(:, k) - W(:, j) * hk1j; end end % orthogonalize to previous Vj for r=1:2 for j = 1:k-1 if d(j)<tol, continue; end hk1j = V(:, j)’ * E * V(:, k); V(:, k) = V(:, k) - V(:, j) * hk1j; end end % normalize d(k) = sqrt(V(:,k)’ * E * V(:,k)); if d(k)>=tol, V(:, k) = V(:, k) / d(k); end end % only keep relevant vectors V = V(:,find(d>tol));

The next step consists in the splitting of the matrix $W=[W_{1};W_{2};W_{3}]$ corresponding to the solution components $x=[x_{1};x_{2};x_{3}]$. Note that even if the columns of $W$ are orthogonal, this will in general no longer be true for the columns of $W_{i}$ and thus some re-orthogonalization is required. For reasons of numerical stability, we here employ the cosine-sine decomposition

[TABLE]

where $U_{1}$, $U_{2}$, and $X$ are orthogonal and $C$, $S$ are diagonal with entries $C_{ii}^{2}+S_{ii}^{2}=1$; this explains the name of the decomposition. Note that the cosine-sine decomposition and the related generalized singular value decomposition can be computed efficiently and stably [21, 44]. The following algorithm additionally takes into account non-standard scalar products.

Algorithm 5.3** (Stable splitting via cosine-sine decomposition).**

% function [W1,W2]=split(W1,W2,M1,M2,tol) % compute cholesky factorizations Mi=Ri*Ri’ R1 = chol(M1); R2 = chol(M2);

  % compute generalized svd
  [U1,U2,X,C,S] = gsvd(R1*W1,R2*W2);

  % eliminate dependent columns
  kc = find(diag(C)>tol); ks = find(diag(S)>tol);
  W1 = R1\U1(:,kc);       W2 = R2\U2(:,ks);

For the problems under investigation, the splitting via the cosine-sine decomposition does not cause a substantial computational overhead but significantly improves the stability compared to the simple splitting $W=[W_{1};W_{2}]$ with subsequent re-orthogonalization; see Section 7.2 for an illustration by numerical tests.

As a final step in the basis construction process, we now apply the modifications to ensure the algebraic compatibility conditions (A1)–(A3) which finally allow us to guarantee the properties (P1)–(P4) also for the reduced models. For this purpose, we use the following implementation.

Algorithm 5.4** (Modifications).**

% function [V1,V2]=modify(W1,W2,M1,M2,o1,nullG,tol) V1 = ortho([W1,o1],[],M1,tol); V2 = M2([G’,N’](([G;N](M2[G’,N’]))[M1*W1;zeros(size(N,1),size(V1,2))])); V2 = ortho([nullG,V2],[],M2,tol);

Note that the matrix $V_{2}$ was defined here following the considerations of Lemma 4.5. This again does not significantly increase the overall complexity but substantially improves validity of the compatibility condition (A2) in the presence of round-off errors.

Summary

The previous considerations allow us to draw the following important conclusions which describe the basic properties of our model reduction approach.

Theorem 5.5**.**

Let (A0)–(A3) be valid and let the projection matrices $V_{i}$ be defined with the Algorithms 5.1–5.4. Then the reduced order model (2.14)–(2.16) is well-posed, conserves mass, dissipates energy, and has exponentially stable steady states, i.e., it satisfies (P1)–(P4).

Proof.

Well-posedness and the existence of unique steady states and thus property (P4) follow from Lemma 3.5. Validity of (P2) is a consequence of the structure-preserving Galerkin projection. The proof of (P1) follows by construction and (P3) can be deduced from (P2) and (P4). The proof of uniform exponential stability (P3) and its independence of the intermediate full order model, however, requires a detailed analysis of the underlying partial differential-algebraic model and its Galerkin approximations which will be presented in the appendix. ∎

Part II: Numerical illustration

In the following two sections, we first demonstrate the importance of the algebraic compatibility conditions and then illustrate our main results with some numerical tests.

6. Comparison of reduced models

In all our experiments, the full order model (2.8)–(2.10) is obtained by the mixed finite element discretization of the system (2.1)–(2.6), as discussed in Section A2 and, therefore, the conditions (A0) are valid. Implicit Runge-Kutta methods are used for integration in time, and the time step $\tau$ is chosen sufficiently small in order to minimize the errors due to the time discretization. Quantities $x^{k}$ will denote the approximations for $x(t^{k})$ at time $t^{k}=k\tau$.

6.1. Test problems

We start with some considerations for the most simple networks consisting of a single pipe of unit length and the same pipe split into two parts. The model parameters are set to $a^{e}=b^{e}=d^{e}=1$.

These test cases are already sufficient to illustrate the necessity of the compatibility conditions (A1)-(A3) for the well-posedness of the reduced models and the validity of properties (P1)–(P4). We will compare the reduced models based on Krylov subspaces with and without modifications and refer to these as

• •

the standard (reduced) model: $\mathbb{V}_{1}=\mathbb{W}_{1}^{L}$, $\mathbb{V}_{2}=\mathbb{W}_{2}^{L}$; and

• •

the improved (reduced) model: $\mathbb{V}_{1}=\mathbb{W}_{1}^{L}+\mathbb{Z}_{1}$, $\mathbb{V}_{2}=\mathbb{W}_{2}^{L}+\mathbb{Z}_{2}$;

respectively. If not stated otherwise, we always set the shift parameter to $s_{0}=0$. Moreover, we only utilize one single input for the subspace construction at the boundary vertex $v_{1}$ on the left side of the pipe. In this case $B_{2}$ consists of a single column.

6.2. Reduction to an ordinary differential equation

For the test case (TP1) of a single pipe, the coupling matrix $N$ and the Lagrange parameter $x_{3}$ have zero dimension. Therefore, the model (2.8)–(2.10) is just a system of ordinary differential equations and the reduced model (2.14)–(2.15) is well-posed for any choice of projection matrices $V_{1}$ and $V_{2}$ having full rank.

This is, however, no longer true for the case (TP2) of two pipes, where $N$ and $x_{3}$ have one row. Well-posedness of the full-order system (2.8)–(2.10) and of the improved reduced model based on subspaces $\mathbb{V}_{i}=\mathbb{W}_{i}+\mathbb{Z}_{i}$ is still ensured by Lemmas 4.1 and 4.4. For the standard reduced model based on subspaces $\mathbb{V}_{i}=\mathbb{W}_{i}$, on the other hand, we deduce from Lemma 4.2 that $\widehat{N}=NV_{2}=0$. Hence (A3) does not hold if $\dim\mathbb{V}_{3}\neq 0$, and the matrix pencil $s\widehat{E}+\widehat{A}$ for the reduced model is singular. Consequently, the linear system (2.14)–(2.16) is not well-posed. The condition $NV_{2}=0$ can, however, be used to eliminate the Lagrange multiplier $z_{3}$ and to obtain the smaller system of ordinary differential equations

[TABLE]

Whenever the standard reduced model is used in the following, we will tacitly eliminate the Lagrange multiplier and the constraints in this way. Since the matrices $\widehat{M}_{i}$ are regular, this problem is clearly well-posed. Note, however, that the Lagrange multiplier $z_{3}$ cannot be recovered uniquely unless assumption (A3) is satisfied. Also, a potential input in the third equation (2.16), which might occur in more general situations, cannot be handled appropriately unless this compatibility condition is valid.

6.3. Choice of initial conditions

Let us briefly discuss the choice of initial conditions for the reduced models. By projection with respect to the energy scalar products, we obtain

[TABLE]

This choice provides the best approximation of the initial conditions with respect to the energy of the problem, and the discrete energy $E_{h}=\frac{1}{2}\left(x_{1}^{\top}M_{1}x_{1}+x_{2}^{\top}M_{2}x_{2}\right)$ at initial time is approximated as good as possible by the energy of the reduced model, which is obtained by replacing $x_{i}=V_{i}z_{i}$. Moreover, the initial energy is not increased by the projection step. The total mass of the full order system is defined as $m_{h}=o_{1}^{\top}M_{1}x_{1}$, and that of the reduced models is obtained by replacing $x_{1}=V_{1}z_{1}$ again.

With the above choice of initial conditions, the improved reduced model will exactly reproduce the initial mass. For the standard reduced model, we can in general not assume that $o_{1}\in\mathbb{V}_{1}$, which may lead to a rather large defect in the initial mass. As a remedy, one may enforce the correct representation of the initial mass by a constraint which, however, leads to a potential increase in the initial energy. In Table 6.1, we display the values for the total mass obtained for initial values $p_{0}=1$ and $q_{0}=0$ with the reduced models using these two strategies.

For the standard reduced model based on spaces $\mathbb{V}_{i}=\mathbb{W}_{i}^{L}$, we observe a substantial miss-specification of the total mass at initial time if the initial conditions are chosen by the energy projection. Exact representation of the total mass via the constraint, on the other hand, leads to an artificial increase of the initial energy. The size of both defects can be reduced by increasing the approximation order $L$ which allows to approximate the initial conditions better and better. The improved reduced model, on the other hand, satisfies (A1) by construction and therefore leads to the exact representation of the mass and a good approximation of the energy at the same time. For the problem under investigation, the energy can even be represented exactly.

6.4. Conservation of energy

The port-Hamiltonian structure of the reduced order system (2.8)–(2.10) automatically leads to exact conservation of the sum of total and dissipated energy

[TABLE]

Here $E_{h}=\frac{1}{2}(x_{1}^{\top}M_{1}x_{1}+x_{2}^{\top}M_{2}x_{2})$ is the total energy and $D_{h}=x_{2}^{\top}Dx_{2}$ the dissipation term. This holds true for the standard and the improved reduced model. An additional numerical dissipation term arises if a stable implicit Runge-Kutta scheme is used for integration in time. If the time step size $\tau$ is chosen sufficiently small, the effect of this artificial numerical dissipation can however be mad arbitrarily small and thus be neglected; more details will be given below.

6.5. Exponential stability

We next consider the influence of the basis construction on the exponential stability of the reduced models. We again set $p_{0}=1$ and $q_{0}=0$ and choose the inputs in the boundary conditions as $u_{1}=0$ and $u_{0}=1$ for $t>0$. In Table 6.2, we display the values of the energy for the full order and the reduced models for a sequence of time steps.

The improved reduced model yields uniform exponential decay of the energy in all cases. Already for $L=3$, the energy is predicted accurately over the whole time interval. The standard reduced model, on the other hand, underestimates the initial energy and does not provide the correct decay rate for small $L$. For the smallest model with $L=1$, we do not even observe any decay in energy at all. As we will see next, this defect may occur for any number of moments.

6.6. Existence of steady states

Consider the standard reduced model: As shown in Section 3.3, the unique solvability of the stationary problem requires $[\widehat{G}^{\top},\widehat{N}^{\top}]$ to be injective. For the case (TP1) of a single pipe, the matrix $N$ has dimension zero since no inner vertex exists. Even in this case, according to Lemma 4.2, surjectivity of $\widehat{G}$ and thus injectivity of $\widehat{G}^{\top}$ can in general only be guaranteed for shift parameter $s_{0}>0$. For $s_{0}=0$, we expect to observe a rank deficiency in the standard reduced model and thus irregularity of the stationary problem.

To see this, we take $s_{0}=0$ and let the initial values for the reduced problem be given by $z^{0}=(z_{1}^{0},z_{2}^{0},0)$ with $z_{0}\in\mathcal{N}(\widehat{A})$ and $\widehat{A}$ denoting the system matrix of the reduced model. Then the solution of the reduced system reads $z(t)=z^{0}$ for all $t\geq 0$, i.e., the damping is completely ineffective and no energy decay takes place. This can already be observed in Table 6.2 for $L=1$. Note that a shift parameter $s_{0}=0$ is the typical choice if one is interested in the long term behavior. The standard reduced model is therefore not exponentially stable for this important case and existence of unique steady states cannot be guaranteed.

7. Numerical tests for the improved reduced model

We now illustrate in more detail the stability and approximation properties of the improved reduced models obtained with the algorithms proposed in the previous sections.

7.1. Mesh independence

As a first example, we consider again the single pipe (TP1). The reduced models are generated for a single input $u_{1}$ at the left vertex $v_{1}$ and we set $s_{0}=0$. In Figure 7.1, we display the basis functions obtained for $L=4$ Krylov iterations.

As explained in detail in Appendix A2, the columns of the matrices $V_{1}$ and $V_{2}$ form orthogonal bases for the subspaces $\mathbb{V}_{1}$ and $\mathbb{V}_{2}$ which can be interpreted as functions on the interval $[0,1]$. In fact, any single Krylov iteration corresponds to the solution of an elliptic boundary value problem which explains why the basis functions are smooth. Also note that the functions look similar to a sequence of orthogonal polynomials of increasing degree. This indicates that the projection onto the Krylov subspaces leads to some sort of higher order approximation.

In the formulation of our algorithms, we payed special attention to a construction that respects the underlying function space setting. As a consequence, the subspaces $\mathbb{V}_{i}$ and even the corresponding bases turn out to be almost independent of the underlying full order model. To illustrate this fact, we display in Figure 7.2 one of the basis functions for velocity and pressure computed with full order models resulting from discretization on different meshes.

Note that the basis functions for different levels coincide almost perfectly, up to discretization errors. This clearly demonstrates the mesh independence of the proposed algorithms. As mentioned before, all algorithms could even be formulated directly for the infinite dimensional problem and the basis functions depicted in Figure 7.2 thus correspond to approximations for the corresponding functions that would be obtained by the Krylov iteration in infinite dimensions.

7.2. Stability of the splitting step

In our basis construction algorithms, we used the cosine-sine decomposition in order to improve the numerical stability of the splitting $W=[W_{1};W_{2}]$. A simple splitting with re-orthogonalization of $W_{1}$ and $W_{2}$, on the other hand, could be realized as follows.

  % function [W1,W2]=simplesplit(W1,W2,M1,M2,tol)
  W1 = ortho(W1,[],M1,tol);
  W2 = ortho(W2,[],M1,tol);

In Figure 7.3, we compare the results obtained by splitting the Krylov bases by this simple strategy with those obtained by means of the cosine-sine decomposition.

Due to the possibility of interpreting the basis vectors as functions on the interval $[0,1]$, one can easily conclude that, already for relatively small dimensions, the standard splitting suffers from severe numerical instabilities. Let us emphasize that this is caused only by the instability of the splitting step and not by the Arnoldi iteration defining the Krylov spaces. The splitting via cosine-sine decomposition, on the other hand, does not suffer from these instabilities and should therefore always be preferred in practice.

7.3. Approximation of the input-output behavior

By Lemma 3.1, we know that the reduced models obtained with our algorithms exactly match the first few moments of the transfer function. This leads to a good overall approximation of the input-output behavior in the frequency domain. With the following tests, we would like to take a closer look also at the approximation in the time domain. For this purpose, we repeat the computations for a single pipe with input functions now given by

[TABLE]

The initial values are set to $p_{0}=0$ and $q_{0}=0$ and we now use the inputs $u_{1}$ and $u_{2}$ at both pipe ends to construct the reduced models. The shift parameter is again set to $s_{0}=0$. For integration of the system in time, we utilize a $\theta$-scheme with uniform time step $\tau$. For $\theta=1$ we obtain the implicit Euler method and for $\theta=\frac{1}{2}+\tau$, the scheme is second order in time. In both cases, the exponential decay is preserved [16].

As can be seen from Figure 7.4, an approximation with only a few moments already leads to a very accurate representation of the input-output behavior in time domain. The oscillations in the initial phase of the output are due to a Gibbs phenomenon. This effect, however, becomes negligible when increasing the dimension of the reduced models. The correct propagation speed and damping of the signal are obtained for all models with $L\geq 5$. For the reduced model with $L=10$ Krylov iterations, which corresponds to $\dim(V_{1})=14$ and $\dim(V_{2})=15$ here, the prediction of the output is almost perfect. The plots displayed in Figure 7.4 also illustrate the exponential decay of the output which becomes faster when the damping is increased.

7.4. Results for a small network

As a final test case, let us now demonstrate that very similar results can also be obtained on more complicated situations. For this purpose we repeat the previous tests for the network depicted in Figure 7.5.

All pipes are chosen to be of unit length $l_{e}=1$ and the model parameters are set constant along every pipe with values

[TABLE]

Here $d_{0}$ is some positive constant that allows us to vary the damping in the whole system by a single factor. These parameters correspond to pipes of different cross-sections; cf. Figure 7.5.

We now repeat the test of the previous section with input defined by (7.1). Since the overall system is substantially larger here, we increase the time horizon by a factor four. As before, we specify a pressure profile at the vertex $v_{1}$ as input and consider as output the resulting mass flux at vertex $v_{2}$ where the pressure is kept at zero. The results are depicted in Figure 7.6.

In comparison to the example with a single pipe, the output function now has a much more complicated structure which is due to multiple pathways through the network and possible reflections at the junctions. Again we observe some Gibbs phenomena for approximations with only a few moments. For $L=20$, which here corresponds to $\dim(V_{1})=40$ and $\dim(V_{2})=47$, we already observe an almost perfect prediction of the input-output behavior. Note that this reduced model amounts to only about $6$ and $7$ degrees of freedom per pipe for pressure and velocity, respectively, which is in good agreement with the experiments for the single pipe.

8. Discussion

Let us briefly summarize the observations made in this paper. Structure-preserving model reduction, as considered for instance in [7, 20, 24, 32, 36], is in principle well suited for the systematic approximation of system of differential-algebraic equations that arise by discretizations of partial differential-algebraic systems modeling of wave propagation phenomena on networks. A proper discretization and subsequent projection onto subspaces allows to preserve the underlying port-Hamiltonian structure and to guarantee passivity of the reduced models. Other important properties, like conservation of mass or exponential stability, are however not inherited automatically. In order to preserve also these properties, some problem specific modifications are required in the subspace construction process. The formulation of appropriate modifications may require a detailed analysis of the underlying mathematical models in infinite dimensions and a detailed understanding of the overall discretization process which in general problem dependent. Therefore, future work has to be devoted to a consideration of further applications, e.g., in elastodynamics or electromagnetics. Another aspect that has to be addressed in future research are nonlinearities in the underlying system.

Acknowledgments

The authors are grateful for financial support by the German Research Foundation (DFG) via grants IRTG 1529 and TRR 154 project B03, C02, C04, as well as by the “Excellence Initiative” of the German Federal and State Governments via the Graduate School of Computational Engineering GSC 233 at Technische Universität Darmstadt.

Appendix

Part III. Functional analytic background

The purpose of this appendix is to show rigorously that the algorithms presented in the previous sections lead to reduced models that satisfy properties (P1)–(P4) uniformly. Most of the following results can in principle be obtained by generalization of those in [17]. To give a complete presentation, we repeat the most important results required for our analysis and provide short proofs where they yield further insight.

Appendix A1 The infinite dimensional problem

Let us start with introducing the relevant notation. We denote by

[TABLE]

the space of square integrable functions over the network with norm

[TABLE]

For convenience of notation, we will sometimes use the symbols $\|\cdot\|_{L^{2}}$ and $\|\cdot\|$ instead. In addition to this basic function space, we will make use of the broken Sobolev space

[TABLE]

consisting of functions that are continuous along edges but may be discontinuous at interior vertices $v\in{\mathcal{V}_{0}}$. The broken derivative of a function $q\in H^{1}(\mathcal{E})$ is denoted by $\partial_{x}^{\prime}q$ defined by

[TABLE]

This allows us to write $H^{1}(\mathcal{E})=\{q\in L^{2}(\mathcal{E}):\partial_{x}^{\prime}q\in L^{2}(\mathcal{E})\}$ with natural norm defined by

[TABLE]

For a piecewise smooth function $q\in H^{1}(\mathcal{E})$, we define at every interior vertex $v\in{\mathcal{V}_{0}}$ the value

[TABLE]

which amounts to the imbalance in the coupling condition (2.3). For a junction of only two pipes, the value $[nq]$ amounts to the jump of $q$ across the junction and the symbol $[nq]$ is the one usually employ in the analysis of discontinuous Galerkin methods. These values of the jumps can be understood as a vector in $\mathbb{R}^{{\mathcal{V}_{0}}}$ and as the scalar product on $\mathbb{R}^{{\mathcal{V}_{0}}}$, we use

[TABLE]

In a similar way, we will denote by $(\lambda,\mu)_{{\mathcal{V}_{\partial}}}$ the corresponding scalar product for $\mathbb{R}^{{\mathcal{V}_{\partial}}}$. We now have the following variational characterization of solutions to our model problem.

Lemma A1.1**.**

Let $(p,q)$ denote a smooth solution of (2.1)–(2.6) and set $\lambda_{v}(t)=p_{v}(t)$. Then

[TABLE]

for all test functions $\tilde{p}\in L^{2}(\mathcal{E})$, $\tilde{q}\in H^{1}(\mathcal{E})$, $\tilde{\lambda}\in\mathbb{R}^{{\mathcal{V}_{0}}}$, and all $t\geq 0$.

Proof.

The proof follows with similar arguments as in [17]. We therefore only sketch the required modifications: The validity of (A1.1) and (A1.3) follows directly from (2.1) and (2.3). Using integration-by-parts on one single edge $e=(v_{1}^{e},v_{2}^{e})$, we get

[TABLE]

Summing over all edges, this yields terms at the inner vertices that can be reordered as

[TABLE]

Using that $p^{e}(v)=\lambda_{v}$ for $v\in{\mathcal{V}_{0}}$ and $p^{e}(v)=u_{v}$ for $v\in{\mathcal{V}_{\partial}}$ and the definition of $[n\tilde{q}]$, this shows the validity of (A1.2) and completes the proof of the lemma. ∎

Remark A1.2** (Well-posedness).**

Together with the initial conditions (2.5), the variational problem (A1.1)–(A1.3) can be shown to admit a unique solution which, for sufficiently regular inputs and initial data, corresponds to the classical solution of (2.1)–(2.6). The variational formulation is thus equivalent to the initial boundary value problem; see [17] for details.

We now proceed by establishing the properties (P1)–(P4) on the continuous level. The results again follow with a slight modification of the arguments given in [17].

A1.1. Conservation of mass

The mass of the fluid contained in a single pipe is given by

[TABLE]

Using the balance equation (2.1) and the conservation condition (2.3), we obtain

Lemma A1.3** (Conservation of mass).**

Let $m(t)=\sum\nolimits_{e\in\mathcal{E}}m^{e}(t)$ denote the total mass contained in the network. Then*

[TABLE]

i.e., the change of mass is caused only by flux across the boundary of the network.

Proof.

Using (2.1), the fundamental theorem of calculus, and (2.4), we get

[TABLE]

The result then follows by using the special form of the output $y_{v}$ given in (2.7). ∎

A1.2. Energy dissipation

Let us first prove property (P1) by deriving an explicit energy dissipation relation. The total acoustic energy contained in a single pipe is given by

[TABLE]

From the differential equations (2.1)–(2.2) and the algebraic continuity conditions (2.3)–(2.4), we can now deduce the following energy dissipation relation.

Lemma A1.4** (Energy dissipation and port-Hamiltonian structure).**

Let $E(t)=\sum\nolimits_{e}E^{e}(t)$ denote the total acoustic energy contained in the network. Then*

[TABLE]

i.e., the change of total energy is caused by power dissipated through the damping mechanism and supplied or drained at the system ports.

Proof.

By elementary calculations and the partial differential equations (2.1)–(2.2), we get

[TABLE]

Integration-by-parts of the second term in the last equation on $e=(v_{1}^{e},v_{2}^{e})$ gives

[TABLE]

Summing over all edges $e$, using the definition of the total energy, and (2.3)–(2.4) leads to

[TABLE]

The result now follows from definition of the in- and output. ∎

A1.3. Exponential stability

Due to linearity of the problem, it suffices to consider the homogeneous case. The energy balance then reveals that kinetic energy is dissipated by the damping mechanism. This, however, also leads to a reduction of the total energy resulting in the exponential stability of the system stated as property (P4).

Lemma A1.5** (Exponential stability).**

Let $u(t)\equiv 0$ for $0\leq t_{1}\leq t\leq t_{2}$, and $E(t)$ be the total energy of the system. Then*

[TABLE]

with positive constants $C,\gamma>0$ that are independent of $t_{1}$, $t_{2}$, $s$, and $t$.

Proof.

The proof is based on energy estimates, some graph theoretic results, and a generalized Poincaré inequality; we refer to [17] for details. ∎

A1.4. Steady states

From the previous result, we obtain convergence to zero steady state in case of homogeneous input $u\equiv 0$. Due to the linearity of the problem, this yields also the existence of unique and stable steady states in the general case.

Lemma A1.6** (Steady states).**

Let $u(t)\equiv const$ for all $t\geq t_{1}$. Then $(p(t),q(t))$ converges to a steady state $(\bar{p},\bar{q})$, which is the unique solution of the corresponding stationary problem.

Proof.

The existence of a unique steady state has been established in [17]. The difference to steady state $(p(t)-\bar{p},q(t)-\bar{q})$ solves (2.1)–(2.5) with $u\equiv 0$ and convergence to steady state thus follows by the exponential stability estimate given in the previous lemma. ∎

Appendix A2 Galerkin approximation

We now extend the discretization strategy proposed in [17] to our setting and review the basic results about the stability of these full order models. Let $P_{h}\subset L^{2}(\mathcal{E})$ and $Q_{h}\subset H^{1}(\mathcal{E})$ be finite dimensional spaces and set $\Lambda_{h}=\mathbb{R}^{{\mathcal{V}_{0}}}$. Let $T>0$ and consider the following conforming Galerkin approximations of the variational principle (A1.1)–(A1.3) as space discretization.

Problem A2.1** (Galerkin approximation and discrete variational problem).**

Find $p_{h}\in H^{1}(0,T;P_{h})$, $q_{h}\in H^{1}(0,T;Q_{h})$, and $\lambda_{h}\in L^{2}(0,T;\Lambda_{h})$ such that

[TABLE]

for all $\tilde{p}_{h}\in P_{h}$ and $\tilde{q}_{h}\in Q_{h}$, and such that the discrete variational equations

[TABLE]

hold for all test functions $\tilde{p}_{h}\in P_{h}$, $\tilde{q}_{h}\in Q_{h}$, $\tilde{\lambda}_{h}\in\Lambda_{h}$, and all $0\leq t\leq T$.

A simple compatibility condition allows to deduce the well-posedness of this problem.

Lemma A2.2** (Discrete well-posedness).**

Assume that $\{1^{e}:e\in\mathcal{E}_{h}\}\subset Q_{h}$, where $1^{e}$ denotes the function in $L^{2}(\mathcal{E})$ which is constant one on the edge $e$ and zero otherwise. Then Problem A2.1 has a unique solution.

Proof.

The condition $\{1^{e}:e\in\mathcal{E}\}\subset Q_{h}$ allows to eliminate the constraint and the Lagrange multiplier $\lambda$ from the system; compare with Section 6.2. By choosing bases for the spaces $P_{h}$ and $Q_{h}$, we may thus obtain linear system of ordinary differential equations. Existence and uniqueness then follow from the Picard-Lindelöf theorem. ∎

Remark A2.3**.**

Once the solution of the discretized problem is found, we can define the corresponding output of the discrete system by $y_{h}(v)=-n^{e}(v)q^{e}_{h}(v)$ for $v\in{\mathcal{V}_{\partial}}$ and $e\in\mathcal{E}(v)$. Let us note that elimination of the constraints (A2.3) yields the problem originally considered in [17]. All results obtained in that paper therefore carry over to the problem considered here, if the condition of the previous Lemma is satisfied, which is called assumption (A3h) below.

We next establish the properties (P1)–(P4) for the Galerkin approximations introduced above. The results again follow with similar arguments as used in [17]. Let us emphasize that additional conditions on the approximation spaces $P_{h}$ and $Q_{h}$ are required for some of the results.

A2.1. Conservation of mass

Mass conservation on the continuous level follows by testing (A1.1) with the function $\tilde{p}\equiv 1$ and some elementary manipulations. A discrete equivalent of this result can be obtained, if $\tilde{p}_{h}\equiv 1$ is contained in the test space.

Lemma A2.4** (Discrete mass conservation).**

Let $m_{h}(t)=\sum\nolimits_{e\in\mathcal{E}}\int_{e}ap_{h}(t)dx$ denote the total mass of the discrete system, and let

• (A1h)

$1\in P_{h}$.

Then the total mass changes only due to flux across the boundary, i.e.,

[TABLE]

Proof.

The assertion follows in the same manner as that of Lemma A1.3. ∎

A2.2. Energy balance

Mimicking the notation used on the continuous level, we may define for every pipe $e\in\mathcal{E}$ the total (discrete) energy content of the pipe by

[TABLE]

With the same arguments as on the continuous level, we then obtain

Lemma A2.5** (Discrete energy balance and port-Hamiltonian structure).**

Let $E_{h}(t)=\sum\nolimits_{e\in\mathcal{E}}E_{h}^{e}(t)$ denote the total discrete energy. Then

[TABLE]

i.e., the energy changes by dissipation and supply or drain via the ports of the network.

Note that no extra condition for the approximation spaces is required for the proof of property (P1) for the Galerkin approximations of the variational formulation (A1.1)–(A1.3).

A2.3. Exponential stability

For input $u\equiv 0$, the energy balance (A2.4) already guarantees that the total energy of the discrete system is non-increasing. Under additional compatibility conditions on the approximation spaces, one can even show the uniform exponential decay.

Lemma A2.6** (Uniform discrete exponential stability).**

Assume that

• (A2h)

$\partial_{x}^{\prime}Q_{h}=P_{h}$;

• (A3h)

$1^{e}\in Q_{h}$* for all $e\in\mathcal{E}$.*

Then for $u\equiv 0$, the discrete energy $E_{h}(t)$ defined in Lemma A2.5 satisfies

[TABLE]

Moreover, the constants $C,\gamma$ can be chosen the same as on the continuous level.

Proof.

Following Remark A2.3, the proof can be deduced from the results of [17]. ∎

A2.4. Steady states

Under the assumptions of the previous lemma, one can also guarantee property (P4), i.e., the existence and uniqueness of discrete steady states.

Lemma A2.7** (Discrete steady states).**

Let $u(t)\equiv const$ for $t\geq t_{0}$ and assume that (A2h)–(A3h) hold. Then for $t\to\infty$, the discrete solution $(p_{h}(t),q_{h}(t))$ converges to a discrete equilibrium $(\bar{p}_{h},\bar{q}_{h})$ which is the unique solution of the corresponding stationary problem.

Proof.

The existence of a unique discrete steady state is established in [17]. Convergence to equilibrium then follows from the energy decay estimate like in Lemma A1.6. ∎

A2.5. A mixed finite element approximation

As a particular Galerkin approximation satisfying the above assumptions, let us briefly discuss the mixed finite element method that is used in our numerical tests. Let $[0,l^{e}]$ be the interval represented by the edge $e$ and denote by $T_{h}(e)=\{T\}$ a uniform mesh of $e$ with subintervals $T$ of length $h^{e}$. The global mesh is then defined as $T_{h}(\mathcal{E})=\{T_{h}(e):e\in\mathcal{E}\}$, and the global mesh size is denoted by $h=\max_{e}h^{e}$. We denote the spaces of piecewise polynomials on $T_{h}(\mathcal{E})$ by

[TABLE]

where $P_{k}(T_{h}(e))=\{v\in L^{2}(e):v|_{T}\in P_{k}(T),\ T\in T_{h}(e)\}$ and $P_{k}(T)$ is the space of polynomials of degree $\leq k$ on the subinterval $T$. Note that $P_{k}(T_{h}(\mathcal{E}))\subset L^{2}(\mathcal{E})$, which is easy to see, but in general $P_{k}(T_{h}(\mathcal{E}))\not\subset H^{1}(\mathcal{E})$. As spaces $V_{h}$ and $Q_{h}$ for the Galerkin approximation presented in the previous sections, we now consider

[TABLE]

This choice of spaces satisfies the compatibility conditions (A1h)–(A3h); see [17] for details.

A2.6. Structure preserving model reduction

As final result of this section, we now give an interpretation of the model reduction approach on the level of function spaces. The system obtained by Galerkin projection onto $P_{h}$, $Q_{h}$ will again be called full order model. The reduced models are obtained by projection onto smaller subspaces $P_{H}\subset P_{h}$ and $Q_{H}\subset Q_{h}$. The compatibility conditions for these coarse subspaces read

• (A1H)

$1\in P_{H}$;

• (A2H)

$\partial_{x}^{\prime}Q_{H}=P_{H}$;

• (A3H)

$1^{e}\in Q_{H}$ for all $e\in\mathcal{E}$.

Note that by construction $P_{H}\subset P_{h}\subset L^{2}(\mathcal{E})$ and $Q_{H}\subset Q_{h}\subset H^{1}(\mathcal{E})$. From the previous results about general Galerkin approximations, we therefore directly deduce the following result.

Lemma A2.8** (Structure preserving model reduction).**

Let $P_{H}\subset P_{h}$ and $Q_{H}\subset Q_{h}$ and assume that (A1H)–(A3H) hold. Then the reduced system satisfies (P1)–(P4) and the assertions of Lemma A2.5–A2.7 hold accordingly.

Remark A2.9**.**

Since the reduced model can be viewed as Galerkin approximation of the infinite dimensional problem (A1.1)–(A1.3), it is clear that the solution $(p_{H}(t),q_{H}(t))$ only depends on the choice of the approximation spaces $P_{H}$ and $Q_{H}$ but not on the spaces $P_{h}$ and $Q_{h}$ used for the generation of the full order model which is only required for computational purposes and has no effect on the quality of the reduced model.

Appendix A3 Reformulation on the algebraic level

We now translate the results of the previous section to the algebraic level. By choosing appropriate bases for the subspaces $P_{h}$ and $Q_{h}$ defining the discrete variational problem (A2.1)–(A2.3), the resulting full order model can be written in algebraic form as follows.

Lemma A3.1** (Equivalent algebraic system).**

Let $\{\phi_{i}\}$ and $\{\psi_{j}\}$ be bases for $P_{h}$ and $Q_{h}$. Then the problem (A2.1)–(A2.3) is equivalent to the system (2.8)–(2.10) with matrices defined by $M_{1}(i,j)=(a\phi_{j},\phi_{i})_{\mathcal{E}}$, $M_{2}(i,j)=(b\psi_{j},\psi_{i})_{\mathcal{E}}$, $G(i,j)=(\partial_{x}^{\prime}\psi_{j},\phi_{i})_{\mathcal{E}}$, $D(i,j)=(d\psi_{j},\psi_{i})_{\mathcal{E}}$, $N(i,j)=[n\psi_{j}](v_{i_{0}})$, and $B_{2}(i,j)=-\delta_{j_{\partial},i}n(v_{j})\psi_{i}(v_{j})$.

Here $(a,b)_{\mathcal{E}}=\sum\nolimits_{e}\int_{e}a(x)b(x)dx$ is the scalar product of $L^{2}(\mathcal{E})$, and $1\leq i_{0}\leq|{\mathcal{V}_{0}}|$ and $1\leq i_{\partial}\leq|{\mathcal{V}_{\partial}}|$ denote the appropriate renumbering of inner and boundary vertices. From the definition of the matrices, we also obtain the following properties.

Lemma A3.2**.**

Let $a^{e}$, $b^{e}$, and $d^{e}$ be positive for all $e\in\mathcal{E}$. Then $M_{1}$, $M_{2}$, $D$ are symmetric and positive definite. If (A2h)–(A3h) hold, then $[G^{\top},N^{\top}]$ is injective and $N$ is surjective on $\mathcal{N}(G)$.

These properties allow us to establish the well-posedness of the linear time-invariant system (2.8)–(2.10) and, following our discussion in Section 3.3, also the unique solvability of the corresponding stationary problem. As a next step, we can now provide an interpretation of the properties (P1)–(P4) on the algebraic level.

A3.1. Conservation of mass

The condition $1\in P_{h}$ is equivalent to

• (A1${}_{h}^{\prime}$)

$\exists o_{1}\in\mathbb{R}^{k_{1}}:\sum\nolimits_{i=1}^{k_{1}}o_{1,i}\phi_{i}=1$ on $\mathcal{E}$.

Here $o_{1}\in\mathbb{R}^{k_{1}}$ is the coordinate vector representing the function $1\in P_{h}$ in the basis $\{\phi_{i}\}$ of the space $P_{h}$. This allows us to express the conservation of mass on the algebraic level as follows.

Lemma A3.3** (Mass conservation).**

Let $p_{h}(\cdot,t)=\sum\nolimits_{i=1}^{k_{1}}x_{1,i}(t)\phi_{i}(\cdot)\subset P_{h}$ and further define $q_{h}(\cdot,t)=\sum\nolimits_{j=1}^{k_{2}}x_{2,j}(t)\psi_{j}(\cdot)\subset Q_{h}$. Then the mass of the discrete system can be expressed as

[TABLE]

With $\hat{o}\in\mathbb{R}^{{\mathcal{V}_{\partial}}}$ denoting the constant one vector and $y=B_{2}^{\top}x_{2}$ the output, we have

[TABLE]

Proof.

The definition of the total mass and (2.8) lead to

[TABLE]

The fact that $o^{\top}_{1}Gx_{2}=\hat{o}^{\top}Bx_{2}=\hat{o}^{\top}y$ can be deduced from the equivalent formulation of the Galerkin approximation in function spaces; cf. Lemma A2.4. ∎

Again, the definition of the total mass and the proof of the mass conservation requires awareness of the underlying problem in function spaces.

A3.2. Energy balance

The energy of the discrete system can be expressed as

[TABLE]

where $(x_{1}(t),x_{2}(t),x_{3}(t))$ is a solution of (2.8)–(2.10) and $p_{h}(\cdot,t)=\sum\nolimits_{i=1}^{k_{1}}x_{1,i}(t)\phi_{i}(\cdot)\subset P_{h}$, $q_{h}(\cdot,t)=\sum\nolimits_{j=1}^{k_{2}}x_{2,j}(t)\phi_{j}(\cdot)$, and $\lambda_{h}(t)=x_{3}(t)$ define the corresponding functions making up the solution of of Problem A2.1. The discrete energy balance of Lemma A2.5 can now be rephrased as

Lemma A3.4** (Energy dissipation and port-Hamiltonian structure).**

Let $(x_{1},x_{2},x_{3})$ denote a solution of the linear system (2.8)–(2.10). Then

[TABLE]

with output defined as $y(t)=B^{\top}x_{2}(t)$. In particular, property (P1) is valid.

Proof.

Let us give a direct derivation of this assertion on the algebraic level. Using the definition of the energy, the symmetry of $M_{i}$, and the algebraic equations, we obtain

[TABLE]

In the last step, we utilized that $x_{2}^{\top}Bu=(B^{\top}x_{2})^{\top}u$ and the definition of the output. ∎

Remark A3.5**.**

As can be seen from the proof, the scalar product and norm induced by the matrices $M_{1}$ and $M_{2}$ are directly associated with the energy of the system and they are the natural ones for the analysis and numerical treatment of the discrete problem in algebraic form.

As shown above, the property (P1) follows directly from the particular form of the algebraic system. As will become clear below, the validity of the remaining properties (P2)–(P4) however requires awareness of the underlying infinite dimensional problem.

A3.3. Exponential stability

In order to translate the assertions of Lemma A2.6 to the algebraic level, we have to describe the meaning of the conditions (A2h)–(A3h).

Lemma A3.6**.**

The compatibility conditions (A2)–(A3) are equivalent to

• (A2${}_{h}^{\prime}$)

for all $x_{1}\in\mathbb{R}^{k_{1}}$ there exists $x_{2}\in\mathbb{R}^{k_{2}}$ such that $\sum\nolimits_{i=1}^{k_{1}}x_{1,i}\phi_{i}=\sum\nolimits_{j=1}^{k_{2}}x_{2,j}\psi_{j}$;

• (A3${}_{h}^{\prime}$)

for all $e\in\mathcal{E}$ there exists $o_{2}^{e}\in\mathbb{R}^{k_{2}}$ such that $\sum\nolimits_{j=1}^{k_{2}}o^{e}_{2,j}\psi_{j}=1^{e}\in Q_{h}$.

As a direct consequence of this characterization, Lemma A2.6, and the equivalence of the algebraic system to the Galerkin approximation, we obtain the following lemma.

Lemma A3.7** (Exponential stability).**

Let (A2${}_{h}^{\prime}$)–(A3${}_{h}^{\prime}$) hold and let $(x_{1}(t),x_{2}(t),x_{3}(t))$ be a solution of (2.8)–(2.10). Then

[TABLE]

with constants $C,\gamma>0$ that can be chosen as in Lemma A1.5 and A2.6.

The derivation of the conditions (A2${}_{h}^{\prime}$)–(A3${}_{h}^{\prime}$) and the proof of the stability estimate can thus be deduced from the underlying Galerkin approximation and the analysis in function spaces.

A3.4. Steady states

The assertion about steady states can finally be translated as follows.

Lemma A3.8** (Steady states).**

Let (A2${}_{h}^{\prime}$)–(A3${}_{h}^{\prime}$) hold. Then for $u(t)\equiv const$, the solutions $(x_{1}(t),x_{2}(t),x_{3}(t))$ of the system (2.8)–(2.10) converge to a steady state $(\bar{x}_{1},\bar{x}_{2},\bar{x}_{3})$ which is the unique solution of the corresponding stationary problem.

Proof.

The result follows again by equivalence to the Galerkin approximation (A2.1)–(A2.3) and the corresponding result stated in Lemma A2.7. ∎

A3.5. Structure preserving model reduction

As a final step of our analysis, we can now provide a proof for Theorem 5.5 by translating the results of Section A2.6 to the algebraic level: Let $\{\phi_{i}\}$ and $\{\psi_{j}\}$ denote bases for $P_{h}$ and $Q_{h}$, and let $V_{1}\in\mathbb{R}^{k_{1}\times K_{1}}$ and $V_{2}\in\mathbb{R}^{k_{2}\times K_{2}}$ be given matrices with linearly independent columns. For $k=1,\ldots,K_{1}$ and $l=1,\ldots,K_{2}$, we define

[TABLE]

which serve as basis functions for low dimensional approximation spaces $P_{H}$ and $Q_{H}$. As a direct consequence of this construction and the previous considerations, we obtain

Lemma A3.9** (Reduced algebraic system).**

Let us define

[TABLE]

Then $P_{H}\subset P_{h}$ and $Q_{H}\subset Q_{h}$, and the corresponding discrete variational problem is equivalent to the reduced order system (2.14)–(2.16) with matrices as defined in Section 2.4.

The property (P1) for the reduced system (2.14)–(2.16) again follows directly from the special algebraic form of the reduced problem. In order to guarantee (P2)–(P4), we require additional compatibility conditions. The following characterization clarifies the picture.

Lemma A3.10** (Algebraic compatibility conditions).**

Let (A1${}_{h}^{\prime}$)–(A3${}_{h}^{\prime}$) hold and assume that the algebraic compatibility conditions (A1)–(A3) are valid. Then $P_{H}\subset P_{h}$ and $Q_{H}\subset Q_{h}$ satisfy the compatibility conditions (A1H)–(A3H).

As a direct consequence of the previous considerations, we now obtain the following result.

Lemma A3.11** (Structure preserving model reduction).**

Let (A1${}_{h}^{\prime}$)–(A3${}_{h}^{\prime}$) hold for the system (2.8)–(2.10) and assume that the algebraic conditions (A1)–(A3) are valid. Then the reduced system (2.14)–(2.16) satisfies (P1)–(P4).

This lemma yields a correct statement of Theorem 5.5 and completes the proof of our assertions.