A Rosenbrock framework for tangential interpolation of port-Hamiltonian descriptor systems

ABSTRACT

We present a new structure-preserving model order reduction (MOR) framework for large-scale port-Hamiltonian descriptor systems (pH-DAEs). Our method exploits the structural properties of the Rosenbrock system matrix for this system class and utilizes condensed forms which often arise in applications and reveal the solution behaviour of a system. Provided that the original system has such a form, our method produces reduced-order models (ROMs) of minimal dimension, which tangentially interpolate the original model’s transfer function and are guaranteed to be again in pH-DAE form. This allows the ROM to be safely coupled with other dynamical systems when modelling large system networks, which is useful, for instance, in electric circuit simulation.

KEYWORDS:

• Port-Hamiltonian systems
• differential-algebraic systems
• structure-preserving model reduction

1. Introduction

The port-Hamiltonian (pH) modelling paradigm provides an energy-based framework for constructing high-fidelity models of complex dynamical systems. The separation between constitutive relations and the interconnection structure enables a modular modelling approach, in which different subsystems are modelled independently and then interconnected via power flows while preserving important physical properties [1]. This flexibility is particularly advantageous when considering interactions between subsystems across different physical domains or time scales [2]. Furthermore, network laws that govern these interactions, such as Kirchhoff’s laws in electrical circuits or position and velocity constraints in mechanical systems, may be incorporated via algebraic constraints, leading to pH differential-algebraic equation systems (pH-DAEs) [3]. PH-DAE systems often naturally emerge in modelling practice, e.g. when modelling electric circuits, incompressible fluid flow, multibody dynamics or pressure waves in gas networks (see [3,26] and the references therein).

If the system at hand increases in complexity, generally, the state-space dimension of its associated pH-DAE model does so too, which makes the simulation and design of model-based controllers computationally challenging. Model order reduction (MOR) may be applied to approximate large-scale full-order models (FOMs) with reduced-order models (ROMs) of substantially smaller dimension to decrease the computational costs. However, this comes with two major challenges. On the one hand, the impacts of the algebraic equations on the model dynamics have to be reflected in the ROM as well, but without significantly increasing the reduced state-space dimension. On the other hand, to facilitate the subsequent coupling of the ROM with other subsystems, the MOR method should be structure-preserving, meaning that the ROM is again a pH-DAE.

Structure-preserving MOR methods for pH-DAEs fall into three major categories: pH-preserving MOR methods, passivity-preserving MOR methods and passivity enforcement techniques. For an overview of existing approaches, the interested reader is referred to [5]. In this work, we focus on pH-preserving, interpolatory MOR techniques which are particularly suited for the reduction of very large-scale models due to their computational efficiency. While these methods are well-established for various system classes, including general DAE systems (see [6] for a comprehensive overview), this is still only partially the case for pH-DAEs. Tangential interpolation of pH ordinary differential equation systems (pH-ODEs) with no algebraic constraints has been addressed in [8–10,11,22] and is well understood. Extensions to the DAE case with algebraic constraints have been proposed in [12,13]. The strategies to deal with algebraic constraints vary depending on the original model’s Kronecker index and how the algebraic constraints are affected by the input. However, in their entirety, these strategies do not cover the whole system class of pH-DAEs and some may not necessarily preserve the pH structure. This motivates the following research question which we address in this work: Is there a unifying interpolatory MOR framework that is applicable to the entire system class of linear, time-invariant pH-DAEs and guarantees a preservation of the pH structure?

A good starting point to answer this question is the staircase forms presented in [13,14]. While staircase forms are generally helpful in analysing the solution behaviour of DAEs, they are also useful for MOR since they allow the impact of algebraic constraints on the system dynamics to be analysed. In general, the transformation of a pH-DAE model to staircase form relies on a series of rank conditions which may be sensitive under perturbations; see e.g [15]. Fortunately, as demonstrated in [3,26] for various physical examples, this can be considered directly in the modelling process such that the resulting model already has (or is close to) staircase form. Since DAEs often include redundant algebraic equations which increase the computational cost for simulation, the task of MOR is to identify a minimal set of equations that describes the dynamical behaviour of the system. In 1970, Rosenbrock [16] proposed to represent dynamical systems with a polynomial matrix, which is beneficial to compute minimal representations and nowadays known as the Rosenbrock system matrix. Interestingly enough, this matrix has not extensively been used in the context of MOR.

In this work, we show that the Rosenbrock system matrix exhibits a particular structure for pH-DAEs which we exploit to derive a novel interpolatory MOR framework that

(i) can be used for linear, time-invariant pH-DAEs in staircase form with arbitrary Kronecker index,

(ii) guarantees ROMs in pH-DAE form with a minimal state-space dimension,

(iii) works with the original (typically sparse) state-space matrices and is therefore computationally efficient, and

(iv) enables a straightforward adaptation and integration of different interpolatory MOR strategies that have been developed for general dynamical systems.

The paper is organized as follows. In Section 2, we use the concept of the Rosenbrock system matrix to introduce the system class of pH-DAEs and its staircase forms and recapitulate the basics of tangential interpolation. In Section 3, we propose a new interpolatory MOR framework that exploits the structural properties of pH-DAEs in staircase form and show how this framework naturally generalizes to several extensions proposed for unstructured systems in Section 4. We conclude with an application of our method to two electric circuits and a discussion of the above claims in Sections 5 and 6.

2. Preliminaries

We consider linear time-invariant (LTI) systems of the form

(1) $\begin{array}{cc}E\dot{x}(t)& =Ax(t)+Bu(t),x(0)=0,\\ y(t)& =Cx(t)+Du(t),\end{array}$(1)

with state vector $x(t)\in {\mathbb{R}}^n$, inputs $u(t)\in {\mathbb{R}}^m$, outputs $y(t)\in {\mathbb{R}}^m$ for all $t\in [0,\mathrm{\infty })$ and constant matrices $E,A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times m}$, $C\in {\mathbb{R}}^{m\times n}$ and $D\in {\mathbb{R}}^{m\times m}$. Systems with a singular descriptor matrix $E$ are referred to as DAE systems, and ODE systems otherwise. In the following, we will assume that the pencil $\lambda E-A$ is regular, i.e. $det(\lambda E-A)\ne 0$ for some $\lambda \in \mathbb{C}$. We denote the ring of polynomials with coefficients in $\mathbb{R}$ by $\mathbb{R}[s]$ and the set of $n\times m$ matrices with entries in $\mathbb{R}[s]$ by $\mathbb{R}[s{]}^{n\times m}$.

2.1. The Rosenbrock system matrix

After a Laplace transformation of the state-space equations in (1), we obtain the following equations in matrix form:

(2) $\mathcal{P}\left(s\right)\left[\begin{array}{c}\xi \left(s\right)\\ \mu \left(s\right)\end{array}\right]=\left[\begin{array}{c}0\\ \psi \left(s\right)\end{array}\right],\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathcal{P}\left(s\right):=\left[\begin{array}{cc}sE-A& -B\\ C& D\end{array}\right]\in \mathbb{R}[s{]}^{\left(n+m\right)\times \left(n+m\right)}.$(2)

where $\xi$, $\mu$ and $\psi$ denote the Laplace transforms of the state, input and output vector, respectively. The polynomial matrix $\mathcal{P}$ is also called the Rosenbrock system matrix (in the following: system matrix) [17]. Assuming regularity, the linear mapping between inputs $\mu$ and outputs $\psi$ that follows from (2) is given by the system’s transfer function

(3) $H(s):=C(sE-A{)}^{-1}B+D.$(3)

Since all transformations of the state-space equations in (1) can be expressed by operations on $\mathcal{P}$, the system matrix proved to be particularly useful for studying the properties of these transformations [16]. In this work, we shall focus on transformations which leave the system’s transfer function $H$ unchanged and are summarized by the notion of strict system equivalence.

Lemma 2.1.

[16] Let $X(s)\in \mathbb{R}[s{]}^{m\times n}$, $Y(s)\in \mathbb{R}[s{]}^{n\times m}$ and define unimodular matrices $L(s),M(s)\in \mathbb{R}[s{]}^{n\times n}$, i.e. their determinants are nonzero constants. Suppose that two system matrices $\mathcal{P}$ and $\tilde{\mathcal{P}}$ are related by the transformation

(4) $\tilde{\mathcal{P}}(s)={\mathcal{T}}_1(s)\mathcal{P}(s){\mathcal{T}}_2(s)=\left[\begin{array}{cc}L(s)& 0\\ X(s)& I_m\end{array}\right]\mathcal{P}(s)\left[\begin{array}{cc}M(s)& Y(s)\\ 0& I_m\end{array}\right].$(4)

Then we shall say that $\mathcal{P}$ and $\tilde{\mathcal{P}}$ are related by strict system equivalence (s.s.e.). The two system matrices give rise to the same transfer function, i.e. $H(s)=\tilde{H}(s)$ for all $s\in \mathbb{C}$.

Proof.

For a proof, we refer the reader to [16, Section 3.1].

The benefits of representing a dynamical system with its system matrix are by no means restricted to system transformations; it may also be useful in the context of model reduction.

2.2. Interpolatory model reduction

The goal of MOR is to find a reduced-order model

(5) $\begin{array}{cc}{E}_r{\dot{x}}_r(t)& =A_r{x}_r(t)+B_{r}u(t),x_r(0)=0,\\ y_r(t)& =C_r{x}_r(t)+D_{r}u(t),\end{array}$(5)

with reduced state vector $x_r(t)\in {\mathbb{R}}^r$ and an associated transfer function $H_r$ such that $r\ll n$ and $y\approx y_r$ for certain $u$. In projection-based MOR, these models are created by means of Petrov-Galerkin projections. Here, we define two matrices $U,V\in {\mathbb{R}}^{n\times r}$ with full column rank. The original state trajectory $x(t)$ is approximated on the column space of $V$, i.e. $x(t)\approx V{x}_r(t)$ for all $t\in [0,\mathrm{\infty })$. The reduced system matrix is obtained by the following operation on the original system matrix $\mathcal{P}$:

(6) ${\mathcal{P}}_r(s)=\left[\begin{array}{cc}{U}^{\mathrm{\top }}& 0\\ 0& I_m\end{array}\right]\mathcal{P}\left[\begin{array}{cc}V& 0\\ 0& I_m\end{array}\right]=\left[\begin{array}{cc}s{U}^{\mathrm{\top }}EV-U^{\mathrm{\top }}AV& -U^{\mathrm{\top }}B\\ CV& D\end{array}\right]=\left[\begin{array}{cc}s{E}_r-A_r& -B_r\\ C_r& D_r\end{array}\right].$(6)

The matrices $U,V$ may be chosen for different purposes, for example, to enforce the preservation of certain system-theoretic or structural properties of the original model. In interpolatory MOR methods, they are used to enforce interpolation conditions between the original and reduced transfer function. In particular, the column space of $V$ may be chosen such that the reduced transfer function $H_r$ tangentially interpolates the original transfer function $H$ at a set of interpolation points or shifts $\{{\sigma }_1,\dots ,{\sigma }_r\}\in \mathbb{C}$:

(7) $H({\sigma }_i)b_i=H_r({\sigma }_i)b_i,i=1,\dots ,r,$(7)

where $b_i\in {\mathbb{C}}^m$ denotes the associated (right) tangential direction. For the sake of simplicity, we assume throughout this work that each shift is distinct while an extension to multiple shifts is straightforward, and we refer the interested reader to [6] for details. In all cases, we require ${\sigma }_{i}E-A$ as well as ${\sigma }_i{E}_r-A_r$ to be nonsingular for all $i=1,\dots ,r$, and that the sets of shifts and associated tangential directions are both closed under complex conjugation. Then, the interpolation conditions in (7) may be enforced by choosing

(8) ${\left({\sigma }_{i}E-A\right)}^{-1}B{b}_i\in range\left(V\right),i=1,\dots ,r,$(8)

where $range\left(V\right)$ denotes the associated column space of $V$. Similarly, $U$ may be used to enforce additional (left) tangential interpolation conditions (see [6] for details). The approximation quality of a ROM is typically assessed by computing the ${\mathcal{H}}_{\mathrm{\infty }}$ or ${\mathcal{H}}_2$ norm of the error function $H-H_r$. These are defined in the Hardy spaces ${\mathcal{R}\mathcal{H}}_{\mathrm{\infty }}^{m\times m}$ (${\mathcal{R}\mathcal{H}}_2^{m\times m}$) of all proper (strictly proper) real-rational $m\times m$ matrices without poles in the closed complex right half-plane (see, e. g., [18] for details).

Remark 1. Note that while the conditions in (8) hold regardless of whether $E$ has full rank or not, as long as ${\sigma }_{i}E-A$ and ${\sigma }_i{E}_r-A_r$ are nonsingular for all $i=1,\dots ,r$, the interpolation of DAE systems poses additional challenges. Unlike in the ODE case, where $H-H_r$ is guaranteed to be strictly proper since ${lim}_{s\to \mathrm{\infty }}H(s)-H_r(s)=D-D_r=0$, this is not necessarily the case for DAE systems. In fact, as we will also see in Section 3, the algebraic constraints may even lead to improper transfer functions with ${lim}_{s\to \mathrm{\infty }}H(s)=\mathrm{\infty }$. It is, therefore, crucial to analyse the impact of the algebraic constraints on $H$ because otherwise, $H-H_r$ may grow unboundedly large for $s\to \mathrm{\infty }$, and the error norms are no longer defined. For an overview of how to deal with this challenge for general LTI systems, the reader is referred to [6, Chapter 9].

2.3. Port-Hamiltonian descriptor systems

The system class of linear port-Hamiltonian descriptor systems with quadratic Hamiltonian was introduced in [19]. In this work, we focus on the subclass of constant-coefficient pH-DAEs, which can be characterized by a special structure of the associated system matrix.

Definition 2.2.

Consider a regular LTI system of the form

(9) $\mathrm{\Sigma }:\left\{\begin{array}{cc}E\dot{x}\left(t\right)& =\left(J-R\right)x\left(t\right)+\left(G-P\right)u\left(t\right),x\left(0\right)=0,\\ y\left(t\right)& =(G+P{)}^{\mathrm{\top }}x\left(t\right)+\left(S+N\right)u\left(t\right),\end{array}\right.$(9)

where $E,J,R\in {\mathbb{R}}^{n\times n}$, $G,P\in {\mathbb{R}}^{n\times m}$, $S,N\in {\mathbb{R}}^{m\times m}$. We call the system a pH-DAE system if its system matrix may be decomposed into the following sum of symmetric and skew-symmetric parts

(10) $\mathcal{P}(s)=s\underset{=:\mathcal{E}}{\underbrace{\left[\begin{array}{cc}E& 0\\ 0& 0\end{array}\right]}}+\underset{=:\mathrm{\Gamma }}{\underbrace{\left[\begin{array}{cc}-J& -G\\ G^{\mathrm{\top }}& N\end{array}\right]}}+\underset{=:W}{\underbrace{\left[\begin{array}{cc}R& P\\ P^{\mathrm{\top }}& S\end{array}\right]}},$(10)

such that

1. the structure matrix $\mathrm{\Gamma }$ is skew-symmetric, i.e. $\mathrm{\Gamma }=-{\mathrm{\Gamma }}^{\mathrm{\top }}$,

2. the dissipation matrix $W$ is symmetric positive semi-definite, i.e. $W=W^{\mathrm{\top }}\ge 0$,

(iii) the extended descriptor matrix $\mathcal{E}$ is symmetric positive semi-definite, i.e. $\mathcal{E}={\mathcal{E}}^{\mathrm{\top }}\ge 0.$

The system has an associated quadratic Hamiltonian $\mathcal{H}(x)=\frac{1}{2}{x}^{\mathrm{\top }}Ex$ and transfer function

$H(s):=(G+P{)}^{\mathrm{\top }}(sE-(J-R){)}^{-1}(G-P)+S+N.$

Note that the definition proposed in [19] appears to be more general since it also allows the representation of systems governed by a quadratic Hamiltonian of the form $\mathcal{H}(x)=\frac{1}{2}{x}^{\mathrm{\top }}{Q}^{\mathrm{\top }}Ex$ with $Q\in {\mathbb{R}}^{n\times n}$ and $Q^{\mathrm{\top }}E=E^{\mathrm{\top }}Q\ge 0$. However, our definition does not impose any additional restrictions since it has been shown, e.g. in [3] that every pH-DAE, as defined in [19], may be reformulated to have the form in (9).

Moreover, it was shown in [13,14] that every pH-DAE may be transformed into staircase form. While the physical interpretation of the states is generally lost during this transformation, staircase forms are useful to study the solution behaviour of pH-DAEs [3] and, as we will show in Section 3, also simplify model reduction. We refer to [14, Algorithm 5] for details on how to perform this transformation. Before we proceed, let us highlight that it may require several subsequent full rank decompositions, which are sensitive to perturbations. Fortunately, due to the structural properties of the system matrix, the number of required decompositions is limited to three in contrast to general (unstructured) DAE systems. Moreover, if this condensed form is directly considered during modelling, fewer steps are required, as discussed in [3,26]. In some practical cases, for example, in the modelling of electric circuits, the staircase form even naturally arises or can be enforced by simple structure-preserving permutations of the system equations, as illustrated in Section 5.

Lemma 2.3.

[13,14] A regular pH-DAE system is in staircase form if it has a partitionedstate vector $x(t)={\left[x_1{(t)}^{\mathrm{\top }},x_2{(t)}^{\mathrm{\top }},x_3{(t)}^{\mathrm{\top }},x_4{(t)}^{\mathrm{\top }}\right]}^{\mathrm{\top }}$, where $x_j(t)\in {\mathbb{R}}^{n_j},n_j\in {\mathbb{N}}_0$ for all $j=1,\dots ,4$ such that

(11) $E:=\left[\begin{array}{cccc}{E}_{11}& 0& 0& 0\\ 0& E_{22}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],J:=\left[\begin{array}{cccc}{J}_{11}& J_{12}& J_{13}& J_{14}\\ J_{21}& J_{22}& J_{23}& 0\\ J_{31}& J_{32}& J_{33}& 0\\ J_{41}& 0& 0& 0\end{array}\right],$(11)

(12) $G:=\left[\begin{array}{c}{G}_1\\ G_2\\ G_3\\ G_4\end{array}\right],P:=\left[\begin{array}{c}{P}_1\\ P_2\\ P_3\\ 0\end{array}\right],R:=\left[\begin{array}{cccc}{R}_{11}& R_{12}& R_{13}& 0\\ R_{21}& R_{22}& R_{23}& 0\\ R_{31}& R_{32}& R_{33}& 0\\ 0& 0& 0& 0\end{array}\right],$(12)

where $E_{11},E_{22}$ are positive definite, and the matrices $J_{41}$ and $J_{33}-R_{33}$ are invertible (if the blocks are nonempty). The Kronecker index $\nu$ of the uncontrolled system satisfies

$\nu =\left\{\begin{array}{c}0ifandonlyif{n}_1=n_4=0and{n}_3=0,\\ 1ifandonlyif{n}_1=n_4=0and{n}_3>0,\\ 2ifandonlyif{n}_1=n_4>0.\end{array}\right.$

As initially stated, it is beneficial to preserve the structural properties of the original pH-DAE model during MOR. Structure-preserving MOR methods, therefore, search for a reduced pH-DAE

(13) ${\mathrm{\Sigma }}_r:\left\{\begin{array}{c}{E}_r{\dot{x}}_r\left(t\right)=\left(J_r-R_r\right)x_r\left(t\right)+\left(G_r-P_r\right)u\left(t\right),\\ y_r\left(t\right)=(G_r+P_r{)}^{\mathrm{\top }}{x}_r\left(t\right)+\left(S_r+N_r\right)u\left(t\right),\end{array}\right.$(13)

with $x_r(t)\in {\mathbb{R}}^r$, $r\ll n$ that fulfils the pH structural constraints, i.e. the associated system matrix may be decomposed as in Definition 2.2 such that

(14) ${\mathcal{P}}_r(s)=s{\mathcal{E}}_r+{\mathrm{\Gamma }}_r+W_r,$(14)

with symmetric positive semi-definite ${\mathcal{E}}_r,W_r$ and skew-symmetric ${\mathrm{\Gamma }}_r$. Note that the system matrix has also recently been used to derive a symplectic MOR method for LTI pH-ODEs without feedthrough in [20].

3. Our approach

In the following, we demonstrate how the concepts of the presented staircase form and the system matrix may be unified to derive a framework for tangential interpolation of pH-DAEs with an arbitrary Kronecker index. For this, we proceed in three steps. First, we apply a transformation under s.s.e. to the original system matrix $\mathcal{P}$ that enables us to decompose the original transfer function into proper parts and improper parts that may originate from algebraic constraints. Since all improper parts have to be preserved in the ROM exactly to keep the error $H-H_r$ bounded (see Remark 1), we propose a new method to efficiently reduce only the proper part in the second step. Third, we show how to reattach the original improper part to the reduced proper transfer function to construct a minimal pH-DAE representation of the ROM in staircase form.

3.1. System matrix decomposition

Let $\mathcal{P}$ denote the system matrix of a full-order pH-DAE system with transfer function $H$ in staircase form as in Lemma 2.3. For the sake of notational simplicity, we use $A=J-R$, $B=G-P$, $C=(G+P{)}^{\mathrm{\top }}$ and $D=S+N$ which are partitioned as in Lemma 2.3: $A_{11}\in {\mathbb{R}}^{n_1\times n_1}$, for example, denotes the upper left block of $J-R$. We define the transformation matrices ${\mathcal{T}}_1,{\mathcal{T}}_2\in {\mathbb{R}}^{(n+m)\times (n+m)}$ with

${\mathcal{T}}_1:=\left[\begin{array}{cc}L& 0\\ X& I_m\end{array}\right]\text{break}=\left[\begin{array}{ccccc}{I}_{n_1}& 0& -A_{13}{A}_{33}^{-1}& 0& 0\\ 0& I_{n_2}& -A_{23}{A}_{33}^{-1}& -(A_{21}-A_{23}{A}_{33}^{-1}{A}_{31})A_{41}^{-1}& 0\\ 0& 0& A_{33}^{-1}& 0& 0\\ 0& 0& 0& A_{41}^{-1}& 0\\ C_4{A}_{14}^{-1}& 0& (C_3-C_4{A}_{14}^{-1}{A}_{13})A_{33}^{-1}& 0& I_m\end{array}\right],$

${\mathcal{T}}_2:=\left[\begin{array}{cc}M& Y\\ 0& I_m\end{array}\right]\text{break}=\left[\begin{array}{ccccc}{I}_{n_1}& 0& 0& 0& -A_{41}^{-1}{B}_4\\ 0& I_{n_2}& 0& 0& 0\\ -A_{33}^{-1}{A}_{31}& -A_{33}^{-1}{A}_{32}& I_{n_3}& 0& -A_{33}^{-1}(B_3-A_{31}{A}_{41}^{-1}{B}_4)\\ 0& A_{14}^{-1}(-A_{12}+A_{13}{A}_{33}^{-1}{A}_{32})& 0& A_{14}^{-1}& 0\\ 0& 0& 0& 0& I_m\end{array}\right].$

It is apparent that ${\mathcal{T}}_1$ and ${\mathcal{T}}_2$ satisfy the conditions in Lemma 2.1 since the determinants of $A_{33},A_{14}$ and $A_{41}$ are constant and nonzero (see Lemma 2.3 and [13,14] for a proof). Note that our definitions of $L$ and $M$ share similarities with the state-space transformation presented in [14, Lemma 6] for autonomous semi-dissipative Hamiltonian DAEs. We obtain a transformed s.s.e. system matrix

(15) $˜\left(s\right)={\mathcal{T}}_1\left(s\right){\mathcal{T}}_2=\left[\begin{array}{ccccc}s{\tilde{E}}_{11}-{\tilde{A}}_{11}& 0& 0& -I_{n_1}& -{\tilde{B}}_1\\ 0& s{E}_p-A_p& 0& 0& -B_p\\ 0& 0& -I_{n_3}& 0& 0\\ -I_{n_1}& 0& 0& 0& 0\\ {\tilde{C}}_1& C_p& 0& 0& D_p+sD\infty \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\end{array}\right],$(15)

with nonsingular ${\tilde{E}}_{11},E_{\mathrm{p}}$. The second and fifth block column and row entries, respectively, are highlighted since these are the only parts that contribute to the transformed transfer function $\tilde{H}$:

(16) $\tilde{H}=\tilde{C}(s\tilde{E}-\tilde{A}{)}^{-1}\tilde{B}+D_{\mathrm{p}}+s{D}_{\mathrm{\infty }}$(16)

(17) $=\left[\begin{array}{cccc}{\tilde{C}}_1& C_{\mathrm{p}}& 0& 0\end{array}\right]\left[\begin{array}{cccc}0& 0& 0& -I_{n_1}\\ 0& {(s{E}_{\mathrm{p}}-A_{\mathrm{p}})}^{-1}& 0& 0\\ 0& 0& -I_{n_3}& 0\\ -I_{n_1}& 0& 0& -(s{\tilde{E}}_{11}-{\tilde{A}}_{11})\end{array}\right]\left[\begin{array}{c}{\tilde{B}}_1\\ B_{\mathrm{p}}\\ 0\\ 0\end{array}\right]+D_{\mathrm{p}}+s{D}_{\mathrm{\infty }}$(17)

(18) $=C_{\mathrm{p}}(s{E}_{\mathrm{p}}-A_{\mathrm{p}}{)}^{-1}{B}_{\mathrm{p}}+D_{\mathrm{p}}+s{D}_{\mathrm{\infty }}$(18)

Consequently, the transfer function of every pH-DAE may be represented as the sum of the transfer function of a proper ODE system with dimension $n_2$ and an additional linear, improper term. Moreover, as shown in [21, Theorem 1], the proper subsystem again satisfies the pH structural constraints in Definition 2.2. It, therefore, admits a pH-ODE representation, which can be found by using its system matrix.

Lemma 3.1.

A pH-ODE representation of the proper subsystem with transfer function $H_{\mathrm{p}}=C_{\mathrm{p}}(s{E}_{\mathrm{p}}-A_{\mathrm{p}}{)}^{-1}{B}_{\mathrm{p}}+D_{\mathrm{p}}$ can be computed by simply decomposing its system matrix

${\mathcal{P}}_{\mathrm{p}}(s):=\left[\begin{array}{cc}s{E}_{\mathrm{p}}-A_{\mathrm{p}}& -B_{\mathrm{p}}\\ C_{\mathrm{p}}& D_{\mathrm{p}}\end{array}\right],$

into symmetric and skew-symmetric parts, respectively.

Proof.

Decomposing ${\mathcal{P}}_{\mathrm{p}}$ yields

${\mathcal{P}}_{\mathrm{p}}(s)=s\left[\begin{array}{cc}{E}_{\mathrm{p}}& 0\\ 0& 0\end{array}\right]+{\mathrm{\Gamma }}_{\mathrm{p}}+W_{\mathrm{p}},$

with

${\mathrm{\Gamma }}_{\mathrm{p}}=\left[\begin{array}{cc}-J_{\mathrm{p}}& -G_{\mathrm{p}}\\ G_{\mathrm{p}}^{\mathrm{\top }}& N_{\mathrm{p}}\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}-A_{\mathrm{p}}+A_{\mathrm{p}}^{\mathrm{\top }}& -B_{\mathrm{p}}-C_{\mathrm{p}}^{\mathrm{\top }}\\ C_{\mathrm{p}}+B_{\mathrm{p}}^{\mathrm{\top }}& D_{\mathrm{p}}-D_{\mathrm{p}}^{\mathrm{\top }}\end{array}\right],$

$W_{\mathrm{p}}=\left[\begin{array}{cc}{R}_{\mathrm{p}}& P_{\mathrm{p}}\\ P_{\mathrm{p}}^{\mathrm{\top }}& S_{\mathrm{p}}\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}-A_{\mathrm{p}}-A_{\mathrm{p}}^{\mathrm{\top }}& -B_{\mathrm{p}}+C_{\mathrm{p}}^{\mathrm{\top }}\\ C_{\mathrm{p}}-B_{\mathrm{p}}^{\mathrm{\top }}& D_{\mathrm{p}}+D_{\mathrm{p}}^{\mathrm{\top }}\end{array}\right].$

The fact that the system is an ODE system follows directly from $E_{\mathrm{p}}=E_{22}>0$. This also proves condition (i) in Definition 2.2. In [21, Theorem 1], it was shown that the sum ${\mathrm{\Gamma }}_{\mathrm{p}}+W_{\mathrm{p}}$ may be obtained from a series of transformations of the original $\mathrm{\Gamma }+W$. These include permutations, Schur complement constructions and congruence-like transformations, which all preserve the positive semi-definiteness of the symmetric part. Therefore, ${\mathrm{\Gamma }}_{\mathrm{p}},W_{\mathrm{p}}$ fulfil the pH structural constraints (ii) and (iii) in Definition 2.2, which completes the proof.

We highlight that the simplicity of this result is a direct consequence of the staircase form and the pH structural constraints. For general DAE systems, this system decomposition approach generally requires the computation of spectral projectors onto the left and right deflating subspaces of the pencil $\lambda E-A$ corresponding to the finite eigenvalues, which are numerically challenging to compute in the large-scale setting [4]. In some applications, such as fluid flow problems or electric circuit simulation where the matrices $E$ and $A$ have a special block structure, the computation of spectral projectors can be done more efficiently or even circumvented, see, e.g [4,23]. However, the proposed interpolatory MOR approaches for general (unstructured) DAE systems [4,6] vary for different Kronecker indices, and their adaptations to pH-DAE systems proposed in [12,13] do not always guarantee that the ROM is again in pH form. In the following, we show how the results in this section enable the construction of a general MOR approach that works irrespective of the original system’s Kronecker index and guarantees to produce minimal ROMs in pH form.

3.2. Tangential interpolation

As discussed in Remark 1, the ROM has to match the improper part $s{D}_{\mathrm{\infty }}$ exactly because otherwise, $H-H_r\notin {\mathcal{R}\mathcal{H}}_{\mathrm{\infty }}^{m\times m}$. We will therefore set this part aside for now and focus on the reduction of the proper part. Since the proper subsystem has only dimension $n_2$, a natural approach would be to reduce the proper system matrix ${\mathcal{P}}_{\mathrm{p}}$ directly. According to (8), this approach requires the computation of solutions $v_i\in {\mathbb{C}}^{n_2}$ of

(19) $({\sigma }_i{E}_{\mathrm{p}}-A_{\mathrm{p}})v_i=B_{\mathrm{p}}{b}_i,$(19)

for all $i=1,\dots r$. From (15), we derive

$E_{\mathrm{p}}=E_{22},A_{\mathrm{p}}=A_{22}-A_{23}{A}_{33}^{-1}{A}_{32}.$

For index-1 pH-DAEs ($n_3>0$) with nonzero $A_{23}{A}_{33}^{-1}{A}_{32}$, the matrix $A_{\mathrm{p}}$ might be dense, and therefore, the solutions in (19) may be more expensive than for the original system with matrices $E$, $A$ and $B$. This is illustrated by the following example.

Figure 1. Sparsity patterns of exemplary full-order matrices $E,A\in {\mathbb{R}}^{n\times n}$ and of the matrices $E_{\mathrm{p}},A_{\mathrm{p}}\in {\mathbb{R}}^{n_2\times n_2}$ of its associated proper subsystem.

Example 3.2.

Consider a pH-DAE with Kronecker index $\nu =1$ in staircase form. The system has dimension $n={10}^4$ with $n_2=n_3=\frac{n}{2}$ and one input-output pair. As shown by the sparsity patterns in Figure 1(a,b), the matrices $E$ and $A$ have a very simple structure with only few nonzero entries (depicted in blue). However, the matrix $A_{\mathrm{p}}$ of the proper subsystem is dense, as shown in Figure 1(d).

To demonstrate the effects on interpolatory MOR, we solved (19) for the proper subsystem and for the original system with matrices $E$, $A$ and $B$ using $100$ random complex shifts ${\sigma }_i$ with MATLAB’s mldivide command. All computations were conducted using MATLAB R2021b (version 9.11.0.1873467) on an Intel® Core^™ i7–8700 CPU (3.20 GHz, 6-Core) with 32 GB RAM. The computations of (19) for the proper subsystem took $2.51$ seconds on average versus only $0.013$ seconds for the original model. Consequently, even though the proper subsystem is significantly smaller in size, the computation of the reduction matrix $V$ takes more than 150 times longer.

Therefore, we propose another approach that works with the original (sparse) system matrix, irrespective of the system’s Kronecker index.

Theorem 3.3.

Given a large-scale pH-DAE in staircase form with system matrix $\mathcal{P}$, as well as interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$, and corresponding right tangential directions $\{b_1,\dots ,b_r\}$, let $V\in {\mathbb{R}}^{n\times r}$ define a reduction matrix such that (8) holds with a decomposition $V^{\mathrm{\top }}=[V_1^{\mathrm{\top }},V_2^{\mathrm{\top }},V_3^{\mathrm{\top }},V_4^{\mathrm{\top }}{]}^{\mathrm{\top }}$ with $V_j\in {\mathbb{R}}^{n_j\times r}$ for all $j=1,\dots ,4$. We define the reduction matrices

(20) $\mathcal{U}:=\left[\begin{array}{cc}0& A_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }}\\ V_2& 0\\ -A_{33}^{-\mathrm{\top }}{A}_{23}^{\mathrm{\top }}{V}_2& A_{33}^{-\mathrm{\top }}\left(C_3^{\mathrm{\top }}-A_{13}^{\mathrm{\top }}{A}_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }}\right)\\ 0& 0\\ 0& I_m\end{array}\right],\mathcal{V}:=\left[\begin{array}{cc}0& A_{14}^{-\mathrm{\top }}{B}_4\\ V_2& 0\\ 0& 0\\ 0& 0\\ 0& I_m\end{array}\right],$(20)

Then, the model associated with the reduced system matrix

${\mathcal{P}}_r:={\mathcal{U}}^{\mathrm{\top }}\mathcal{P}\mathcal{V},$

satisfies the tangential interpolation conditions (7), and its proper subsystem fulfils the pH structural constraints. The reduced system matrix admits a decomposition

${\mathcal{P}}_r=s\left[\begin{array}{cc}{E}_{\mathrm{p},r}& 0\\ 0& D_{\mathrm{\infty }}\end{array}\right]+\left[\begin{array}{cc}-J_{\mathrm{p},r}& -G_{\mathrm{p},r}\\ G_{\mathrm{p},r}^{\mathrm{\top }}& N_{\mathrm{p},r}\end{array}\right]+\left[\begin{array}{cc}{R}_{\mathrm{p},r}& P_{\mathrm{p},r}\\ P_{\mathrm{p},r}^{\mathrm{\top }}& S_{\mathrm{p},r}\end{array}\right].$

Proof.

Using the decomposition of $V$, the original system matrix $\mathcal{P}$ may initially be reduced in the following way:

(21) ${\mathcal{P}}_2(s)=\left[\begin{array}{cc}{U}^{\mathrm{\top }}& 0\\ 0& I_m\end{array}\right]\mathcal{P}\left[\begin{array}{cc}U& 0\\ 0& I_m\end{array}\right],U:=\left[\begin{array}{cccc}{I}_{n_1}& 0& 0& 0\\ 0& V_2& 0& 0\\ 0& 0& I_{n_3}& 0\\ 0& 0& 0& I_{n_4}\end{array}\right].$(21)

The reduced model associated with ${\mathcal{P}}_2$ satisfies the tangential interpolation conditions in (7) since $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(V)\subseteq \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(U)$. Note that this approach, which was also proposed in [13] for index-1 pH-DAEs, produces ROMs that are still comparatively large: the ROM has dimension $n_1+r+n_3+n_4$. However, since the ROM is still a pH-DAE in staircase form, we can obtain a minimal representation by applying a transformation under s.s.e., as in Section 3.1, and extracting the proper system matrix without changing its transfer function. A combination of the reduction in (21), the transformation and extraction of proper parts yields the reduction matrices $\mathcal{U},\mathcal{V}$. Moreover, simple algebraic manipulations show that

${\mathcal{P}}_r={\mathcal{U}}^{\mathrm{\top }}\mathcal{P}\mathcal{V}=\left[\begin{array}{cc}{V}_2^{\mathrm{\top }}& 0\\ 0& I_m\end{array}\right]\left[\begin{array}{cc}s{E}_{\mathrm{p}}-A_{\mathrm{p}}& -B_{\mathrm{p}}\\ C_{\mathrm{p}}& D_{\mathrm{p}}+s{D}_{\mathrm{\infty }}\end{array}\right]\left[\begin{array}{cc}{V}_2& 0\\ 0& I_m\end{array}\right],$

and consequently, the proper part of ${\mathcal{P}}_r$ fulfils the pH structural constraints and may be decomposed to obtain a pH representation, which completes the proof.

To improve the numerical stability of Krylov subspace methods, the matrix $V$ is usually orthogonalized such that $V^{\mathrm{\top }}V=I_r$. This orthogonalization does not change the moment matching conditions in (8) since these only depend on the subspace that is spanned by the column vectors of $V$, not the basis itself. However, note that in our case, even if $V$ is orthogonal, this is generally not the case for its submatrix $V_2$. To improve the numerical stability, we employ the cosine-sine decomposition, as discussed in [24]. For this, we split $V$ into two parts: $V_2$ and the remaining submatrices $V_{\mathrm{r}\mathrm{e}\mathrm{m}}^{\mathrm{T}}=[V_1^{\mathrm{\top }},V_3^{\mathrm{\top }},V_4^{\mathrm{\top }}{]}^{\mathrm{\top }}\in {\mathbb{R}}^{n_{\mathrm{r}\mathrm{e}\mathrm{m}}\times r}$. We then compute the decomposition

(22) $\left[\begin{array}{c}{V}_2\\ V_{\mathrm{r}\mathrm{e}\mathrm{m}}\end{array}\right]=\left[\begin{array}{cc}{\bar{V}}_2& 0\\ 0& {\bar{V}}_{\mathrm{r}\mathrm{e}\mathrm{m}}\end{array}\right]\left[\begin{array}{c}{C}_s\\ S_s\end{array}\right]X_s^{\mathrm{\top }},$(22)

with orthogonal ${\bar{V}}_2\in {\mathbb{R}}^{n_2\times r}$, ${\bar{V}}_{\mathrm{r}\mathrm{e}\mathrm{m}}\in {\mathbb{R}}^{n_{\mathrm{r}\mathrm{e}\mathrm{m}}\times r}$, and $X_s\in {\mathbb{R}}^{r\times r}$, as well as diagonal $C_s,S_s\in {\mathbb{R}}^{r\times r}$ such that $C_s^{\mathrm{\top }}{C}_s+S_s^{\mathrm{\top }}{S}_s=I_r$. Replacing $V_2$ in (20) with ${\bar{V}}_2$ yields the final reduced system matrix ${\mathcal{P}}_r$.

3.3. Minimal pH-DAE representation

To find a pH-DAE representation for ${\mathcal{P}}_r$, we have to incorporate the improper part $s{D}_{\mathrm{\infty }}$ that has been separated back into the model. Two different methods have been proposed in [21] and [25] for this purpose. Since the method in [25] only leads to a minimal ROM representation if $D_{\mathrm{\infty }}$ has full rank, we proceed similarly as in [21]. We have that

(23) $D_{\mathrm{\infty }}=-C_4{A}_{14}^{-1}{E}_{11}{A}_{41}^{-1}{B}_4=G_4^{\mathrm{\top }}{A}_{41}^{-\mathrm{\top }}{E}_{11}{A}_{41}^{-1}{G}_4=D_{\mathrm{\infty }}^{\mathrm{\top }}\ge 0.$(23)

Consequently, there exists a rank-revealing factorization $D_{\mathrm{\infty }}=L_{\mathrm{\infty }}{L}_{\mathrm{\infty }}^{\mathrm{\top }}$ with $L_{\mathrm{\infty }}\in {\mathbb{R}}^{m\times q}$. A minimal ROM representation ${\mathrm{\Sigma }}_r$ in staircase form can be found with

$E_r=\left[\begin{array}{ccc}{I}_q& 0& 0\\ 0& E_{\mathrm{p},r}& 0\\ 0& 0& 0\end{array}\right],J_r=\left[\begin{array}{ccc}0& 0& I_q\\ 0& J_{\mathrm{p},r}& 0\\ -I_q& 0& 0\end{array}\right],R_r=\left[\begin{array}{ccc}0& 0& 0\\ 0& R_{\mathrm{p},r}& 0\\ 0& 0& 0\end{array}\right],$

$G_r=\left[\begin{array}{c}0\\ G_{\mathrm{p},r}\\ L_{\mathrm{\infty }}^{\mathrm{\top }}\end{array}\right],P_r=\left[\begin{array}{c}0\\ P_{\mathrm{p},r}\\ 0\end{array}\right],S_r=S_{\mathrm{p},r},N_r=N_{\mathrm{p},r}.$

We combine the obtained results of this section in Algorithm 1.

Algorithm 1:

Tangential Interpolation of pH-DAEs

Input : Large-scale pH-DAE $\mathrm{\Sigma }$ in staircase form with system matrix $\mathcal{P}$; set of interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$ and corresponding right tangential directions $\{b_1,\dots ,b_r\}$ (both closed under complex conjugation).

Output: Reduced pH-DAE ${\mathrm{\Sigma }}_r$ with system matrix ${\mathcal{P}}_r$.

1 Compute $V={\left[V_1^{\mathrm{\top }},V_2^{\mathrm{\top }},V_3^{\mathrm{\top }},V_4^{\mathrm{\top }}\right]}^{\mathrm{\top }}\in {\mathbb{R}}^{n\times r}$ such that (8) holds.

2 Orthogonalize $V_2$ via the cosine-sine decomposition in (22) such that $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(V_2)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}({\bar{V}}_2)$ with ${\bar{V}}_2^{\mathrm{\top }}{\bar{V}}_2=I_r$.

3 Compute the reduction matrices $\mathcal{U},\mathcal{V}$ as in Section 3.2:

$\mathcal{U}=\left[\begin{array}{cc}0& A_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }}\\ {\bar{V}}_2& 0\\ -A_{33}^{-\mathrm{\top }}{A}_{23}^{\mathrm{\top }}{\bar{V}}_2& A_{33}^{-\mathrm{\top }}(C_3^{\mathrm{\top }}-A_{13}^{\mathrm{\top }}{A}_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }})\\ 0& 0\\ 0& I_m\end{array}\right],\mathcal{V}=\left[\begin{array}{cc}0& A_{14}^{-\mathrm{\top }}{B}_4\\ {\bar{V}}_2& 0\\ 0& 0\\ 0& 0\\ 0& I_m\end{array}\right].$

4 Compute and decompose the reduced system matrix

${\mathcal{P}}_r={\mathcal{U}}^{\mathrm{\top }}\mathcal{P}\mathcal{V}=s\left[\begin{array}{cc}{E}_{\mathrm{p},r}& 0\\ 0& D_{\mathrm{\infty }}\end{array}\right]+\left[\begin{array}{cc}-J_{\mathrm{p},r}& -G_{\mathrm{p},r}\\ G_{\mathrm{p},r}^{\mathrm{\top }}& N_{\mathrm{p},r}\end{array}\right]+\left[\begin{array}{cc}{R}_{\mathrm{p},r}& P_{\mathrm{p},r}\\ P_{\mathrm{p},r}^{\mathrm{\top }}& S_{\mathrm{p},r}\end{array}\right].$

5 Compute a rank-revealing factorization $D_{\mathrm{\infty }}=L_{\mathrm{\infty }}{L}_{\mathrm{\infty }}^{\mathrm{\top }}$ with $L_{\mathrm{\infty }}\in {\mathbb{R}}^{m\times q}$.

6 if $q>0$ then

7 Construct ROM ${\mathrm{\Sigma }}_r$ as in (13) with

$E_r=\left[\begin{array}{ccc}{I}_q& 0& 0\\ 0& E_{\mathrm{p},r}& 0\\ 0& 0& 0\end{array}\right],J_r=\left[\begin{array}{ccc}0& 0& I_q\\ 0& J_{\mathrm{p},r}& 0\\ -I_q& 0& 0\end{array}\right],R_r=\left[\begin{array}{ccc}0& 0& 0\\ 0& R_{\mathrm{p},r}& 0\\ 0& 0& 0\end{array}\right],$

$G_r=\left[\begin{array}{c}0\\ G_{\mathrm{p},r}\\ L_{\mathrm{\infty }}^{\mathrm{\top }}\end{array}\right],P_r=\left[\begin{array}{c}0\\ P_{\mathrm{p},r}\\ 0\end{array}\right],S_r=S_{\mathrm{p},r},N_r=N_{\mathrm{p},r}.$

8 else

9 Construct ROM ${\mathrm{\Sigma }}_r$ as in (13) directly from ${\mathcal{P}}_r$.

10 end

3.4. Discussion

In the following, we would like to briefly discuss the differences between the proposed Rosenbrock framework and other interpolatory MOR approaches that have been proposed for pH-DAEs. On the one hand, the approaches get more complex if the system’s Kronecker index $\nu$ increases. On the other hand, as revealed in this section, if the system’s transfer function has improper parts, these require special care since they have to be preserved exactly in the ROM. As shown in Lemma 2.3, the Kronecker index of LTI pH-DAEs is at most two, and improper parts may only occur if the Kronecker index is two (see (23)). Therefore, in the context of MOR, we may identify six different system categories. An overview of existing methods for each category, to the best of the author’s knowledge, is given in Table 1.

Table 1. Proposed tangential interpolation methods for pH-DAEs. The six different categories are derived from the original system’s Kronecker index $\nu$ and whether its transfer function has improper parts.

Category	$\nu$	$n_2$	$n_3$	$n_4/n_1$	$D_{\mathrm{\infty }}$	References
Index-0	$0$	$\ne 0$	$0$	$0$	$0$	[22, Theorem 7]
Index-1	$1$	$\ne 0$	$\ne 0$	$0$	$0$	[13, Theorems 1 and 2]
						[12, Theorem 5]
Proper Index-2	$2$	$\ne 0$	$0$	$\ne 0$	$0$	[13, Theorem 3]
Improper Index-2	$2$	$\ne 0$	$0$	$\ne 0$	$\ne 0$	[13, Theorem 4]
Proper Index-1–2	$2$	$\ne 0$	$\ne 0$	$\ne 0$	$0$	n/a
Improper Index-1–2	$2$	$\ne 0$	$\ne 0$	$\ne 0$	$\ne 0$	n/a

The tangential interpolation of pH-ODEs ($\nu =0$), as initially proposed in [22], is well understood and leads to minimal ROM representations in pH-ODE form. For these systems, our approach is equivalent since the reduction matrices simplify to $\mathcal{U}=\mathcal{V}=\left[\begin{array}{cc}{\bar{V}}_2& 0\\ 0& I_m\end{array}\right]$. For index-1 pH-DAEs, methods in [12,13] have been proposed, which rely on different semi-explicit representations of the FOM. In [13, Theorem 1], the feedthrough matrix of the ROM is modified to match the polynomial part of the FOM, which in the index-1 case, is constant. However, this method does not guarantee that the ROM fulfils the pH structural constraints. A remedy to this problem is to preserve the algebraic constraints of the FOM as proposed in [13, Theorem 2] and [12, Theorem 5]. However, this does not generally yield minimal ROMs since redundant algebraic equations cannot be removed, as discussed in [13, Remark 2]. In contrast, our method does not impose additional assumptions since it also only requires a semi-explicit FOM representation but guarantees the preservation of the pH form and yields minimal ROM representations.

For systems with $\nu =2$ that do not have an index-1 part ($n_3=0$), we may distinguish between the proper and improper case. The proper index-2 case is comparable to the pH-ODE case if the system is in semi-explicit form and the method proposed in [13] yields minimal ROMs in pH form. Our method achieves the same goals, but may require one additional transformation to identify the full-rank matrix $J_{41}$. The improper index-2 case is treated in [13, Theorem 4]. However, this method does not guarantee that the ROM fulfils the pH structural constraints, which is generally the case for our method, again under the assumption that the system is in staircase form. To the best of the author’s knowledge, methods for systems with index-1 and index-2 parts ($n_3\ne 0$, $n_4\ne 0$) have not been proposed yet. These two categories are also covered by our framework, which applies to all pH-DAEs in staircase form.

4. ${\mathcal{H}}_{\mathbf{2}}$- and ${\mathcal{H}}_{\mathrm{\infty }}$-inspired tangential interpolation

The question that remains is how to choose the parameters in Algorithm 1, i.e. the set of interpolation points and tangential directions. For unstructured ODE and DAE systems, different ${\mathcal{H}}_2$- and ${\mathcal{H}}_{\mathrm{\infty }}$-inspired strategies have been proposed (see, e.g. [7, 28, 29, 31, 32]) which have also partly been adapted to the pH-ODE and pH-DAE cases in [13,22].

In the following, we demonstrate how these ideas may be incorporated into the proposed Rosenbrock framework, and refer to the cited work in each section for implementation details.

4.1. Interpolatory ${\mathcal{H}}_{\mathbf{2}}$ approximation

In ${\mathcal{H}}_2$-inspired interpolation methods for general DAE systems, the transfer function of the FOM is typically decomposed into the sum $H(s)=H_{\mathrm{s}\mathrm{p}}(s)+H_{\mathrm{p}\mathrm{o}\mathrm{l}}$ with a strictly proper part $H_{\mathrm{s}\mathrm{p}}$ satisfying ${lim}_{s\to \mathrm{\infty }}{H}_{\mathrm{s}\mathrm{p}}(s)=0$ and a polynomial part $H_{\mathrm{p}\mathrm{o}\mathrm{l}}$ that potentially grows polynomially for $s\to \mathrm{\infty }$. For the ${\mathcal{H}}_2$ error $\parallel H-H_r{\parallel }_{{\mathcal{H}}_2}$ to be bounded, we require $H-H_r\in {\mathcal{R}\mathcal{H}}_2^{m\times m}$. Assuming a similar decomposition of $H_r$ into $H_r(s)=H_{\mathrm{s}\mathrm{p},\mathrm{r}}(s)+H_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{r}}(s)$, this is only the case if $H_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{r}}=H_{\mathrm{p}\mathrm{o}\mathrm{l}}$. The necessary conditions for locally ${\mathcal{H}}_2$-optimal ROMs may then be formulated as interpolation conditions.

Lemma 4.1.

[6, Theorem 9.2.1] Let the transfer function $H(s)=H_{\mathrm{s}\mathrm{p}}(s)+H_{\mathrm{p}\mathrm{o}\mathrm{l}}$ be decomposed into strictly proper and polynomial parts. Let the ROM transfer function $H_r$ have an analogous decomposition $H_r(s)=H_{\mathrm{s}\mathrm{p},\mathrm{r}}(s)+H_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{r}}(s)$ with strictly proper part $H_{\mathrm{s}\mathrm{p},\mathrm{r}}(s)=C_{\mathrm{s}\mathrm{p},\mathrm{r}}(s{E}_{\mathrm{s}\mathrm{p},\mathrm{r}}-A_{\mathrm{s}\mathrm{p},\mathrm{r}}{)}^{-1}{B}_{\mathrm{s}\mathrm{p},\mathrm{r}}$ such that $A_{\mathrm{s}\mathrm{p},\mathrm{r}}\in {\mathbb{R}}^{r\times r}$, nonsingular $E_{\mathrm{s}\mathrm{p},\mathrm{r}}\in {\mathbb{R}}^{r\times r}$ and $B_{\mathrm{s}\mathrm{p},\mathrm{r}},C_{\mathrm{s}\mathrm{p},\mathrm{r}}^{\mathrm{\top }}\in {\mathbb{R}}^{r\times m}$. If $H_r$ minimizes the ${\mathcal{H}}_2$ error $\parallel H-H_r{\parallel }_{{\mathcal{H}}_2}$ over all ROMs with an $r$-th order strictly proper part, then $H_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{r}}=H_{\mathrm{p}\mathrm{o}\mathrm{l}}$ and $H_{\mathrm{s}\mathrm{p},\mathrm{r}}$ minimizes the ${\mathcal{H}}_2$ error $\parallel H_{\mathrm{s}\mathrm{p}}-H_{\mathrm{s}\mathrm{p},\mathrm{r}}{\parallel }_{{\mathcal{H}}_2}$. Let $H_{\mathrm{s}\mathrm{p},\mathrm{r}}$ be represented by its pole-residue expansion $H_{\mathrm{s}\mathrm{p},\mathrm{r}}(s)=\sum _{i=1}^r\frac{l_i{r}_i^{\mathrm{\top }}}{s-{\lambda }_i}$ where $r_i\in {\mathbb{C}}^m,l_i\in {\mathbb{C}}^m$ and with simple poles ${\lambda }_i\in \mathbb{C}$. Then, the tangential interpolation conditions

(24a) $H(-{\lambda }_i)r_i=H_r(-{\lambda }_i)r_i,$(24a)

(24b) $l_i^{\mathrm{\top }}H(-{\lambda }_i)=l_i^{\mathrm{\top }}{H}_r(-{\lambda }_i),$(24b)

(24c) $l_i^{\mathrm{\top }}{H}^{\prime }(-{\lambda }_i)r_i=l_i^{\mathrm{\top }}{H}_r^{\prime }(-{\lambda }_i)r_i,$(24c)

hold for all $i=1,\dots ,r$.

This connection between interpolatory and ${\mathcal{H}}_2$-optimal model reduction is the motivation behind the well-known iterative rational Krylov algorithm (IRKA) [7], which utilizes a fixed-point iteration to enforce the necessary ${\mathcal{H}}_2$ optimality conditions for general DAE systems. Since for pH systems, the matrix $U$ in (6) is typically chosen to enforce structure preservation, fewer degrees of freedom are available for interpolation, and it is generally only possible to fulfil a subset of (24), e.g. the conditions in (24a). Enforcing these conditions in an iterative manner using (8) leads to the IRKA-PH algorithm proposed for pH-ODE systems in [22].

Embedding the IRKA-PH algorithm into the proposed pH-DAE framework is straightforward. In each IRKA-PH iteration, we first compute the reduced system matrix ${\mathcal{P}}_r$, using some initial interpolation data in the first iteration. We directly obtain the reduced strictly proper transfer function

$H_{\mathrm{s}\mathrm{p},\mathrm{r}}(s)=(G_{\mathrm{p},r}+P_{\mathrm{p},r}{)}^{\mathrm{\top }}(s{E}_{\mathrm{p},r}-(J_{\mathrm{p},r}-R_{\mathrm{p},r}){)}^{-1}(G_{\mathrm{p},r}-P_{\mathrm{p},r}).$

Suppose that the pencil $\lambda E_{\mathrm{p},r}-(J_{\mathrm{p},r}-R_{\mathrm{p},r})$ has simple eigenvalues ${\lambda }_i$ and let $t_i\in {\mathbb{C}}^{\mathrm{r}}$ denote a left eigenvector associated with ${\lambda }_i$, i.e.

(25) $t_i^{\mathrm{\top }}({\lambda }_i{E}_{\mathrm{p},r}-(J_{\mathrm{p},r}-R_{\mathrm{p},r}))=0.$(25)

The (right) residual direction is then given by $r_i=(G_{\mathrm{p},r}-P_{\mathrm{p},r}{)}^{\mathrm{\top }}{t}_i$ and to enforce the interpolation conditions in (24a), we set ${\sigma }_i=-{\lambda }_i$ and $b_i=r_i$ for all $i=1,\dots ,r$. This procedure is repeated, and upon convergence, the ROM satisfies the subset (24a) of ${\mathcal{H}}_2$ optimality conditions. Afterwards, we may attach the polynomial part $H_{\mathrm{p}\mathrm{o}\mathrm{l}}(s)=D_{\mathrm{p}}+s{D}_{\mathrm{\infty }}$ to the strictly proper ROM as described in Section 3.3. One disadvantage of this approach, which is summarized in Algorithm 2, is that a new ROM is computed in each iteration. In [27], an adaptation named CIRKA-PH was proposed, which has the potential to significantly accelerate IRKA-PH, especially in large-scale settings for which interpolatory methods are particularly powerful. Embedding CIRKA-PH into the pH-DAE framework works in a similar way.

Algorithm 2:

IRKA-PH for pH-DAEs (based on [22])

Input :Large-scale pH-DAE $\mathrm{\Sigma }$ in staircase form; set of interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$ and corresponding right tangential directions $\{b_1,\dots ,b_r\}$ (both closed under complex conjugation).

Output: Reduced pH-DAE ${\mathrm{\Sigma }}_r$.

• 1 while not converged do

• 2 Perform steps 1-4 of Algorithm 1.

• 3 Compute $t_i\in {\mathbb{C}}^r,{\lambda }_i\in \mathbb{C}$ solving (25) for all $i=1,\dots ,r$.

• 4 ${\sigma }_i\leftarrow -{\lambda }_i$ and $b_i\leftarrow (G_{\mathrm{p},r}-P_{\mathrm{p},r}{)}^{\mathrm{\top }}{t}_i\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}i=1,\dots ,r$

• 5 end

• 6 Perform steps 5-10 of Algorithm 1.

4.2. Adaptive interpolation

In practice, besides the computational expense of IRKA-PH, it may sometimes be difficult to determine a suitable reduced order $r$ in the first place. In [28,29], an adaptive approach was proposed that tackles this problem by iteratively adding new interpolation data in a complex region $\mathcal{S}$ where the approximation quality of the ROM is still poor. For pH-DAEs, this approach initially requires the computation of the proper system matrix ${\mathcal{P}}_{\mathrm{p}}$, as described in Section 3.1. The approximation quality of a ROM generated in steps 1–4 of Algorithm 1 can then be assessed at points $\mu \in \mathcal{S}$ with the following residual matrix [29]:

(26) $\zeta (\mu ):=(A_{\mathrm{p}}-\mu E_{\mathrm{p}}){\bar{V}}_2(J_{\mathrm{p},r}-R_{\mathrm{p},r}-\mu E_{\mathrm{p},r}{)}^{-1}(G_{\mathrm{p},r}-P_{\mathrm{p},r})-B_{\mathrm{p}}.$(26)

In each iteration, a new interpolation point ${\sigma }_{r+1}$ is added at the point in $\mathcal{S}$ where the norm of this residual matrix reaches its maximum, and a similar approach is taken to compute a new corresponding tangential direction $b_{r+1}$. Since ${\sigma }_{r+1}$ and $b_{r+1}$ are generally complex, their complex conjugates ${\bar{\sigma }}_{r+1}$ and ${\bar{b}}_{r+1}$ are also added to the interpolation data to keep it closed under conjugation. This way, the ROM dimension $r$ increases in each iteration until its transfer function does not significantly change between two subsequent iterations or the predefined maximum dimension $r_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is reached. For strategies on how to choose $b_{r+1}$ and update the complex region $\mathcal{S}$, the interested reader is referred to [29]. The general approach is summarized in Algorithm 3.

Algorithm 3:

TRKSM-PH for pH-DAEs (based on [28,29])

Input: Large-scale pH-DAE $\mathrm{\Sigma }$ in staircase form; set of interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$ and corresponding right tangential directions $\{b_1,\dots ,b_r\}$ (both closed under complex conjugation); maximum reduced order $r_{\mathrm{m}\mathrm{a}\mathrm{x}}>r_0$; initial complex region $\mathcal{S}$.

Output:Reduced pH-DAE ${\mathrm{\Sigma }}_r$.

• 1 Compute and decompose ${\mathcal{P}}_{\mathrm{p}}$ as in Section 3.1.

• 2 while not converged and $r<r_{\mathrm{m}\mathrm{a}\mathrm{x}}$ do

• 3 Perform steps 1-4 of Algorithm 1.

• 4 Solve ${\sigma }_{r+1}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}{}_{\mu \in \mathcal{S}}\parallel \zeta (\mu )\parallel$.

• 5 Solve $b_{r+1}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}{}_{\parallel d\parallel =1}\parallel \zeta ({\sigma }_{r+1})d\parallel$.

• 6 Add $({\sigma }_{r+1},{\bar{\sigma }}_{r+1})$ and $(b_{r+1},{\bar{b}}_{r+1})$ to the interpolation data.

• 7 Update the complex region $\mathcal{S}$.

• 8 end

• 9 Perform steps 5-10 of Algorithm 1.

4.3. Interpolatory ${\mathcal{H}}_{\mathrm{\infty }}$ approximation

So far, we have enforced that $S_{\mathrm{p},r}=S_{\mathrm{p}}$ and $N_{\mathrm{p},r}=N_{\mathrm{p}}$ to keep the ${\mathcal{H}}_2$ error bounded. For ${\mathcal{H}}_{\mathrm{\infty }}$-inspired MOR, we only require $H-H_r\in {\mathcal{R}\mathcal{H}}_{\mathrm{\infty }}^{m\times m}$ and thus, the reduced feedthrough matrices pose additional degrees of freedom that may be exploited in a similar manner as proposed in [30–32] for unstructured ODE systems. We may add structure-preserving perturbations to the feedthrough matrices $N_{\mathrm{p},r}$ and $S_{\mathrm{p},r}$ in the following way:

(27) ${\hat{N}}_{\mathrm{p},r}=N_{\mathrm{p},r}+{\mathrm{\Delta }}_N,{\mathrm{\Delta }}_N=-{\mathrm{\Delta }}_N^{\mathrm{\top }},$(27)

(28) ${\hat{S}}_{\mathrm{p},r}=S_{\mathrm{p},r}+{\mathrm{\Delta }}_S,{\mathrm{\Delta }}_S={\mathrm{\Delta }}_S^{\mathrm{\top }}\ge 0.$(28)

Simply adding these perturbations would only change the direct feedthrough of the ROM but not the dynamics and is therefore not expected to improve the ${\mathcal{H}}_{\mathrm{\infty }}$ approximation quality significantly. However, as shown in [30], it is also possible to perturb the feedthrough matrix of the ROM while retaining predefined interpolation conditions with the FOM – and the same holds for pH-DAEs.

Lemma 4.2.

Assume that we have obtained a reduced system matrix ${\mathcal{P}}_r=s{\mathcal{E}}_r+{\mathrm{\Gamma }}_r+W_r$ in steps 1–4 of Algorithm 1 using a set of interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$ and corresponding right tangential directions $\{b_1,\dots ,b_r\}$, both closed under complex conjugation. Let $F\in {\mathbb{R}}^{n_2\times m}$ be a solution to

$F^{\mathrm{\top }}{\bar{V}}_2=[b_1,\dots ,b_r]T_v,$

with $T_v\in {\mathbb{C}}^{r\times r}$ such that ${\bar{V}}_2=V_2{T}_v$. If the system matrix ${\mathcal{P}}_r$ is perturbed such that

${\hat{\mathcal{P}}}_r=s{\mathcal{E}}_r+{\hat{\mathrm{\Gamma }}}_r+{\hat{W}}_r,$

with

(29) ${\hat{\mathrm{\Gamma }}}_r=\left[\begin{array}{cc}-{\hat{J}}_{\mathrm{p},r}& -{\hat{G}}_{\mathrm{p},r}\\ {\hat{G}}_{\mathrm{p},r}^{\mathrm{T}}& {\hat{N}}_{\mathrm{p},r}\end{array}\right]={\mathrm{\Gamma }}_r+\left[\begin{array}{c}{\bar{V}}_2^{\mathrm{\top }}F\\ -I_m\end{array}\right]{\mathrm{\Delta }}_N{\left[\begin{array}{c}{\bar{V}}_2^{\mathrm{\top }}F\\ -I_m\end{array}\right]}^{\mathrm{\top }},$(29)

(30) ${\hat{W}}_r=\left[\begin{array}{cc}{\hat{R}}_{\mathrm{p},r}& {\hat{P}}_{\mathrm{p},r}\\ {\hat{P}}_{\mathrm{p},r}^{\mathrm{\top }}& {\hat{S}}_{\mathrm{p},r}\end{array}\right]=W_r+\left[\begin{array}{c}{\bar{V}}_2^{\mathrm{\top }}F\\ -I_m\end{array}\right]{\mathrm{\Delta }}_S{\left[\begin{array}{c}{\bar{V}}_2^{\mathrm{\top }}F\\ -I_m\end{array}\right]}^{\mathrm{\top }},$(30)

then the perturbed ROM ${\hat{\mathrm{\Sigma }}}_r$ with transfer function ${\hat{H}}_r$ obtained by steps 5–10 in Algorithm 1 is a pH-DAE system, and it holds that

(31) $H({\sigma }_i)b_i=H_r({\sigma }_i)b_i={\hat{H}}_r({\sigma }_i)b_i,$(31)

for all $i=1,\dots ,r$ and for any ${\mathrm{\Delta }}_N=-{\mathrm{\Delta }}_N^{\mathrm{\top }}$ and ${\mathrm{\Delta }}_S={\mathrm{\Delta }}_S^{\mathrm{\top }}\ge 0$.

Proof.

The fact that the perturbed system ${\hat{\mathrm{\Sigma }}}_r$ fulfils the pH structural constraints follows directly from the properties of ${\mathrm{\Delta }}_N$ and ${\mathrm{\Delta }}_S$. The proof of (31) follows the proof of Theorem 3 in [30] for general LTI systems and is therefore omitted here.

This result enables us to optimize the new degrees of freedom ${\mathrm{\Delta }}_N,{\mathrm{\Delta }}_S$ in an ${\mathcal{H}}_{\mathrm{\infty }}$-inspired way. Optimizing the ${\mathcal{H}}_{\mathrm{\infty }}$ error directly is challenging since ${\mathcal{H}}_{\mathrm{\infty }}$ norm computations are computationally taxing, and the ${\mathcal{H}}_{\mathrm{\infty }}$ norm depends non-smoothly on ${\mathrm{\Delta }}_N,{\mathrm{\Delta }}_S$ (see [33] Sect 3.2.1]). Instead, the SOBMOR algorithm, proposed in [33], may be employed. Therein, the functions $vtu\left(\cdot \right)$ and $vtu\left(\cdot \right)$ are introduced, which map vectors row-wise to appropriately sized upper triangular and strictly upper triangular matrices, respectively. The function names are abbreviations for vector-to-upper and vector-to-strictly-upper, respectively. Using these functions, we may define parameter vectors ${\theta }_N\in {\mathbb{R}}^{m(m-1)/2}$ and ${\theta }_S\in {\mathbb{R}}^{m(m+1)/2}$ and design ${\mathrm{\Delta }}_N,{\mathrm{\Delta }}_S$ in the following way:

${\mathrm{\Delta }}_N\left({\theta }_N\right):=vtsu{\left({\theta }_N\right)}^{\mathrm{\top }}-vtsu\left({\theta }_N\right),$

${\mathrm{\Delta }}_S\left({\theta }_S\right):=vtsu{\left({\theta }_S\right)}^{\mathrm{\top }}-vtsu\left({\theta }_S\right).$

Finally, a levelled least-squares approach can be taken to optimize the error $\parallel H-{\hat{H}}_r(\cdot ,\theta ){\parallel }_{{\mathcal{H}}_{\mathrm{\infty }}}$ with ${\theta }^{\hspace{0.17em}}:=[{\theta }_N^{\mathrm{\top }},{\theta }_S^{\mathrm{\top }}{]}^{\mathrm{\top }}\in {\mathbb{R}}^{m^2}$, as described in [33,34], and to which we refer for implementation details. We summarize the approach in Algorithm 4.

Algorithm 4:

IHA-PH for pH-DAEs (based on [32])

Input : Large-scale pH-DAE $\mathrm{\Sigma }$ in staircase form; set of interpolation points $\{{\sigma }_1,\dots ,{\sigma }_r\}$ and corresponding right tangential directions $\{b_1,\dots ,b_r\}$ (both closed under complex conjugation); initial parameter vector $\theta \in {\mathbb{R}}^{m^2}$.

Output: Perturbed reduced pH-DAE ${\hat{\mathrm{\Sigma }}}_r$

• 1 Compute the reduced system matrix ${\mathcal{P}}_r$ with steps 1-5 in Algorithm 2.

• 2 Solve ${\theta }^{\ast }=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{}_{\theta \in {\mathbb{R}}^{m^2}}\parallel H-{\hat{H}}_r(\cdot ,\theta ){\parallel }_{{\mathcal{H}}_{\mathrm{\infty }}}$ using the approach in [33,34]

• 3 Compute ${\hat{\mathcal{P}}}_r({\theta }^{\ast })$ as in Lemma 4.2.

• 4 Construct ${\hat{\mathrm{\Sigma }}}_r({\theta }^{\ast })$ with steps 5-10 in Algorithm 1.

4.4. Suitable representations for MOR

Since we apply Galerkin projections to preserve the pH structure, transformations of the FOM under s.s.e. will not change its transfer function, but they will have an impact on MOR, which raises the question of how to find suitable representations of the FOM that yield better approximations. This was recently examined in [35] for explicit ODE systems, and may be incorporated into our framework as follows.

Assume that we have computed the proper system matrix

${\mathcal{P}}_{\mathrm{p}}(s)=\left[\begin{array}{cc}s{E}_{\mathrm{p}}-A_{\mathrm{p}}& -B_{\mathrm{p}}\\ C_{\mathrm{p}}& D_{\mathrm{p}}\end{array}\right],$

of the FOM as in Section 3.1. In Lemma 3.1, we obtained a pH representation of this subsystem by simply decomposing the system matrix into symmetric and skew-symmetric parts. However, there are other ways, and the family of pH representations for this system is parameterized by the Kalman-Yakubovich-Popov (KYP) inequality, as shown in [36]. If the proper system is behaviourally observable, i.e. $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}[(s{E}_{\mathrm{p}}-A_{\mathrm{p}}{)}^{\mathrm{\top }},C_{\mathrm{p}}^{\mathrm{\top }}]=n_2$ for all $s\in \mathbb{C}$, which we assume in the following, the KYP inequality

(32) $\left[\begin{array}{cc}-A_{\mathrm{p}}^{\mathrm{\top }}X-X^{\mathrm{\top }}{A}_{\mathrm{p}}& C_{\mathrm{p}}^{\mathrm{\top }}-X^{\mathrm{\top }}{B}_{\mathrm{p}}\\ C_{\mathrm{p}}-B_{\mathrm{p}}^{\mathrm{\top }}X& D_{\mathrm{p}}+D_{\mathrm{p}}^{\mathrm{\top }}\end{array}\right]\ge 0,X^{\mathrm{T}}{E}_{\mathrm{p}}=E_{\mathrm{p}}^{\mathrm{\top }}X\ge 0,$(32)

has solutions $X\in {\mathbb{R}}^{n_2\times n_2}$ that are bounded such that

$0<E_{\mathrm{p}}^{\mathrm{\top }}{X}_{-}\le E_{\mathrm{p}}^{\mathrm{\top }}X\le E_{\mathrm{p}}^{\mathrm{\top }}{X}_{+}$

with minimal and maximal solutions $X_{-}$ and $X_{+}$, respectively (see [36,Theorem 1]). Now assume that we apply another transformation under s.s.e. on ${\mathcal{P}}_{\mathrm{p}}$ using any $X$ satisfying (32) in the following way

${\tilde{\mathcal{P}}}_{\mathrm{p}}=\left[\begin{array}{cc}{X}^{\mathrm{\top }}& 0\\ 0& I_m\end{array}\right]{\mathcal{P}}_{\mathrm{p}}\left[\begin{array}{cc}{I}_{n_2}& 0\\ 0& I_m\end{array}\right],$

then we obtain the decomposition

$P_{\text{p}}=s{ℰ}_{\text{p}}+{\Gamma }_{\text{p}}+W_{\text{p}},$

with matrices

(33) ${\tilde{E}}_{\mathrm{p}}=\left[\begin{array}{cc}{X}^{\mathrm{\top }}{E}_{\mathrm{p}}& 0\\ 0& 0\end{array}\right],$(33)

(34) ${\tilde{\mathrm{\Gamma }}}_{\mathrm{p}}=\left[\begin{array}{cc}{A}_{\mathrm{p}}^{\mathrm{\top }}X-X^{\mathrm{\top }}{A}_{\mathrm{p}}& -(C_{\mathrm{p}}^{\mathrm{\top }}+X^{\mathrm{\top }}{B}_{\mathrm{p}})\\ C_{\mathrm{p}}+B_{\mathrm{p}}^{\mathrm{\top }}X& D_{\mathrm{p}}-D_{\mathrm{p}}^{\mathrm{T}}\end{array}\right],$(34)

(35) ${\tilde{W}}_{\mathrm{p}}=\left[\begin{array}{cc}-A_{\mathrm{p}}^{\mathrm{\top }}X-X^{\mathrm{\top }}{A}_{\mathrm{p}}& C_{\mathrm{p}}^{\mathrm{\top }}-X^{\mathrm{\top }}{B}_{\mathrm{p}}\\ C_{\mathrm{p}}-B_{\mathrm{p}}^{\mathrm{\top }}X& D_{\mathrm{p}}+D_{\mathrm{p}}^{\mathrm{\top }}\end{array}\right],$(35)

which clearly fulfils the pH structural constraints due to (32). The Hamiltonian of the system associated with ${\tilde{\mathcal{P}}}_{\mathrm{p}}$ changes to

${\tilde{\mathcal{H}}}_{\mathrm{p}}({\tilde{x}}_{\mathrm{p}})=\frac{1}{2}{\tilde{x}}_{\mathrm{p}}^{\mathrm{\top }}{X}^{\mathrm{\top }}{E}_{\mathrm{p}}{\tilde{x}}_{\mathrm{p}}.$

where ${\tilde{x}}_{\mathrm{p}}$ denotes the new state vector of the transformed proper subsystem. As discussed in [35], choosing the minimal solution $X_{-}$ for the transformation is particularly suited for MOR and may significantly improve the approximation quality. In our framework, we can include this transformation change by simply replacing $\mathcal{U}$ in step 3 of Algorithm 1 by

${\mathcal{U}}_{-}=\left[\begin{array}{cc}0& A_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }}\\ X_{-}{\bar{V}}_2& 0\\ -A_{33}^{-\mathrm{\top }}{A}_{23}^{\mathrm{\top }}{X}_{-}{\bar{V}}_2& A_{33}^{-\mathrm{\top }}(C_3^{\mathrm{\top }}-A_{13}^{\mathrm{\top }}{A}_{14}^{-\mathrm{\top }}{C}_4^{\mathrm{\top }})\\ 0& 0\\ 0& I_m\end{array}\right].$

Note that this does not affect the tangential interpolation conditions since we retain $\mathcal{V}$ and also does not affect the matrices $S_{\mathrm{p},r}$, $N_{\mathrm{p},r}$ or $D_{\mathrm{\infty }}$ since we keep the second block column of $\mathcal{U}$ unchanged. However, it does indeed have an effect on the strictly proper part $H_{\mathrm{s}\mathrm{p},r}$ of the ROM’s transfer function between the interpolation points, which is illustrated in the next section by numerical examples. Note that, especially if the matrices $E_{\mathrm{p}},A_{\mathrm{p}}$ are very large or even dense (see Example 3.2), the solution of (32) may be computationally taxing. Therefore, efficient low-rank and/or sparse approximations of $X_{-}$ are required that ideally originate from the sparse FOM matrices, which is an open research problem.

5. Numerical experiments

As initially mentioned, pH-DAE models naturally arise in different engineering fields. A prominent example are electrical RCL networks consisting of linear resistors, capacitors, inductors as well as voltage and current sources. These networks are used, for instance, for the simulation of VLSI interconnect systems or transmission lines. RCL networks are typically modelled using modified nodal analysis (MNA), which is also used by simulation software such as SPICE. This directly yields models of the form

$\begin{array}{cc}E\dot{x}(t)& =(J-R)x(t)+Gu(t),\\ y(t)& =G^{\mathrm{\top }}x(t),\end{array}$

which satisfy the pH structural constraints in Definition 2.2 (see, e.g [37]).

To evaluate our MOR framework, we consider two RCL ladder networks RCL-1M¹ and RCL-12S² whose structure is depicted in Figure 2. The number of RCL loops in the network is denoted by $\bar{n}$ and the number of state-space equations obtained by MNA is $n=3\bar{n}+2m$ where $m$ denotes the number of voltage sources in the network. Both models can be transformed to pH-DAE staircase form via simple permutations of the state-space equations. For this, we use the port-Hamiltonian benchmark collection.³ The staircase dimensions and sigma plots of RCL-1M and RCL-12S are provided in Figure 3.

Figure 2. Structure of the RCL ladder networks modelled by RCL-1M and RCL-12S.

Figure 3. Model parameters and sigma plots for the RCL ladder network models RCL-1M and RCL-12S in pH-DAE staircase form.

RCL-1M is a multiple-input multiple-output (MIMO) version with $\bar{n}=10000$, where the inputs are the voltages of both voltage sources and the outputs are the associated currents. For this model, the red box contains a resistor $R_0$ which leads to a Kronecker index $\nu =1$ and an approximately constant input-output gain for higher frequencies (see Figure 3). In the single-input single-output (SISO) configuration used for RCL-12S, we replace the second voltage source by a wire and only consider the input-to-output behaviour from $u_1$ to $y_1$. Here, we use $\bar{n}=500$ and the red box contains a capacitor $C_0$. This leads to a Kronecker index $\nu =2$ and an improper transfer function. The value of $D_{\mathrm{\infty }}$, i.e. the slope of the transfer function in the high-frequency region is determined by its capacitance.

Let us first consider the reduction of RCL-1M. We reduce the original model to dimension $r=40$ using the algorithms IRKA-PH, TRKSM-PH and IHA-PH. The sigma plots of the resulting ROM transfer functions and respective error systems $H-H_r$ are plotted in Figure 4. For this example, the proper matrices $E_{\mathrm{p}},A_{\mathrm{p}}$ are sparse, and therefore, the residual $\zeta$ in TRKSM-PH can be evaluated very efficiently. This makes TRKSM-PH a computationally efficient alternative to IRKA-PH since it yields a comparable performance while requiring significantly fewer solutions to large-scale linear systems of equations. The additional degrees of freedom in IHA-PH, on the other hand, enable more accurate results in small frequency regions at the expense of a constant error gain at higher frequencies, which results from the perturbation of the reduced feedthrough matrix.

Figure 4. Reduction of the model RCL-1M to order $r=40$ using different interpolatory MOR methods. Given are the sigma plots of the FOM and ROMs (left) and of the respective error systems (right).

For model RCL-12S, we compute ROMs of different dimensions ranging from $r=2$ to $r=20$. The ${\mathcal{H}}_2$ and ${\mathcal{H}}_{\mathrm{\infty }}$ errors are plotted in Figure 5 for each MOR method presented in Section 4. For all methods, the errors decrease for increasing reduced orders $r$, which is expected since more interpolation conditions can be enforced. TRKSM-PH again yields similar ${\mathcal{H}}_2$ errors as IRKA-PH for larger ROM dimensions. For IHA-PH, only the ${\mathcal{H}}_{\mathrm{\infty }}$ errors are plotted in Figure 5 since it produces unbounded ${\mathcal{H}}_2$ errors due to the perturbation of the feedthrough matrix. For most reduced orders, the ${\mathcal{H}}_{\mathrm{\infty }}$ errors are only marginally smaller than those produced by IRKA-PH since the model has only one input-output pair and consequently, $\theta \in \mathbb{R}$. However, as shown for $r=18$, even one additional parameter may improve the ${\mathcal{H}}_{\mathrm{\infty }}$ approximation quality. Finally, we also illustrate the importance of choosing a suitable reduction matrix $\mathcal{U}$. Replacing the matrix $\mathcal{U}$ in IRKA-PH by ${\mathcal{U}}_{-}$, as described in Section 4.4, yields significantly smaller errors both in the ${\mathcal{H}}_2$ and ${\mathcal{H}}_{\mathrm{\infty }}$ norm. Note that a similar basis change could, of course, also be applied to the IHA-PH and TRKSM-PH methods which is expected to yield similar improvements but is omitted here.

Figure 5. Reduction of the model RCL-12S to different reduced orders $r\in \{2,4,\dots ,20\}$. Plotted are the ${\mathcal{H}}_2$ errors (left) and ${\mathcal{H}}_{\mathrm{\infty }}$ errors (right) for different interpolatory MOR methods.

When large-scale RCL networks are reduced, it is crucial that the ROM retains the passivity property of the original model to couple the ROM with other (possibly non-linear) parts of the system. Interpolatory reduction methods that preserve the passivity or MNA structure of RCL models are given by the PRIMA [38] and SPRIM [37] algorithms, respectively. The developed MOR framework in this work extends these methods since the passivity property is inherently preserved by enforcing a pH-DAE structure of the ROM. The models RCL-1M and RCL-12S were also used for an evaluation of optimization-based MOR methods for pH-DAEs in [39] and [21], respectively, to which we refer for a comparison. Similar to these methods, all ROMs in our experiments fulfill the pH structural constraints. This is an advantage compared to passivity-preserving methods for general LTI systems (see, e.g [35,40]) which yield models of the form (5). A subsequent transformation to pH-DAE form requires the solution of a (reduced-order) KYP inequality which may lead to minor violations of the pH structural constraints caused by numerical inaccuracies (see [39, Remark 4]). We also highlight that our MOR framework directly yields ROMs with minimal state-space dimension which may require the solution of discrete-time projected Lyapunov equations for the passivity-preserving approaches proposed in [40].

6. Conclusion

Port-Hamiltonian descriptor systems are particularly suited for the energy-based modelling of multiphysical systems and naturally emerge in different applications. When large-scale pH-DAE models are reduced for simulation or controller design, it is desired to preserve their structural properties and the system dynamics which may be subject to algebraic constraints. In this work, we showed that the Rosenbrock system matrix exhibits a particular structure for pH-DAE models that can be exploited for model reduction. We have deduced a novel interpolatory MOR framework for pH-DAEs in staircase form, which, compared to other structure-preserving interpolation methods, guarantees minimal ROMs in pH-DAE form irrespective of the original system’s Kronecker index. Moreover, its simple structure allows the incorporation of different strategies for choosing suitable interpolation data which were originally proposed for unstructured DAE systems. In applications where pH-DAE models naturally arise, our approach can be considered as an alternative to stability- or passivity-preserving MOR methods since these properties directly follow from the pH structure. As a numerical example, we considered the reduction of electrical circuits which are used, for instance, in VLSI design.

Acknowledgments

We sincerely appreciate the help of Maximilian Bonauer, Nora Reinbold and Paul Schwerdtner with the implementation of the algorithms in Section 4.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) [project number 418612884].

Notes

1. The FOM state-space matrices of RCL-1 M are available at https://doi.org/10.5281/zenodo.6497076.

2. The FOM state-space matrices of RCL-12S are available at https://doi.org/10.5281/zenodo.6602125.

3. https://port-hamiltonian.io.