Model Reduction of Parametric Differential-Algebraic Systems by Balanced Truncation

ABSTRACT

We deduce a procedure to apply balanced truncation to parameter-dependent differential-algebraic systems. For that we solve multiple projected Lyapunov equations for different parameter values to compute the Gramians that are required for the truncation procedure. As this process would lead to high computational costs if we perform it for a large number of parameters, we combine this approach with the reduced basis method that determines a reduced representation of the Lyapunov equation solutions for the parameters of interest. Residual-based error estimators are then used to evaluate the quality of the approximations. After introducing the procedure for a general class of differential-algebraic systems, we turn our focus to systems with a specific structure, for which the method can be applied particularly efficiently. We illustrate the efficiency of our approach on several models from fluid dynamics and mechanics.

KEYWORDS:

• Model order reduction
• reduced basis method
• balanced truncation
• error estimation
• parameter-dependency
• differential-algebraic equations

1. Introduction

In the modelling of various industrial processes one often obtains systems of differential-algebraic equations (DAEs) that are of high dimension. Typical examples are electrical circuits, thermal and diffusion processes, multibody systems, and specific discretized partial differential equations (Brenan, Campbell, and Petzold 1995; Campbell 1980). Often, these systems are further dependent on parameters that can describe physical properties such as material constants and which are often subject to optimization. The resulting parameter-dependent differential-algebraic systems have the form

(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a)

(1b) $\mathbf{y}(\mathit{t})=\mathbf{C}(\mathit{\mu })\mathbf{z}(\mathit{t}),$(1b)

where $\mathbf{E}(\mathit{\mu }),\mathbf{A}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$, $\mathbf{B}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},\mathit{m}}$, and $\mathbf{C}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{p},\mathit{N}}$ are dependent on parameters $\mathit{\mu }\in \mathcal{D}$, where $\mathcal{D}\subset {\mathbb{R}}^d$. The input, the state, and the output are denoted by $\mathbf{u}(\mathit{t})\in {\mathbb{R}}^m$, $\mathbf{z}(\mathit{t})\in {\mathbb{R}}^N$ and $\mathbf{y}(\mathit{t})\in {\mathbb{R}}^p$ for $\mathit{t}\in [0,\mathrm{\infty })$. The matrix $\mathbf{E}(\mathit{\mu })$ can be a singular. However, throughout this paper, we assume uniform regularity, i.e., the pencil $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$, which is a matrix polynomial in the variable $\mathit{s}$, is assumed to be regular for all $\mathit{\mu }\in \mathcal{D}$, that is $det(\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu }))$ is not the zero polynomial for all $\mathit{\mu }\in \mathcal{D}$. Moreover, we assume that the system under consideration is uniformly asymptotically stable, i.e., the system (or the pencil $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$) is asymptotically stable for all $\mathit{\mu }\in \mathcal{D}$, meaning that all finite eigenvalues of the pencil $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$ have negative real part for all $\mathit{\mu }\in \mathcal{D}$. Throughout the derivations provided in this work, we assume that $\mathbf{E}(\cdot )$, $\mathbf{A}(\cdot )$, $\mathbf{B}(\cdot )$, and $\mathbf{C}(\cdot )$ are affine in the parameter $\mathit{\mu }$, i.e.,

(2) $\begin{array}{cc}\mathbf{E}(\mathit{\mu })& ={\mathbf{E}}_0+\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{E}}}{\mathrm{\Theta }}_k^{\mathbf{E}}(\mathit{\mu }){\mathbf{E}}_k,\mathbf{A}(\mathit{\mu })={\mathbf{A}}_0+\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{A}}}{\mathrm{\Theta }}_k^{\mathbf{A}}(\mathit{\mu }){\mathbf{A}}_k,\\ \\ \mathbf{B}(\mathit{\mu })& ={\mathbf{B}}_0+\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{B}}}{\mathrm{\Theta }}_k^{\mathbf{B}}(\mathit{\mu }){\mathbf{B}}_k,\mathbf{C}(\mathit{\mu })={\mathbf{C}}_0+\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{C}}}{\mathrm{\Theta }}_k^{\mathbf{C}}(\mathit{\mu }){\mathbf{C}}_k,\end{array}$(2)

where all ${\mathrm{\Theta }}_k^{\mathbf{E}},{\mathrm{\Theta }}_k^{\mathbf{A}},{\mathrm{\Theta }}_k^{\mathbf{B}},{\mathrm{\Theta }}_k^{\mathbf{C}}:\mathcal{D}\to \mathbb{R}$ are parameter-dependent continuous functions and ${\mathbf{E}}_k$, ${\mathbf{A}}_k$, ${\mathbf{B}}_k$, ${\mathbf{C}}_k$ are parameter-independent matrices. For reasons of computational efficiency, we always assume that $n_{\mathbf{E}},n_{\mathbf{A}},n_{\mathbf{B}},n_{\mathbf{C}}\ll \mathit{N}$ and we also assume that the number of inputs and outputs is small compared to the dimension of the states, i.e., $\mathit{m},\mathit{p}\ll \mathit{N}$. For ease of notation, it is sometimes handy to put the matrices ${\mathbf{E}}_0$, ${\mathbf{A}}_0$, ${\mathbf{B}}_0$, and ${\mathbf{C}}_0$ into the sums in (3) and multiply them with the factors ${\mathrm{\Theta }}_0^{\mathbf{E}}={\mathrm{\Theta }}_0^{\mathbf{A}}={\mathrm{\Theta }}_0^{\mathbf{B}}={\mathrm{\Theta }}_0^{\mathbf{C}}\equiv 1$. The model (1) is also referred to as the full-order model (FOM). Often, one also considers the input/output mapping of the system in the frequency domain. This relation is typically expressed by the transfer function $\mathcal{G}(\mathit{\mu },\mathit{s}):=\mathbf{C}(\mathit{\mu })(\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu }){)}^{-1}\mathbf{B}(\mathit{\mu }).$

Because of the high state-space dimension $\mathit{N}$ of (1) it is useful to apply reduction methods to extract the essential information of the system and its solution. More precisely, we want to determine a (uniformly regular and asymptotically stable) reduced-order model (ROM)

(3) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{E}}_{\mathrm{R}}(\mathit{\mu }){\mathbf{z}}_{\mathrm{R}}(\mathit{t})& ={\mathbf{A}}_{\mathrm{R}}(\mathit{\mu }){\mathbf{z}}_{\mathrm{R}}(\mathit{t})+{\mathbf{B}}_{\mathrm{R}}(\mathit{\mu })\mathbf{u}(\mathit{t}),\\ \\ {\mathbf{y}}_{\mathrm{R}}(\mathit{t})& ={\mathbf{C}}_{\mathrm{R}}(\mathit{\mu }){\mathbf{z}}_{\mathrm{R}}(\mathit{t}),\end{array}$(3)

with ${\mathbf{E}}_{\mathrm{R}}(\mathit{\mu }),{\mathbf{A}}_{\mathrm{R}}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{r},\mathit{r}}$, ${\mathbf{B}}_{\mathrm{R}}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{r},\mathit{m}}$, and ${\mathbf{C}}_{\mathrm{R}}\in {\mathbb{R}}^{\mathit{p},\mathit{r}}$, where $\mathit{r}\ll \mathit{N}$. The reduced state and output are ${\mathbf{z}}_{\mathrm{R}}(\mathit{t})\in {\mathbb{R}}^r$ and ${\mathbf{y}}_{\mathrm{R}}(\mathit{t})\in {\mathbb{R}}^p$. The aim is to determine the surrogate model in such a way that the output of the ROM well-approximates the output of the FOM for all admissible inputs $\mathbf{u}(\cdot )$ and all parameters $\mathit{\mu }\in \mathcal{D}$. In terms of the transfer function, this means that the error $∥\mathcal{G}(\mathit{\mu },\cdot )-{\mathcal{G}}_{\mathrm{R}}(\mathit{\mu },\cdot )∥$ is sufficiently small in an appropriate norm for all $\mathit{\mu }\in \mathcal{D}$, where ${\mathcal{G}}_{\mathrm{R}}(\mathit{\mu },\mathit{s})$ denotes the transfer function of the ROM.

For parameter-independent problems, i.e.,

$\mathit{E}(\mathit{\mu })\equiv \mathit{E},\mathit{A}(\mathit{\mu })\equiv \mathit{A},\mathit{B}(\mathit{\mu })\equiv \mathit{B},\mathit{C}(\mathit{\mu })\equiv \mathit{C},\mathit{G}(\mathit{\mu },\mathit{s})\equiv \mathit{G}(\mathit{s})$

there are several classes of model order reduction techniques. Examples are singular value-based approaches like balanced truncation (Benner et al. 2017; Moore 1981; Tombs and Postlethwaite 1987) and Hankel norm approximations (Glover 1984). Additionally, there are Krylov subspace-based methods such as the iterative rational Krylov algorithm (IRKA) (Benner et al. 2017; Flagg, Beattie, and Gugercin 2012; Gugercin, Antoulas, and Beattie 2008; Gugercin, Stykel, and Wyatt 2013) and moment matching as well as data-driven methods such as the Loewner framework (Beattie, Gugercin, and Mehrmann 2022; Gosea, Poussot-Vassal, and Antoulas 2021; Mayo and Antoulas 2007). Corresponding methods for parameter-dependent systems are introduced in (Baur, Benner, and Feng 2014; Benner et al. 2016; Geuss, Panzer, and Lohmann 2013).

In this article, we focus on balanced truncation which is one of the most popular reduction techniques. This is mainly due to the guaranteed asymptotic stability of the ROM, the existence of an error bound, and appealing numerical techniques for the involved Lyapunov equations (Benner and Schneider 2013; Druskin, Simoncini 2011; Schmidt, Wittwar, and Haasdonk 2020; Son and Stykel 2017). Within balanced truncation, Lyapunov equations corresponding to the original system need to be solved. The solutions of these equations are called Gramians. They describe the input-to-state and state-to-output behavior and are used to construct projection matrices for the reduction. The multiplication of the system matrices with these projection matrices then results in a ROM.

All the above-mentioned methods focus on the case $\mathbf{E}=I_N$ and are, however, not directly applicable in the DAE case. Even if the problem is not parameter-dependent, there are several challenges that one has to face (here we assume for simplicity, that the problem is parameter-independent):

a) Since the matrix $\mathbf{E}$ is typically singular, the transfer function $\mathcal{G}(\mathit{s})$ is possibly improper, i.e., it may have a polynomial part which is unbounded for growing values of $|\mathit{s}|$. If the reduced transfer function $\mathbf{z}(\mathit{t})={\mathbf{z}}_{\mathrm{p}}(\mathit{t})+{\mathbf{z}}_{\mathrm{i}}(\mathit{t})$ does not match this polynomial part, then the error transfer function $N_{\mathrm{f}}$ is not an element of the rational Hardy spaces $\mathcal{R}{\mathcal{H}}_{\mathrm{\infty }}^{p,m}$ or $\mathcal{R}{\mathcal{H}}_2^{p,m}$ (see, e.g., Zhou, Doyle, and Glover 1996) for a definition) and the output error cannot be bounded by the input norm. Thus, a model reduction scheme for DAEs must preserve the polynomial part of its transfer function. This is addressed in (Gugercin, Stykel, and Wyatt 2013; Mehrmann and Stykel 2005; Stykel 2004).

b) In balanced truncation, one has to solve two large-scale Lyapunov equations. In the DAE case, one has to consider a generalization of these – so-called projected Lyapunov equations. Thereby, the Lyapunov equations are projected onto the left and right deflating subspaces corresponding to the finite or infinite eigenvalues of the matrix pencil $\mathit{s}\mathbf{E}-\mathbf{A}$. This involves specific projectors which are hard to form explicitly and that might destroy the sparsity structure of the coefficient matrices. However, for DAE systems of special structure, the projectors are of a particular form which can be exploited from the numerical point of view. For details, we refer to (Benner, Saak, and Uddin 2016; Heinkenschloss, Sorensen, and Sun 2008; Saak and Voigt 2018; Stykel 2008).

Besides the approaches mentioned above, some of the issues in model reduction of DAEs are also addressed in (Alì et al. 2014; Banagaaya, Ali, and Schilders 2016). Here, the authors introduce an index-aware model order reduction scheme that splits the DAE system into a differential and an algebraic part. However, in this method, it is often difficult to maintain sparsity of the system matrices, and hence, it is in general not feasible for many systems with large dimensions. Also data-driven approaches have been recently extended to differential-algebraic systems, see, e.g. (Antoulas, Gosea, and Heinkenschloss 2020; Gosea, Zhang, and Antoulas 2020). Finally, we would like to mention (Benner and Stykel 2017) which provides an overview of model order reduction for DAE systems.

Different approaches for parametric model reduction have been developed as well. Besides others, there exist two main classes of such methods. First, one could determine global bases for the projection spaces for the entire parametric model, e.g., by concatenating the local bases for different parameter values, see, e.g. (Baur et al. 2011). Second, there exist approaches for constructing local bases for a specific parameter value by combining local basis information of several prereduced models at different parameter values, e.g., by interpolating between different models (Eid et al. 2011; Geuss, Panzer, and Lohmann 2013). We also refer to (Benner, Gugercin, and Willcox 2015) for an extensive survey on such methods.

However, to the best of our knowledge, none of the techniques presented in (Benner, Gugercin, and Willcox 2015) has already been extended to the case of DAE systems. Therefore, the motivation of this paper is to develop a first method for this purpose based on balanced truncation. If we aim to reduce a parameter-dependent system by balanced truncation in a naive fashion, we would have to solve the Lyapunov equations for each individual parameter of interest, which is computationally prohibitive. To avoid these computational costs, for the case $\mathbf{E}=I_N$, the authors in (Son et al. 2021) apply some interpolation techniques to approximate the Gramians for all admissible parameters. The authors in (Son and Stykel 2017), on the other hand, apply the reduced basis method which is a well-established method to reduce parameter-dependent partial differential equations (Hesthaven, Rozza, and Stamm 2016; Kunkel and Mehrmann 2006). Using this method, the authors determine the approximate subspaces in which the solutions of the corresponding Lyapunov equations are contained in.

In this paper, we extend the reduced basis method to the case of DAE systems where the parametric projected Lyapunov equations are only solved for few sampling parameters. Then, based on these solutions, a reduced subspace in which the Lyapunov equation solutions for all $\mathit{\mu }\in \mathcal{D}$ approximately live, is constructed. The latter steps form the computationally expensive offline phase. After that, using the reduced basis representation, the Lyapunov equations can be solved much more efficiently for all $\mathit{\mu }\in \mathcal{D}$ in the online phase. A crucial question in the offline phase is the choice of the sample parameters. Usually, a grid of test parameters is selected. For these, the error is quantified using a posteriori error estimators. Then, new samples are taken at those parameters at which the error estimator estimates the largest error.

In this paper, we generalize the reduced basis balanced truncation method of (Son and Stykel 2017) to differential-algebraic systems with a focus on structured systems from specific application areas. The main problems we solve in this paper are the following:

a) We derive error estimators for parameter-dependent projected Lyapunov equations. These can be improved if the given system is uniformly strictly dissipative. Since this condition is not always satisfied, we discuss transformations to achieve this for systems of a particular index 3 structure.

b) We discuss in detail how the polynomial part of parameter-dependent transfer functions of DAE systems can be extracted by efficiently approximating their Markov parameters.

c) We apply this approach to several application problems and illustrate its effectiveness.

Our method also follows the paradigm of decomposing the procedure into an offline and online phase. In the offline phase, we approximate the image spaces of controllability and observability Gramians that are the solutions of projected Lyapunov equations. These approximate image spaces are then used in the online phase, to determine approximations of the Gramians, and hence, the reduced order model, efficiently for all required parameters. The online phase can be performed for any parameter $\mathit{\mu }\in \mathcal{D}$. Due to the reduced dimension obtained by the offline phase, this step is very fast and can be performed repeatedly for all required parameters. Moreover, in this work, we introduce new error estimators, which are used in the offline phase to assess the quality of the reduction.

This paper is organized as follows. In Section 2, two model problems are introduced that motivate the method presented in this paper. In Section 3, we review existing concepts and approaches that will be used later in this paper. Thereby, in Subsection 3.1, we recall projection techniques to decouple the differential and algebraic equations in (1(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a) ). Afterwards, in Subsection 3.2, we consider model order reduction by balanced truncation for DAE systems. We further address the numerical challenges that arise in computing the solutions of the required Lyapunov equations. Since solving Lyapunov equations for every requested parameter leads to high computational costs, in Section 4 the reduced basis method is presented, which was first applied to standard Lyapunov equations in (Son and Stykel 2017). We derive two error estimators for our method, one of them is motivated by the estimator from (Son and Stykel 2017), the other one is an adaption of the estimator presented in (Schmidt, Wittwar, and Haasdonk 2020). Here, we also discuss the treatment of the algebraic part of the system that may lead to polynomial parts in the transfer function. Finally, in Section 5, we evaluate the method of this paper by applying it to our model problems from Section 2.

2. Model problems

2.1. Problem 1: Stokes-like DAEs

The first example is the system

(4) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{cc}\mathit{E}(\mathit{\mu })& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}\mathit{x}(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right]& =\left[\begin{array}{cc}\mathit{A}(\mathit{\mu })& \mathit{G}(\mathit{\mu })\\ \mathit{G}{(\mathit{\mu })}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}\mathit{x}(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right]+\left[\begin{array}{c}{B}_1(\mathit{\mu })\\ B_2(\mathit{\mu })\end{array}\right]\mathit{u}(\mathit{t}),\\ \\ \mathit{y}(\mathit{t})& =\left[\begin{array}{cc}{C}_1(\mathit{\mu })& C_2(\mathit{\mu })\end{array}\right]\left[\begin{array}{c}\mathit{x}(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right],\end{array}$(4)

which arises, e.g., if we discretize the incompressible Navier–Stokes equation. The parameter-independent version is presented in (Mehrmann and Stykel 2005; Stykel 2004). The system matrices are dependent on a parameter $\mathit{\mu }\in \mathcal{D}$, where $\mathcal{D}\subset {\mathbb{R}}^d$ is the parameter domain. For fluid-flow problems, the matrices $\mathit{E}(\mathit{\mu }),\mathit{A}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{n},\mathit{n}}$ represent the masses and the discretized Laplace operator. Naturally, it holds that $\mathit{E}(\mathit{\mu })=\mathit{E}(\mathit{\mu }{)}^{\mathrm{T}}>0$ and $\mathit{A}(\mathit{\mu })=\mathit{A}(\mathit{\mu }{)}^{\mathrm{T}}<0$ for all $\mathit{\mu }\in \mathcal{D}$, where $\mathit{H}>0$ ($\mathit{H}\ge 0$) denotes a positive (semi)definite matrix $\mathit{H}$ and correspondingly $\mathit{H}<0$ ($\mathit{H}\le 0$) denotes a negative (semi)definite matrix. The discrete gradient is given by $\mathit{G}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{n},\mathit{q}}$ which we assume to be of full column rank for all $\mathit{\mu }\in \mathcal{D}$. The matrices $B_1(\mathit{\mu })\in {\mathbb{R}}^{\mathit{n},\mathit{m}},B_2(\mathit{\mu })\in {\mathbb{R}}^{\mathit{q},\mathit{m}}$ and $C_1(\mathit{\mu })\in {\mathbb{R}}^{\mathit{p},\mathit{n}},C_2(\mathit{\mu })\in {\mathbb{R}}^{\mathit{p},\mathit{q}}$ are the input and the output matrices, respectively. The state of the system at time $\mathit{t}$ is given by $\mathit{x}(\mathit{t})\in {\mathbb{R}}^n$ and $\mathit{\lambda }(\mathit{t})\in {\mathbb{R}}^q$. The vectors $\mathit{u}(\mathit{t})\in {\mathbb{R}}^m$ and $\mathit{y}(\mathit{t})\in {\mathbb{R}}^p$ are the input and output of the system.

In view of (2) we assume that the symmetry and definiteness of $\mathit{E}(\cdot )$ and $\mathit{A}(\cdot )$ is represented in the affine parameter structure as

$\mathit{E}(\mathit{\mu })=\sum _{\mathit{k}=0}^{n_E}{\Theta }_k^E(\mathit{\mu })E_k\mathrm{and}\mathit{A}(\mathit{\mu })=\sum _{\mathit{k}=0}^{n_A}{\Theta }_k^A(\mathit{\mu })A_k,$

where ${\Theta }_k^E(\cdot )$, ${\Theta }_k^A(\cdot )$ are positive parameter-dependent continuous functions with the convention ${\mathrm{\Theta }}_0^E={\mathrm{\Theta }}_0^A\equiv 1$. Furthermore, $E_0=E_0^{\mathrm{T}}>0$, $A_0=A_0^{\mathrm{T}}<0$, $E_k=E_k^{\mathrm{T}}\ge 0$, $A_k=A_k^{\mathrm{T}}\le 0$ for $\mathit{k}\ge 1$ are parameter-independent matrices. Moreover, for reasons of computational efficiency, we always assume that $n_E,n_A\ll \mathit{n}$.

2.2. Problem 2: Mechanical systems

The second system we consider is of the form

(5) $\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{ccc}{I}_{n_x}& 0& 0\\ 0& \mathit{M}(\mathit{\mu })& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1(\mathit{t})\\ x_2(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right]=\left[\begin{array}{ccc}0& I_{n_x}& 0\\ -\mathit{K}(\mathit{\mu })& -\mathit{D}(\mathit{\mu })& \mathit{G}(\mathit{\mu })\\ \mathit{G}{(\mathit{\mu })}^{\mathrm{T}}& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1(\mathit{t})\\ x_2(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right]+\left[\begin{array}{c}0\\ B_x(\mathit{\mu })\\ 0\end{array}\right]\mathit{u}(\mathit{t}),\text{break}\mathit{y}(\mathit{t})=\left[\begin{array}{ccc}{C}_x(\mathit{\mu })& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1(\mathit{t})\\ x_2(\mathit{t})\\ \mathit{\lambda }(\mathit{t})\end{array}\right],\end{array}$(5)

which results from a linearization of the spring-mass-damper model presented in (Mehrmann and Stykel 2005). The mass, damping, and stiffness matrices $\mathit{M}(\mathit{\mu }),\mathit{D}(\mathit{\mu }),\mathit{K}(\mathit{\mu })\in {\mathbb{R}}^{n_x,n_x}$ are assumed to be symmetric and positive definite for all $\mathit{\mu }\in \mathcal{D}$. The matrices $B_x(\mathit{\mu })\in {\mathbb{R}}^{n_x,\mathit{m}}$ and $C_x(\mathit{\mu })\in {\mathbb{R}}^{p,n_x}$ are the input and the output matrices. The matrix $\mathit{G}(\mathit{\mu })\in {\mathbb{R}}^{n_x,\mathit{q}}$ reflects algebraic constraints on the system and is of full column rank. In this example, the state includes the displacement $x_1(\mathit{t})\in {\mathbb{R}}^{n_x}$, the velocity $x_2(\mathit{t})\in {\mathbb{R}}^{n_x}$, and Lagrange multiplier $\mathit{\lambda }(\mathit{t})\in {\mathbb{R}}^q$.

For convenience, we also define

(6) $\begin{array}{cc}\mathit{E}(\mathit{\mu })& :=\left[\begin{array}{cc}{I}_{n_x}& 0\\ 0& \mathit{M}(\mathit{\mu })\end{array}\right],\mathit{A}(\mathit{\mu }):=\left[\begin{array}{cc}0& I_{n_x}\\ -\mathit{K}(\mathit{\mu })& -\mathit{D}(\mathit{\mu })\end{array}\right],\mathit{B}(\mathit{\mu }):=\left[\begin{array}{c}0\\ B_x(\mathit{\mu })\end{array}\right],\\ \mathit{C}(\mathit{\mu })& :=\left[\begin{array}{cc}{C}_x(\mathit{\mu })& 0\end{array}\right],\mathit{n}:=2n_x.\end{array}$(6)

Then with this notation, we can write the mechanical system similarly as in (4) with the difference that the off-diagonal blocks of the state matrix are not the transposes of each other.

Similarly to the first model problem, we assume that $\mathit{M}(\cdot )$, $\mathit{D}(\cdot )$, and $\mathit{K}(\cdot )$ can be written as

$\mathit{M}(\mathit{\mu })=\sum _{\mathit{k}=0}^{n_M}{\Theta }_k^M(\mathit{\mu })M_k,\mathit{D}(\mathit{\mu })=\sum _{\mathit{k}=0}^{n_D}{\Theta }_k^D(\mathit{\mu })D_k,\mathrm{and}\mathit{K}(\mathit{\mu })=\sum _{\mathit{k}=0}^{n_K}{\Theta }_k^K(\mathit{\mu })K_k,$

where ${\Theta }_k^M(\cdot )$, ${\Theta }_k^D(\cdot )$, ${\Theta }_k^K(\cdot )$ are positive parameter-dependent continuous functions using the convention ${\mathrm{\Theta }}_0^M={\mathrm{\Theta }}_0^D={\mathrm{\Theta }}_0^K\equiv 1$. Furthermore, $M_0=M_0^{\mathrm{T}}>0$, $D_0=D_0^{\mathrm{T}}>0$, $K_0=K_0^{\mathrm{T}}>0$ and $M_k=M_k^{\mathrm{T}}\ge 0$, $D_k=D_k^{\mathrm{T}}\ge 0$, $K_k=K_k^{\mathrm{T}}\ge 0$ for $\mathit{k}\ge 1$ are parameter-independent matrices with $n_M,n_D,n_K\ll n_x$.

3. Preliminaries

3.1. Decoupling of differential and algebraic equations

In this section, we briefly describe the projection technique for decoupling the algebraic constraints in DAEs of the form (1) from the remaining differential equations for a fixed parameter $\mathit{\mu }\in \mathcal{D}$ and under the assumptions of regularity and asymptotic stability. Therefore, for simplicity of presentation, we omit $\mathit{\mu }$ in this section. The details of this can be found in (Heinkenschloss, Sorensen, and Sun 2008; Saak and Voigt 2018; Stykel 2004, 2006).

As shown in (Berger, Ilchmann, and Trenn 2012) there exist regular matrices $\mathbf{W},\mathbf{T}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$ that transform the system matrices into quasi-Weierstraß form, that is

(7) $\mathbf{E}=\mathbf{W}\left[\begin{array}{cc}{I}_{{\mathit{N}}_{\mathrm{f}}}& 0\\ 0& \mathbf{N}\end{array}\right]\mathbf{T},\mathbf{A}=\mathbf{W}\left[\begin{array}{cc}\mathbf{J}& 0\\ 0& I_{N_{\mathrm{\infty }}}\end{array}\right]\mathbf{T},$(7)

where $\mathbf{J}\in {\mathbb{R}}^{{\mathit{N}}_{\mathrm{f}},N_{\mathrm{f}}}$ and $\mathbf{N}\in {\mathbb{R}}^{N_{\mathrm{\infty }},N_{\mathrm{\infty }}}$ is nilpotent with nilpotency index $\mathit{\nu }$ that coincides with the (Kronecker) index of the DAE system (1). Here, $N_{\mathrm{f}}$ is the number of the finite eigenvalues of the matrix pencil $\mathit{s}\mathbf{E}-\mathbf{A}$ and $N_{\mathrm{\infty }}$ is the number of the infinite ones. Accordingly, spectral projectors onto the left and right deflating subspace of $\mathit{s}\mathbf{E}-\mathbf{A}$ corresponding to the $N_{\mathrm{f}}$ finite eigenvalues are defined as

${\mathbf{\Pi }}_{\mathrm{l}}:=\mathbf{W}\left[\begin{array}{cc}{I}_{{\mathit{N}}_{\mathrm{f}}}& 0\\ 0& 0\end{array}\right]{\mathbf{W}}^{-1},{\mathbf{\Pi }}_{\mathrm{r}}:={\mathbf{T}}^{-1}\left[\begin{array}{cc}{I}_{{\mathit{N}}_{\mathrm{f}}}& 0\\ 0& 0\end{array}\right]\mathbf{T}.$

We define

${\mathbf{F}}_{\mathbf{J}}(\mathit{t}):={\mathbf{T}}^{-1}\left[\begin{array}{cc}{\mathrm{e}}^{\mathbf{J}t}& 0\\ 0& 0\end{array}\right]{\mathbf{W}}^{-1},\mathit{t}\in \mathbb{R}\mathrm{and}$

${\mathbf{F}}_{\mathbf{N}}(\mathit{k}):={\mathbf{T}}^{-1}\left[\begin{array}{cc}0& 0\\ 0& -{\mathbf{N}}^k\end{array}\right]{\mathbf{W}}^{-1},\mathit{k}=0,\dots ,\mathit{\nu }-1$

to derive the states

${\mathbf{z}}_{\mathrm{p}}(\mathit{t}):={\int }_0^t{\mathbf{F}}_{\mathbf{J}}(\mathit{t}-\mathit{\tau })\mathbf{Bu}(\mathit{\tau })\mathrm{d}\mathit{\tau },{\mathbf{z}}_{\mathrm{i}}(\mathit{t}):=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{u}}^{(\mathit{k})}(\mathit{t}),$

under the assumption that the input $\mathbf{u}(\cdot )$ is sufficiently smooth and that the consistency conditions are satisfied, i.e., we assume that ${\mathbf{z}}_{\mathrm{i}}(0):={\sum }_{k=0}^{\nu -1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{u}}^{(\mathit{k})}(0)$. They satisfy $\mathbf{z}(\mathit{t})={\mathbf{z}}_{\mathrm{p}}(\mathit{t})+{\mathbf{z}}_{\mathrm{i}}(\mathit{t})$ such that $\mathbf{z}(\mathit{t})$ solves Equation (1a)(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a) . We refer to ${\mathbf{z}}_{\mathrm{p}}(\mathit{t})$ and ${\mathbf{z}}_{\mathrm{i}}(\mathit{t})$ as differential and algebraic states. These states are associated with the proper and improper controllability Gramians

(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8)

respectively. They span the reachable subspaces of the differential and algebraic states. Analogously, the proper and improper observability Gramians

${\mathbf{Q}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathrm{d}\mathit{t},{\mathbf{Q}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}^{\mathrm{T}}(\mathit{k}){\mathbf{C}}^{\mathrm{T}}\mathbf{C}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})$

span the corresponding observability spaces. The proper Gramians are obtained by solving the projected continuous-time Lyapunov equations

(9) $\mathbf{E}{\mathbf{P}}_{\mathrm{p}}{\mathbf{A}}^{\mathrm{T}}+\mathbf{A}{\mathbf{P}}_{\mathrm{p}}{\mathbf{E}}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}},{\mathbf{P}}_{\mathrm{p}}={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{p}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(9)

(10) ${\mathbf{E}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{A}+{\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{E}=-{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}{\mathbf{\Pi }}_{\mathrm{r}},{\mathbf{Q}}_{\mathrm{p}}={{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}{\mathbf{\Pi }}_{\mathrm{l}}$(10)

and the improper ones by solving the discrete-time Lyapunov equations

(11) $\mathbf{A}{\mathbf{P}}_{\mathrm{i}}{\mathbf{A}}^{\mathrm{T}}-\mathbf{E}{\mathbf{P}}_{\mathrm{i}}{\mathbf{E}}^{\mathrm{T}}=(I_N-{\mathbf{\Pi }}_{\mathrm{l}})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}(I_N-{\mathbf{\Pi }}_{\mathrm{l}}{)}^{\mathrm{T}},0={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{i}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(11)

(12) ${\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}\mathbf{A}-{\mathbf{E}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}\mathbf{E}=(I_N-{\mathbf{\Pi }}_{\mathrm{r}}{)}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}(I_N-{\mathbf{\Pi }}_{\mathrm{r}}),0={{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}{\mathbf{\Pi }}_{\mathrm{l}}.$(12)

In Subsection 3.2.2 we will describe how these projected Lyapunov equations can be solved.

In practice, the projectors ${\mathbf{\Pi }}_{\mathrm{l}}$ and ${\mathbf{\Pi }}_{\mathrm{r}}$ are computationally expensive to form and even more expensive to compute with as they are dense matrices of large dimension in general. To determine them, the derivations from (Kunkel and Mehrmann 2006) can be applied, however they are not numerically stable since a Jordan canonical form needs to be computed. On the other hand, the methods presented in (Banagaaya 2014; Campbell, Meyer, and Rose 1976), or (März 1996) might be used to generate the differential and algebraic state spaces, however, they involve the determination of matrix kernels which are based on matrix decompositions. Suitable parts of the matrices involved in these decompositions are then used to determine the projection matrices ${\mathbf{\Pi }}_{\mathrm{l}}$ and ${\mathbf{\Pi }}_{\mathrm{r}}$. They are in general not of a sparse structure and, therefore, not applicable if we consider matrices of large dimensions. Hence, if possible, it is advantageous to determine the projection matrices analytically. In many application fields, this is possible because the matrices $\mathbf{E}$ and $\mathbf{A}$ are of special structure that can be utilized for this purpose. The examples in Section 2 have such a numerically advantageous structure. Further examples can be found in (März 1996; Stykel 2008).

3.2. Model reduction of differential-algebraic systems

The aim of this section is to review methods to reduce the DAE system (1(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a) ) for a fixed parameter. We utilize balanced truncation modified for projected systems which is presented in Subsection 3.2.1. Afterwards, in Subsection 3.2.2, we summarize the ADI and Smith methods to solve the involved projected Lyapunov equations. Finally, in Subsection 3.2.3, we present an alternative approach to approximate the algebraic parts of the system (potentially leading to polynomial parts in the transfer function) which is inspired by (Schwerdtner et al. 2023).

3.2.1. Balanced truncation

The aim of this subsection is to find a reduced system

(13) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{E}}_{\mathrm{R}}{\mathbf{z}}_{\mathrm{R}}(\mathit{t})& ={\mathbf{A}}_{\mathrm{R}}{\mathbf{z}}_{\mathrm{R}}(\mathit{t})+{\mathbf{B}}_{\mathrm{R}}\mathbf{u}(\mathit{t}),\\ {\mathbf{y}}_{\mathrm{R}}(\mathit{t})& ={\mathbf{C}}_{\mathrm{R}}{\mathbf{z}}_{\mathrm{R}}(\mathit{t}),\end{array}$(13)

with ${\mathbf{E}}_{\mathrm{R}},{\mathbf{A}}_{\mathrm{R}}\in {\mathbb{R}}^{\mathit{r},\mathit{r}},{\mathbf{B}}_{\mathrm{R}}\in {\mathbb{R}}^{\mathit{r},\mathit{m}},{\mathbf{C}}_{\mathrm{R}}\in {\mathbb{R}}^{\mathit{p},\mathit{r}}$, $\mathit{r}\ll \mathit{N}$ which approximates the input/output behaviour of the original system (1(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a) ) (for a fixed parameter).

Since the different Gramians are positive semidefinite, there exist factorizations

${\mathbf{P}}_{\mathrm{p}}={\mathbf{R}}_{\mathrm{p}}{\mathbf{R}}_{\mathrm{p}}^{\mathrm{T}},{\mathbf{P}}_{\mathrm{i}}={\mathbf{R}}_{\mathrm{i}}{\mathbf{R}}_{\mathrm{i}}^{\mathrm{T}},{\mathbf{Q}}_{\mathrm{p}}={\mathbf{S}}_{\mathrm{p}}{\mathbf{S}}_{\mathrm{p}}^{\mathrm{T}},{\mathbf{Q}}_{\mathrm{i}}={\mathbf{S}}_{\mathrm{i}}{\mathbf{S}}_{\mathrm{i}}^{\mathrm{T}}$

with factors ${\mathbf{R}}_{\mathrm{p}},{\mathbf{R}}_{\mathrm{i}},{\mathbf{S}}_{\mathrm{p}},{\mathbf{S}}_{\mathrm{i}}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$. We consider the singular value decompositions

${\mathbf{S}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{E}{\mathbf{R}}_{\mathrm{p}}={\mathbf{U}}_{\mathrm{p}}{\mathrm{\Sigma }}_{\mathrm{p}}{\mathbf{V}}_{\mathrm{p}}^{\mathrm{T}}=\left[\begin{array}{cc}{\mathbf{U}}_{\mathrm{p},1}& {\mathbf{U}}_{\mathrm{p},2}\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathrm{p},1}& 0\\ 0& {\mathrm{\Sigma }}_{\mathrm{p},2}\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathrm{p},1}^{\mathrm{T}}\\ {\mathbf{V}}_{\mathrm{p},2}^{\mathrm{T}}\end{array}\right],$

(14) ${\mathbf{S}}_{\mathrm{i}}^{\mathrm{T}}\mathbf{A}{\mathbf{R}}_{\mathrm{i}}={\mathbf{U}}_{\mathrm{i}}{\mathrm{\Sigma }}_{\mathrm{i}}{\mathbf{V}}_{\mathrm{i}}^{\mathrm{T}}=\left[\begin{array}{cc}{\mathbf{U}}_{\mathrm{i},1}& {\mathbf{U}}_{\mathrm{i},2}\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathrm{i},1}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathrm{i},1}^{\mathrm{T}}\\ {\mathbf{V}}_{\mathrm{i},2}^{\mathrm{T}}\end{array}\right],$(14)

where the matrix ${\mathrm{\Sigma }}_{\mathrm{p}}=\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_N)$ is a diagonal matrix with nonincreasing nonnegative entries that are the proper Hankel singular values and ${\mathrm{\Sigma }}_{\mathrm{i},1}$ is a diagonal matrix including the improper nonzero Hankel singular values. With the matrix ${\mathrm{\Sigma }}_{\mathrm{p},1}$ which contains the $r_{\mathrm{p}}$ largest Hankel singular values and ${\mathrm{\Sigma }}_{\mathrm{i},1}\in {\mathbb{R}}^{{\mathit{r}}_{\mathrm{i}},r_{\mathrm{i}}}$ we construct the left and right projection matrices

${\mathbf{W}}_{\mathrm{R}}:=\left[\begin{array}{cc}{\mathbf{S}}_{\mathrm{p}}{\mathbf{U}}_{\mathrm{p},1}{\mathrm{\Sigma }}_{\mathrm{p},1}^{-12}& {\mathbf{S}}_{\mathrm{i}}{\mathbf{U}}_{\mathrm{i},1}{\mathrm{\Sigma }}_{\mathrm{i},1}^{-12}\end{array}\right]\mathrm{and}{\mathbf{T}}_{\mathrm{R}}:=\left[\begin{array}{cc}{\mathbf{R}}_{\mathrm{p}}{\mathbf{V}}_{\mathrm{p},1}{\mathrm{\Sigma }}_{\mathrm{p},1}^{-12}& {\mathbf{R}}_{\mathrm{i}}{\mathbf{V}}_{\mathrm{i},1}{\mathrm{\Sigma }}_{\mathrm{i},1}^{-12}\end{array}\right]$

that balance the system, reduce the differential subsystem, and generate a minimal realization of the algebraic subsystem.

We obtain the reduced system (13) by setting

(15) $\begin{array}{cc}{\mathbf{E}}_{\mathrm{R}}:=& {\mathbf{W}}_{\mathrm{R}}^{\mathrm{T}}\mathbf{E}{\mathbf{T}}_{\mathrm{R}}:=\left[\begin{array}{cc}{I}_{{\mathit{r}}_{\mathrm{p}}}& 0\\ 0& {\mathbf{N}}_{\mathrm{R}}\end{array}\right],{\mathbf{A}}_{\mathrm{R}}:={\mathbf{W}}_{\mathrm{R}}^{\mathrm{T}}\mathbf{A}{\mathbf{T}}_{\mathrm{R}}:=\left[\begin{array}{cc}{\mathbf{J}}_{\mathrm{R}}& 0\\ 0& I_{{\mathit{r}}_{\mathrm{i}}}\end{array}\right],\\ {\mathbf{B}}_{\mathrm{R}}:=& \left[\begin{array}{c}{\mathbf{B}}_{\mathrm{R},1}\\ {\mathbf{B}}_{\mathrm{R},2}\end{array}\right]:={\mathbf{W}}_{\mathrm{R}}^{\mathrm{T}}\mathbf{B},{\mathbf{C}}_{\mathrm{R}}:=\left[\begin{array}{cc}{\mathbf{C}}_{\mathrm{R},1}& {\mathbf{C}}_{\mathrm{R},2}\end{array}\right]:=\mathbf{C}{\mathbf{T}}_{\mathrm{R}}.\end{array}$(15)

The generated reduced system has decoupled differential and algebraic states.

If ${\sigma }_{{\mathit{r}}_{\mathrm{p}}}>{\sigma }_{{\mathit{r}}_{\mathrm{p}}+1}$, then the ROM is asymptotically stable and one can estimate the output error of the reduced system by

(16) $\parallel \mathbf{y}-{\mathbf{y}}_{\mathrm{R}}{\parallel }_{L^2([0,\mathrm{\infty }),{\mathbb{R}}^p)}\le \parallel \mathcal{G}-{\mathcal{G}}_{\mathrm{R}}{\parallel }_{\mathrm{\infty }}\parallel \mathbf{u}{\parallel }_{L^2([0,\mathrm{\infty }),{\mathbb{R}}^m)}\le \left(2\sum _{\mathit{j}={\mathit{r}}_{\mathrm{p}}+1}^{\mathit{n}}{\sigma }_j\right)\parallel \mathbf{u}{\parallel }_{L^2([0,\mathrm{\infty }),{\mathbb{R}}^m)},$(16)

where $\mathcal{G}(\mathit{s}):=\mathbf{C}{\left(\mathit{s}\mathbf{E}-\mathbf{A}\right)}^{-1}\mathbf{B}\in \mathcal{R}{\mathcal{H}}_{\mathrm{\infty }}^{p,m}$ and ${\mathcal{G}}_{\mathrm{R}}(\mathit{s}):={\mathbf{C}}_{\mathrm{R}}{\left(\mathit{s}{\mathbf{E}}_{\mathrm{R}}-{\mathbf{A}}_{\mathrm{R}}\right)}^{-1}{\mathbf{B}}_{\mathrm{R}}\in \mathcal{R}{\mathcal{H}}_{\mathrm{\infty }}^{p,m}$ are the transfer functions of the original and the reduced system [Antoulas 2005, Theorem 7.9]. The ${\mathcal{H}}_{\mathrm{\infty }}$-norm is defined as $\parallel \mathcal{G}{\parallel }_{\mathrm{\infty }}:={sup}_{\mathit{\omega }\in \mathbb{R}}{\sigma }_{max}(\mathcal{G}(\mathrm{i}\mathit{\omega }))$ where ${\sigma }_{max}(\cdot )$ denotes the maximum singular value of its matrix argument.

Remark 1. In the balancing procedure outlined above, it is important that the number of inputs and outputs is very small compared to the state-space dimension, i.e., $\mathit{m},\mathit{p}\ll \mathit{N}$, which we assume throughout this work. In this case, the Gramians will typically have low numerical rank (Sorenson and Zhou 2002), which ensures that a low-rank approximation of the Gramians is possible and numerically efficient and that the FOM can be well-approximated by a ROM of low order. Balanced truncation model reduction is also possible in the case that one of the dimensions $\mathit{m}$ or $\mathit{p}$ is large while the other one is small, see, e.g. (Benner and Schneider 2013; Pulch, Narayan, and Stykel 2021), However, in this case, the numerical procedure is based on the approximation of an integral by quadrature. Our method cannot be generalized to this more general situation in a straightforward manner, but we believe that it would be an interesting topic for future research.

3.2.2. Solving projected Lyapunov equations

The aim of this subsection is to present numerical techniques to solve the projected Lyapunov equations (10(9) $\mathbf{E}{\mathbf{P}}_{\mathrm{p}}{\mathbf{A}}^{\mathrm{T}}+\mathbf{A}{\mathbf{P}}_{\mathrm{p}}{\mathbf{E}}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}},{\mathbf{P}}_{\mathrm{p}}={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{p}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(9) ) and (12(12) ${\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}\mathbf{A}-{\mathbf{E}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}\mathbf{E}=(I_N-{\mathbf{\Pi }}_{\mathrm{r}}{)}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}(I_N-{\mathbf{\Pi }}_{\mathrm{r}}),0={{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{i}}{\mathbf{\Pi }}_{\mathrm{l}}.$(12) ) in order to approximate the Gramians of system (1(1a) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{E}(\mathit{\mu })\mathbf{z}(\mathit{t})=\mathbf{A}(\mathit{\mu })\mathbf{z}(\mathit{t})+\mathbf{B}(\mathit{\mu })\mathbf{u}(\mathit{t}),$(1a) ). For standard systems with $\mathbf{E}=I_N$, there are several methods to solve the corresponding Lyapunov equations. If small-scale matrices are considered, the Bartels-Stewart algorithm (Bartels and Stewart 1972) or Hammarling’s method (Hammarling 1982) are used. These methods, however, are inefficient in the case of large matrices. In typical applications of practical relevance, the solution of a Lyapunov equation is often of very low numerical rank. So it is desired to compute low-rank factors of these solutions directly. State-of-the-art methods are the ADI method (Kürschner 2016; Penzl 2000), the sign function method (Benner et al. 2017) or Krylov subspace methods (Simoncini and Druskin 2009). In the literature, there are also several extensions to projected Lyapunov equations such as (Heinkenschloss, Sorensen, and Sun 2008; Stykel 2008).

In this section, we utilize the ADI method to solve the projected continuous-time Lyapunov equation (10(9) $\mathbf{E}{\mathbf{P}}_{\mathrm{p}}{\mathbf{A}}^{\mathrm{T}}+\mathbf{A}{\mathbf{P}}_{\mathrm{p}}{\mathbf{E}}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}},{\mathbf{P}}_{\mathrm{p}}={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{p}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(9) ) and the generalized Smith method to solve the discrete-time Lyapunov equation (12(11) $\mathbf{A}{\mathbf{P}}_{\mathrm{i}}{\mathbf{A}}^{\mathrm{T}}-\mathbf{E}{\mathbf{P}}_{\mathrm{i}}{\mathbf{E}}^{\mathrm{T}}=(I_N-{\mathbf{\Pi }}_{\mathrm{l}})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}(I_N-{\mathbf{\Pi }}_{\mathrm{l}}{)}^{\mathrm{T}},0={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{i}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(11) ). Here, we follow the ideas from (Stykel 2008).

We can use the following lemma from (Stykel 2008) to derive an equation that is equivalent to the projected continuous-time Lyapunov equation (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ).

Lemma 3.1. Let the matrix pencil $\mathit{s}\mathbf{E}-\mathbf{A}$ with $\mathbf{E},\mathbf{A}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$ be regular and asymptotically stable. Let further the matrix $\mathbf{A}$ be nonsingular and $\mathbf{B}\in {\mathbb{R}}^{\mathit{N},\mathit{m}}$. Assume that the left and right spectral projectors onto the finite spectrum of $\mathit{s}\mathbf{E}-\mathbf{A}$ are denoted by ${\mathbf{\Pi }}_{\mathrm{l}},{\mathbf{\Pi }}_{\mathrm{r}}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$. If $\mathit{p}\in \mathbb{C}$ with $\mathrm{Re}(\mathit{p})<0$ is not an eigenvalue of the pencil $\mathit{s}\mathbf{A}-\mathbf{E}$, then the projected discrete-time Lyapunov equation

(17) ${\mathbf{P}}_{\mathrm{p}}=\mathbf{S}(\mathit{p})\mathbf{R}(\mathit{p}){\mathbf{P}}_{\mathrm{p}}\mathbf{R}(\mathit{p}{)}^{\mathrm{H}}\mathbf{S}(\mathit{p}{)}^{\mathrm{H}}-2\mathrm{Re}(\mathit{p})\mathbf{S}(\mathit{p}){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}\mathbf{S}(\mathit{p}{)}^{\mathrm{H}},{\mathbf{P}}_{\mathrm{p}}={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{p}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}}$(17)

with

(18) $\mathbf{S}(\mathit{p}):=(\mathbf{E}+\mathit{p}\mathbf{A}{)}^{-1}\mathrm{and}\mathbf{R}(\mathit{p}):=\mathbf{E}-\bar{p}\mathbf{A}$(18)

is equivalent to the projected continuous-time Lyapunov Equation (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ), i.e. their solution sets coincide.

Lemma 3.1 motivates the ADI iteration

(19) $\begin{array}{cc}{\mathbf{P}}_0& :=0,\\ {\mathbf{P}}_k& :=\mathbf{S}(p_k)\mathbf{R}(p_k){\mathbf{P}}_{\mathit{k}-1}\mathbf{R}{(p_k)}^{\mathrm{H}}\mathbf{S}{(p_k)}^{\mathrm{H}}-2\mathrm{Re}(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}\mathbf{S}{(p_k)}^{\mathrm{H}},\mathit{k}=1,2,\dots .\end{array}$(19)

As shown in (Stykel 2008) the following proposition provides the convergence of the iteration (19).

Proposition 3.2. Let the matrix pencil $\mathit{s}\mathbf{E}-\mathbf{A}$ with $\mathbf{E},\mathbf{A}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$ be regular and asymptotically stable (i.e., all its finite eigenvalues have negative real part) and let $\mathbf{B}\in {\mathbb{R}}^{\mathit{N},\mathit{m}}$. Assume that the left and right spectral projectors onto the finite spectrum of $\mathit{s}\mathbf{E}-\mathbf{A}$ are denoted by ${\mathbf{\Pi }}_{\mathrm{l}},{\mathbf{\Pi }}_{\mathrm{r}}\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$. Suppose that we have given a sequence of shift parameters $(p_k{)}_{\mathit{k}\ge 0}$ with $\mathrm{Re}(p_k)<0$ for all $\mathit{k}\ge 0$ and with $p_{\mathit{k}+\mathrm{\ell }}=p_k$ for some $\mathrm{\ell }\ge 1$ and all $\mathit{k}=0,1,2,\dots$. Then the iteration (19) converges to the solution ${\mathbf{P}}_{\mathrm{p}}$ of the projected Lyapunov equation (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ).

Remark 2. The periodicity of the shift parameters ensures that they fulfill a non-Blaschke condition that is sufficient for the convergence of the ADI iteration (Massoudi, Opmeer, and Reis 2017).

To work with the ADI iteration more efficiently, our aim is to compute low-rank factors of ${\mathbf{P}}_{\mathrm{p}}$, i.e., we aim to compute a tall and skinny matrix $\mathbf{Z}$ such that ${\mathbf{P}}_{\mathrm{p}}\approx \mathbf{Z}{\mathbf{Z}}^{\mathrm{H}}$. We can represent the iteration (20) by the low-rank factors of ${\mathbf{P}}_k={\mathbf{Z}}_k{\mathbf{Z}}_k^{\mathrm{H}}$ with

${\mathbf{Z}}_k=\left[\begin{array}{cc}\mathit{\kappa }(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}& \mathbf{S}(p_k)\mathbf{R}(p_k){\mathbf{Z}}_{\mathit{k}-1}\end{array}\right]=[\mathit{\kappa }(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}\mathit{\kappa }(p_{\mathit{k}-1})\mathbf{S}(p_k)\mathbf{R}(p_k)\mathbf{S}(p_{\mathit{k}-1}){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}\text{break}\dots \mathit{\kappa }(p_1)\mathbf{S}(p_k)\mathbf{R}(p_k)\cdot \dots \cdot \mathbf{S}(p_2)\mathbf{R}(p_2)\mathbf{S}(p_1){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}],$

where $\mathit{\kappa }(p_k)=\sqrt{-\mathrm{Re}(p_k)}$ and ${\mathbf{Z}}_0$ is chosen as the empty matrix in ${\mathbb{R}}^{\mathit{N},0}$. We note that the following properties hold for all $\mathit{j},\mathit{k}=0,1,2,\dots$:

(20) $\begin{array}{cc}\mathit{\hspace{0.17em}}& \mathbf{S}(p_k)\mathbf{AS}(p_j)=\mathbf{S}(p_j)\mathbf{AS}(p_k),\mathbf{R}(p_k){\mathbf{A}}^{-1}\mathbf{R}(p_j)=\mathbf{R}(p_j){\mathbf{A}}^{-1}\mathbf{R}(p_k),\\ \mathit{\hspace{0.17em}}& \mathbf{S}(p_k)\mathbf{R}(p_j)={\mathbf{A}}^{-1}\mathbf{R}(p_j)\mathbf{S}(p_k)\mathbf{A}.\end{array}$(20)

We further define

${\mathbf{B}}_0:=\mathit{\kappa }(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}\mathrm{and}{\mathbf{F}}_j:=\frac{\kappa (p_j)}{\mathit{\kappa }(p_{\mathit{j}-1})}\mathbf{S}(p_j)\mathbf{R}(p_{\mathit{j}+1}),\mathit{j}=1,\dots ,\mathit{k}.$

Using (20), we obtain

${\mathbf{Z}}_k=\left[\begin{array}{cccc}{\mathbf{B}}_0& {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_0& \dots & {\mathbf{F}}_1\cdot \dots \cdot {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_0\end{array}\right].$

It remains to solve the discrete-time Lyapunov equation (11(10) ${\mathbf{E}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{A}+{\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{E}=-{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}{\mathbf{\Pi }}_{\mathrm{r}},{\mathbf{Q}}_{\mathrm{p}}={{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}{\mathbf{\Pi }}_{\mathrm{l}}$(10) ). Under the assumption that $\mathbf{A}$ is nonsingular, (11) is equivalent to the transformed discrete-time Lyapunov equation

${\mathbf{P}}_{\mathrm{i}}-{\mathbf{A}}^{-1}\mathbf{E}{\mathbf{P}}_{\mathrm{i}}{\mathbf{E}}^{\mathrm{T}}{\mathbf{A}}^{-\mathrm{T}}={\mathbf{A}}^{-1}(I_N-{\mathbf{\Pi }}_{\mathrm{l}})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}(I_N-{\mathbf{\Pi }}_{\mathrm{l}}{)}^{\mathrm{T}}{\mathbf{A}}^{-\mathrm{T}},0={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{i}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}}.$

This equation could be solved with the Smith method (Stykel 2008). Since ${\mathbf{A}}^{-1}(I_N-{\mathbf{\Pi }}_{\mathrm{l}})=(I_N-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}$ and the matrix $(I_N-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}\mathbf{E}={\mathbf{A}}^{-1}\mathbf{E}(I_N-{\mathbf{\Pi }}_{\mathrm{r}})$ is nilpotent with the nilpotency index $\mathit{\nu }$ of the matrix $\mathbf{N}$ of the quasi-Weierstrass form of the pencil $\mathit{s}\mathbf{E}-\mathbf{A}$, the iteration stops after $\mathit{\nu }$ steps. This leads to the unique solution

${\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}({\mathbf{A}}^{-1}\mathbf{E}{)}^k(I_N-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{A}}^{-\mathrm{T}}(I_N-{\mathbf{\Pi }}_{\mathrm{r}}{)}^{\mathrm{T}}(({\mathbf{A}}^{-1}\mathbf{E}{)}^{\mathrm{T}}{)}^k.$

We note that we have to solve multiple linear systems with the matrix $\mathbf{A}$ to obtain the solution ${\mathbf{P}}_{\mathrm{i}}$. Therefore, we can utilize the sparse LU decomposition to do this efficiently. Instead of computing the full matrix ${\mathbf{P}}_{\mathrm{i}}$ we can also generate the low-rank factors ${\mathbf{P}}_{\mathrm{i}}=\mathbf{Y}{\mathbf{Y}}^{\mathrm{T}}$ as

$\mathbf{Y}=\left[\begin{array}{cccc}(\mathit{I}-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}\mathbf{B}& {\mathbf{A}}^{-1}\mathbf{E}(\mathit{I}-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}\mathbf{B}& \dots & {({\mathbf{A}}^{-1}\mathbf{E})}^{\mathit{\nu }-1}(\mathit{I}-{\mathbf{\Pi }}_{\mathrm{r}}){\mathbf{A}}^{-1}\mathbf{B}\end{array}\right].$

Remark 3. In many relevant applications, the solution ${\mathbf{P}}_{\mathrm{i}}$ can be explicitly calculated by making use of specific block structures in $\mathit{s}\mathbf{E}-\mathbf{A}$. Moreover, the improper Hankel singular values are often zero (namely, when the FOM’s transfer function is strictly proper) in which case we do not have to solve the improper Lyapunov equations at all, see also (Stykel 2006).

3.2.3. Alternative approach to determine the algebraic parts

Since the reduced basis method presented later in this paper is applied only to the differential part of the system, in this section we discuss ways to construct reduced representations of the algebraic part. One way is to use the reduction provided by the improper Gramians as explained in Subsection 3.2.1. This is a suitable choice for parameter-independent problems. However, it is difficult to build a reduced representation for parameter-dependent problems, if the improper part depends on parameters. Then, the discrete-time Lyapunov equations have to be solved for each parameter value of interest in the online phase which can be too expensive. For this reason, in this subsection, we present an alternative approach that is also viable in the parametric setting and allows an efficient evaluation in the online phase (see also Subsection 4.3).

Our method is inspired by (Schwerdtner et al. 2023) in which the polynomial part of the transfer functions of DAE systems of index at most two is approximated, but the methodology can also be extended to systems of even higher index. Consider the transfer function $\mathcal{G}(\mathit{s})$ of system (1) (for a fixed parameter) which can be decomposed into its strictly proper part and a polynomial part

$\mathcal{G}(\mathit{s})={\mathcal{G}}_{\mathrm{sp}}(\mathit{s})+{\mathcal{G}}_{\mathrm{poly}}(\mathit{s}),$

where the polynomial part is given as

${\mathcal{G}}_{\mathrm{poly}}(\mathit{s})=-\sum _{\mathit{k}=0}^{\mathit{\nu }-1}\mathbf{C}{\mathbf{T}}^{-1}\left[\begin{array}{cc}0& 0\\ 0& {\mathbf{N}}^k\end{array}\right]{\mathbf{W}}^{-1}\mathbf{B}{s}^k,$

with $\mathbf{T}$, $\mathbf{W}$, and $\mathbf{N}$ with index of nilpotency $\mathit{\nu }$ as in (8), see (Stykel 2004) for the details. Hence, for ${\mathbf{W}}^{-1}\mathbf{B}=:\left[\begin{array}{c}{\mathbf{B}}_1\\ {\mathbf{B}}_2\end{array}\right]$ and $\mathbf{C}{\mathbf{T}}^{-1}=:\left[\begin{array}{cc}{\mathbf{C}}_1& {\mathbf{C}}_2\end{array}\right]$ partitioned conforming to the partioning of (8), the polynomial part can be rewritten as

${\mathcal{G}}_{\mathrm{poly}}(\mathit{s})=-\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{C}}_2{\mathbf{N}}^k{\mathbf{B}}_2{s}^k,$

where ${\mathbf{M}}_k=-{\mathbf{C}}_2{\mathbf{N}}^k{\mathbf{B}}_2$, $\mathit{k}=0,\dots ,\mathit{\nu }-1$ denote the $\mathit{k}$-th Markov parameters of the transfer function. As shown in (Antoulas, Gosea, and Heinkenschloss 2020), the strictly proper part ${\mathcal{G}}_{\mathrm{sp}}(\mathrm{i}\mathit{\omega })$ converges to zero for $\mathit{\omega }\to \mathrm{\infty }$. For this reason we can approximate the Markov parameters if we consider sufficiently large values of $\mathit{\omega }\in \mathbb{R}$. We exemplify the methodology for a few cases.

The index one case In this case, we simply have

${\mathbf{M}}_0\approx \mathcal{G}(\mathrm{i}\mathit{\omega })$

for a sufficiently large value of $\mathit{\omega }\in \mathbb{R}$ and moreover, ${\mathbf{M}}_k=0$ for $\mathit{k}\ge 1$.

The index two case In this case, ${\mathbf{M}}_k=0$ for $\mathit{k}\ge 2$ and hence ${\mathcal{G}}_{\mathrm{poly}}(\mathit{s})={\mathbf{M}}_0+{\mathbf{M}}_1\mathit{s}$. By inserting $\mathrm{i}\mathit{\omega }\in \mathrm{i}\mathbb{R}$ we obtain that

$\mathrm{Re}\left({\mathcal{G}}_{\mathrm{poly}}(\mathrm{i}\mathit{\omega })\right)={\mathbf{M}}_0,\mathrm{Im}\left({\mathcal{G}}_{\mathrm{poly}}(\mathrm{i}\mathit{\omega })\right)=\mathit{\omega }{\mathbf{M}}_1.$

Hence, for sufficiently large values $\mathit{\omega }\in \mathbb{R}$ we can approximate

${\mathbf{M}}_0\approx \mathrm{Re}\left(\mathcal{G}(\mathrm{i}\mathit{\omega })\right)\mathrm{and}{\mathbf{M}}_1\approx \frac{\mathrm{Im}\left(\mathcal{G}(\mathrm{i}\omega )\right)}{\omega }.$

The index three case In the index three case, we approximate the polynomial part of the form ${\mathcal{G}}_{\mathrm{poly}}(\mathit{s})={\mathbf{M}}_0+{\mathbf{M}}_1\mathit{s}+{\mathbf{M}}_2{s}^2.$ Inserting $\mathrm{i}\mathit{\omega }\in \mathrm{i}\mathbb{R}$ into the polynomial part of the transfer function ${\mathcal{G}}_{\mathrm{poly}}$ provides a real and an imaginary part

$\mathrm{Re}\left({\mathcal{G}}_{\mathrm{poly}}(\mathrm{i}\mathit{\omega })\right)={\mathbf{M}}_0-{\omega }^2{\mathbf{M}}_2,\mathrm{Im}\left({\mathcal{G}}_{\mathrm{poly}}(\mathrm{i}\mathit{\omega })\right)=\mathit{\omega }{\mathbf{M}}_1.$

Then the Markov parameters are approximated by

${\mathbf{M}}_2\approx \frac{\mathrm{Re}\left(\mathcal{G}(\mathrm{i}{\omega }_1)-\mathcal{G}(\mathrm{i}{\omega }_2)\right)}{{\omega }_2^2-{\omega }_1^2},{\mathbf{M}}_1\approx \frac{\mathrm{Im}\left(\mathcal{G}(\mathrm{i}\omega )\right)}{\omega },{\mathbf{M}}_0\approx \mathrm{Re}\left(\mathcal{G}(\mathrm{i}\mathit{\omega })\right)+{\omega }^2{\mathbf{M}}_2$

for sufficiently large $\mathit{\omega }$, ${\omega }_1$, and ${\omega }_2$ with ${\omega }_1\ne {\omega }_2$.

System realization from Markov parameters If the polynomial part ${\mathcal{G}}_{\mathrm{poly}}(\mathit{s})={\sum }_{k=0}^{\nu -1}{\mathbf{M}}_k{s}^k$ has been extracted and all Markov parameters ${\mathbf{M}}_k$, $\mathit{k}=0,\dots ,\mathit{\nu }-1$ have been determined we can realize it by a DAE system as follows. Assume that ${\mathbf{M}}_k={\mathbf{U}}_k{\mathrm{\Sigma }}_k{\mathbf{V}}_k^{\mathrm{T}}$ are compact SVDs of ${\mathbf{M}}_k$ with positive definite ${\mathrm{\Sigma }}_k\in {\mathbb{R}}^{r_k,r_k}$ for all $\mathit{k}=1,\dots ,\mathit{\nu }-1$. Then for $\mathit{k}=1,\dots ,\mathit{\nu }-1$ we define the matrices

${\hat{\mathbf{B}}}_{\mathit{k}}:=\left[\begin{array}{c}0\\ ⋮\\ 0\\ {\mathbf{V}}_k^{\mathrm{T}}\end{array}\right]\in {\mathbb{R}}^{k\cdot r_k,\mathit{m}},{\hat{\mathbf{C}}}_{\mathit{k}}:=\left[\begin{array}{cccc}-{\mathbf{U}}_k& 0& \cdots & 0\end{array}\right]\in {\mathbb{R}}^{p,k\cdot r_k}$

and

${\hat{\mathbf{E}}}_{\mathit{k}}:=\left[\begin{array}{cccc}0& {\mathrm{\Sigma }}_k^{\frac{1}{k}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \ddots & \ddots & \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \ddots & {\mathrm{\Sigma }}_k^{\frac{1}{k}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& 0\end{array}\right]\in {\mathbb{R}}^{k\cdot r_k,\mathit{k}\cdot r_k},\mathrm{sothat}\hat{\mathbf{E}}\mathit{kk}=\left[\begin{array}{cccc}0& \cdots & 0& {\mathrm{\Sigma }}_k\\ \mathit{\hspace{0.17em}}& \ddots & \mathit{\hspace{0.17em}}& 0\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \ddots & ⋮\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& 0\end{array}\right].$

This yields

${\hat{\mathbf{C}}}_{\mathit{k}}{\left(\mathit{s}{\hat{\mathbf{E}}}_{\mathit{k}}-I_{k\cdot r_k}\right)}^{-1}{\hat{\mathbf{B}}}_{\mathit{k}}=-\sum _{\mathit{j}=0}^k{\hat{\mathbf{C}}}_{\mathit{k}}\hat{\mathbf{E}}\mathit{kj}{\hat{\mathbf{B}}}_{\mathit{k}}{s}^j=-{\hat{\mathbf{C}}}_{\mathit{k}}\hat{\mathbf{E}}\mathit{kk}{\hat{\mathbf{B}}}_{\mathit{k}}{s}^k={\mathbf{M}}_k{s}^k.$

Now we obtain an overall realization of the algebraic part of the system as

$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{ccc}{\hat{\mathbf{E}}}_1& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \ddots & \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& {\hat{\mathbf{E}}}_{\mathit{\nu }-1}\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_1(\mathit{t})\\ ⋮\\ {\mathbf{z}}_{\mathit{\nu }-1}(\mathit{t})\end{array}\right]=\left[\begin{array}{c}{\mathbf{z}}_1(\mathit{t})\\ ⋮\\ {\mathbf{z}}_{\mathit{\nu }-1}(\mathit{t})\end{array}\right]+\left[\begin{array}{c}{\hat{\mathbf{B}}}_1\\ ⋮\\ {\hat{\mathbf{B}}}_{\mathit{k}}\end{array}\right]\mathbf{u}(\mathit{t}),$

$\hat{\mathbf{y}}(\mathit{t})=\left[\begin{array}{ccc}{\hat{\mathbf{C}}}_1& \dots & {\hat{\mathbf{C}}}_{\mathit{k}}\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_1(\mathit{t})\\ ⋮\\ {\mathbf{z}}_{\mathit{\nu }-1}(\mathit{t})\end{array}\right]+{\mathbf{M}}_0\mathbf{u}(\mathit{t}).$

4. Reduced basis method

In this section, we return to the problem of reducing parameter-dependent systems. Since the solution of the Lyapunov equations is the most expensive part of balanced truncation, we aim to limit the number of Lyapunov equation solves. To achieve this, the reduced basis method presented by (Son and Stykel 2017) can be applied. Besides the regularity and asymptotic stability of the matrix pencil $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$ for all $\mathit{\mu }\in \mathcal{D}$ we impose the following assumption on the problem (1) for the rest of the paper.

Assumption 4.1. We assume that $\mathcal{D}\subset {\mathbb{R}}^d$ is a bounded and connected domain. We assume further that the matrix-valued functions $\mathbf{E}(\cdot )$, $\mathbf{A}(\cdot )$, $\mathbf{B}(\cdot )$, $\mathbf{C}(\cdot )$ in (1) as well as the left and right spectral projector functions ${\mathbf{\Pi }}_{\mathrm{l}}(\cdot )$ and ${\mathbf{\Pi }}_{\mathrm{r}}(\cdot )$ on the deflating subspace associated with the finite eigenvalues of $(\mathit{s}\mathbf{E}-\mathbf{A})(\cdot )$ depend continuously on the parameter $\mathit{\mu }$.

Assumption 4.1 ensures that the solutions of the parameter-dependent projected Lyapunov equations also depend continuously on $\mathit{\mu }\in \mathcal{D}$. Then, also the ROM depends continuously on $\mathit{\mu }$ for a fixed reduced order.

In practice, many examples are continuously parameter-dependent, e.g., in electrical circuits, the resistances can be continuously dependent on the temperature. Moreover, inductances and capacitances may be assumed continuous if one wishes to tune them within a parameter optimization. In other examples, for instance, mechanical systems, the parameters may represent (continuously varying) material constants. On the other hand, there are also examples where the dependence on the parameters is not continuous. This is, for example, the case for electrical circuits in which the network topology changes by switching. These discontinuities are often known since, in practice, the user has insight into the considered problem. In this case, the parameter domain could be decomposed into several subdomains, on which the reduction can be performed separately.

In the following, we will focus on the reduction of the differential part of the DAE system (1). The main idea consists of finding a reduced representation of the solutions of the Lyapunov equations (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ) of the form

(21) ${\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\approx \mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathrm{with}\mathbf{Z}(\mathit{\mu })=\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu })\mathrm{\forall }\mathit{\mu }\in \mathcal{D},$(21)

where $\mathbf{V}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},{\mathit{r}}_{\mathbf{V}}}$ has orthonormal columns with ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\mathbf{V}(\mathit{\mu })=\mathbf{V}(\mathit{\mu })$ and where ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$ denotes the projector onto the right deflating subspace of $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$ for the parameter value $\mathit{\mu }$. Under the assumption that $\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{E}(\mathit{\mu })\mathbf{V}(\mathit{\mu })$ is invertible and the pencil $\mathit{s}\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{E}(\mathit{\mu })\mathbf{V}(\mathit{\mu })-\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{A}(\mathit{\mu })\mathbf{V}(\mathit{\mu })$ is regular and asymptotically stable (This condition can be guaranteed a priori by the strict dissipativity assumptions in Subsection 4.2), the low-rank factor $\tilde{\mathbf{Z}}(\mathit{\mu })$ solves the reduced Lyapunov equation

$\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{E}(\mathit{\mu })\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{V}(\mathit{\mu })$

(22) $+\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{A}(\mathit{\mu })\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{V}(\mathit{\mu })=-\mathbf{V}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}\mathbf{V}(\mathit{\mu }).$(22)

The equation ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\mathbf{V}(\mathit{\mu })=\mathbf{V}(\mathit{\mu })$ replaces the condition

${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })(\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }))(\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }){)}^{\mathrm{H}}{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}=(\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }))(\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu }){)}^{\mathrm{H}}.$

Algorithm 1 Reduced basis method for projected Lyapunov equations (offline phase)

Input: Dynamical system (1), test parameter set ${\mathcal{D}}_{\mathrm{Test}}$, tolerance $\mathrm{Tol}$.

Output: Reduced basis matrix ${\mathbf{V}}_{\mathrm{glob}}$

1: Choose any ${\mu }_1\in {\mathcal{D}}_{\mathrm{Test}}$.

2: Solve the projected Lyapunov equation (24(23) $\mathbf{V}(\mathit{\mu })=\mathrm{orth}({\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}}),$(23) ) for $\mathit{\mu }={\mu }_1$ to obtain $\mathbf{Z}({\mu }_1)$.

3: Set $\mathcal{M}:=\{{\mu }_1\}$.

4: Set ${\mathbf{V}}_{\mathrm{glob}}:=\mathrm{orth}(\mathbf{Z}({\mu }_1))$.

5: Set $\hat{\mu }:=\arg \mathrm{ma}{\mathrm{x}}_{\mathit{\mu }\in {\mathcal{D}}_{\mathrm{Test}}\mathrm{\setminus }\mathcal{M}}{\mathrm{\Delta }}_{\mathbf{V}}(\mathit{\mu })$.

6: Set ${\mathrm{\Delta }}_{\mathbf{V}}^{max}:={\mathrm{\Delta }}_{\mathbf{V}}(\hat{\mu })$.

7: while ${\mathrm{\Delta }}_{\mathbf{V}}^{max}>\mathrm{Tol}$ do

8: Solve the projected Lyapunov equation (24(23) $\mathbf{V}(\mathit{\mu })=\mathrm{orth}({\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}}),$(23) ) for $\mathit{\mu }=\hat{\mu }$ to obtain $\mathbf{Z}(\hat{\mu })$.

9: Set $\mathcal{M}:=\mathcal{M}\cup \{\hat{\mu }\}$.

10: Set ${\mathbf{V}}_{\mathrm{glob}}:=\mathrm{orth}([{\mathbf{V}}_{\mathrm{glob}},\mathbf{Z}(\hat{\mu })])$.

11: Set $\hat{\mu }:=\arg \mathrm{ma}{\mathrm{x}}_{\mathit{\mu }\in {\mathcal{D}}_{\mathrm{Test}}\mathrm{\setminus }\mathcal{M}}{\mathrm{\Delta }}_{\mathbf{V}}(\mathit{\mu })$.

12: Set ${\mathrm{\Delta }}_{\mathbf{V}}^{max}:={\mathrm{\Delta }}_{\mathbf{V}}(\hat{\mu })$.

13: end while

In practice, we aim to determine a matrix ${\mathbf{V}}_{\mathrm{glob}}\in {\mathbb{R}}^{\mathit{N},{\mathit{r}}_{\mathrm{glob}}}$ with $\mathit{N}\gg r_{\mathrm{glob}}$ that contains the information of $\mathbf{V}(\mathit{\mu })$ for all $\mathit{\mu }\in \mathcal{D}$ globally. When we have determined such a ${\mathbf{V}}_{\mathrm{glob}}$, we set

(23) $\mathbf{V}(\mathit{\mu })=\mathrm{orth}({\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}}),$(23)

where $\mathrm{orth}({\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}})$ denotes the orthonormalized columns of ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}}$. If the matrix ${\mathbf{V}}_{\mathrm{glob}}$ has been found, then the ROM for a particular parameter $\mathit{\mu }\in \mathcal{D}$ can be determined very efficiently in the so-called online phase, where one simply solves (22) with (21) and (23) to determine the two Gramians and then projects the system as in Subsection 3.2.1.

The computation of the reduced basis ${\mathbf{V}}_{\mathrm{glob}}$ is done beforehand in the offline phase. There, we solve the full-order Lyapunov equation only for those parameters, where their solutions are currently approximated worst to enrich the current reduced basis. A posteriori error estimators are employed on a test parameter set to find these points efficiently.

We determine ${\mathbf{V}}_{\mathrm{glob}}$ by considering a test parameter set ${\mathcal{D}}_{\mathrm{Test}}\subset \mathcal{D}$. We begin by computing a low-rank factor $\mathbf{Z}({\mu }_1)$ in ${\mu }_1\in {\mathcal{D}}_{\mathrm{Test}}$ such that ${\mathbf{P}}_{\mathrm{p}}({\mu }_1)=\mathbf{Z}({\mu }_1)\mathbf{Z}({\mu }_1{)}^{\mathrm{H}}$ solves the projected Lyapunov equation

(24) $\begin{array}{rl}& \mathbf{E}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{A}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}},\\ & {\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}={\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\end{array}$(24)

for $\mathit{\mu }={\mu }_1$.

The first reduced basis is then given by ${\mathbf{V}}_{\mathrm{glob}}:=\mathrm{orth}(\mathbf{Z}({\mu }_1))$. Next, we determine the test parameter ${\mu }_2\in {\mathcal{D}}_{\mathrm{Test}}$ for which the solution ${\mathbf{P}}_{\mathrm{p}}({\mu }_2)$ of the Lyapunov equation (24(24) $\begin{array}{rl}& \mathbf{E}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{A}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}},\\ & {\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}={\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\end{array}$(24) ) is approximated worst by using (22), (23), and (24). For that, one of the error estimators ${\mathrm{\Delta }}_{\mathbf{V}}(\mathit{\mu })$ presented in Subsection 4.1 is utilized. For this parameter ${\mu }_2$, we solve the projected Lyapunov equation (24(24) $\begin{array}{rl}& \mathbf{E}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{A}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}},\\ & {\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}={\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\end{array}$(24) ) with $\mathit{\mu }={\mu }_2$ to obtain the low-rank factor $\mathbf{Z}({\mu }_2)$. Then, we define the new reduced basis by setting ${\mathbf{V}}_{\mathrm{glob}}:=\mathrm{orth}([{\mathbf{V}}_{\mathrm{glob}},\mathbf{Z}({\mu }_2)])$. This procedure is repeated until the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}$ is smaller than a prescribed tolerance for every test parameter in ${\mathcal{D}}_{\mathrm{Test}}$. This method results in Algorithm 1.

Remark 4. Under the assumption that the projected Lyapunov equations (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ), (10(10) ${\mathbf{E}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{A}+{\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}\mathbf{E}=-{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}}{\mathbf{C}}^{\mathrm{T}}\mathbf{C}{\mathbf{\Pi }}_{\mathrm{r}},{\mathbf{Q}}_{\mathrm{p}}={{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}}{\mathbf{Q}}_{\mathrm{p}}{\mathbf{\Pi }}_{\mathrm{l}}$(10) ) are solved exactly for all parameters $\mathit{\mu }\in \mathcal{D}$ using the reduced Lyapunov equations of the form (22), the ${\mathcal{H}}_{\mathrm{\infty }}$ error bound (16) is also valid in the entire parameter domain. Unfortunately, this bound is no longer valid if numerical approximation errors affect the numerically computed Gramians which is already the case in the parameter-independent case. A general and in-depth analysis of the effect of such errors on the Hankel singular values is still an open problem. On the other hand, the numerically computed Hankel singular values can still be used to obtain an approximate error bound.

4.1. Error estimation

We want to estimate the error $\mathfrak{E}(\mathit{\mu }):={\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })-\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}$ with the help of the residual

$\mathrm{\Re }(\mathit{\mu }):=\mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}+{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}.$

Lemma 4.2. Let $\mathit{\mu }\in \mathcal{D}$ be given. Let the matrix pencil $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$ with $\mathbf{E}(\mathit{\mu }),\mathbf{A}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$ be regular and asymptotically stable. Assume that the left and right spectral projectors onto the finite spectrum of $\mathit{s}\mathbf{E}(\mathit{\mu })-\mathbf{A}(\mathit{\mu })$ are denoted by ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }),{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$. Let the controllability Gramian ${\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},\mathit{N}}$ be defined as in (8) and let it be approximated as described in Equation (21)(21) ${\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\approx \mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathrm{with}\mathbf{Z}(\mathit{\mu })=\mathbf{V}(\mathit{\mu })\tilde{\mathbf{Z}}(\mathit{\mu })\mathrm{\forall }\mathit{\mu }\in \mathcal{D},$(21) . Then the error $\mathfrak{E}(\mathit{\mu }):=P_{\text{p}}\left(\mathit{\mu }\right)-\mathit{Z}\left(\mathit{\mu }\right)\mathit{Z}{(\mathit{\mu })}^{\text{H}}$ solves the projected error equation

(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25)

Proof. The condition

$\mathfrak{E}\left(\mathit{\mu }\right)={\Pi }_{\text{r}}\left(\mathit{\mu }\right)\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}$

follows from ${\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })={\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}$ and $\mathbf{Z}(\mathit{\mu })={\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })$. For the remaining equation, we consider its left-hand side

$\mathbf{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}=\mathbf{A}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{E}(\mathit{\mu }){\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}-\mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}-\mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}-\mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}-\mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}.$

With ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{A}(\mathit{\mu })=\mathbf{A}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$, ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{E}(\mathit{\mu })=\mathbf{E}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$ and $\mathbf{Z}(\mathit{\mu })={\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })$ we obtain

$\mathbf{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}-{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}$

which proves the statement.

Similarly to (Son and Stykel 2017), we consider the linear system

$\mathcal{L}\left(\mathit{\mu }\right)\mathbf{x}\left(\mathit{\mu }\right)=\mathbf{b}\left(\mathit{\mu }\right),\mathbf{x}\left(\mathit{\mu }\right)={\mathbf{\Pi }}_{\mathrm{r}}^{\otimes }\left(\mathit{\mu }\right)\mathbf{x}\left(\mathit{\mu }\right)$

with

(26) $\mathit{ℒ}(\mathit{\mu })\text{ }:=-\mathit{A}(\mathit{\mu })\otimes \mathit{E}(\mathit{\mu })-\mathit{E}(\mathit{\mu })\otimes \mathit{A}(\mathit{\mu }),\mathit{x}(\mathit{\mu })\text{ }:=\text{vec}\left(P_{\text{p}}(\mathit{\mu })\right),\mathit{b}(\mathit{\mu }):=\text{vec}\left({\Pi }_{\text{l}}(\mathit{\mu })\mathit{B}{B}^{\text{T}}{\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}}\right),{\Pi }_{\text{r}}^{\otimes }\left(\mathit{\mu }\right):={\Pi }_{\text{r}}(\mathit{\mu })\otimes {\Pi }_{\text{r}}(\mathit{\mu }),$(26)

where the operator $\otimes$ denotes the Kronecker product and $\mathrm{vec}$ the vectorization operator that stacks the columns of the matrix on top of one another. This linear system is equivalent to the projected Lyapunov equation (24(23) $\mathbf{V}(\mathit{\mu })=\mathrm{orth}({\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }){\mathbf{V}}_{\mathrm{glob}}),$(23) ). We note that $\mathcal{L}(\mathit{\mu })$ is a singular matrix and that the unique solvability is enforced by the additional restriction given by the projection with ${\Pi }_{\text{r}}^{\otimes }\left(\mathit{\mu }\right)$.

Using the linear system representation and the abbreviations $\mathbf{e}(\mathit{\mu }):=\mathrm{vec}(\mathfrak{E}(\mathit{\mu }))$, $\mathbf{r}(\mathit{\mu }):=\mathrm{vec}(\mathrm{\Re }(\mathit{\mu }))$, we can rewrite the error Equation (25)(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25) as

$\mathcal{L}\left(\mathit{\mu }\right)\mathbf{e}\left(\mathit{\mu }\right)=\mathbf{r}\left(\mathit{\mu }\right),\mathbf{e}\left(\mathit{\mu }\right)={\mathbf{\Pi }}_{\mathrm{r}}^{\otimes }\left(\mathit{\mu }\right)\mathbf{e}\left(\mathit{\mu }\right).$

This equation is equivalent to

(27) ${\mathbf{\Pi }}_{\mathrm{l}}^{\otimes }\left(\mathit{\mu }\right)\mathcal{L}\left(\mathit{\mu }\right){\mathbf{\Pi }}_{\mathrm{r}}^{\otimes }\left(\mathit{\mu }\right)\mathbf{e}\left(\mathit{\mu }\right)={\mathbf{\Pi }}_{\mathrm{l}}^{\otimes }\left(\mathit{\mu }\right)\mathbf{r}\left(\mathit{\mu }\right)=\mathbf{r}\left(\mathit{\mu }\right),\mathrm{where}{\mathbf{\Pi }}_{\mathrm{l}}^{\otimes }\left(\mathit{\mu }\right):={\mathbf{\Pi }}_{\mathrm{l}}\left(\mathit{\mu }\right)\otimes {\mathbf{\Pi }}_{\mathrm{l}}\left(\mathit{\mu }\right).$(27)

Due to the unique solvability of the latter, we obtain

$\mathit{e}(\mathit{\mu })={\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{ℒ}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })\mathit{r}(\mathit{\mu }),$

where $(\cdot {)}^{†}$ denotes the Moore-Penrose inverse. Then, we can estimate

(28) $\parallel \text{}\mathfrak{E}(\mathit{\mu }){\parallel }_{\text{F}}=\parallel \text{}\mathit{e}(\mathit{\mu })\text{}\parallel 2\le \parallel \text{ }{\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{L}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu }){\parallel }_2\parallel \mathit{r}(\mathit{\mu })\text{}\parallel 2=1{\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{L}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })\parallel \text{}\mathit{\Re }(\mathit{\mu })\text{}\parallel \text{F}\le \frac{1}{\mathit{\alpha }(\mathit{\mu })}\parallel \text{}\mathit{\Re }(\mathit{\mu })\text{}\parallel \text{F}=:{\Delta }_V^{(1)}(\mathit{\mu }).$(28)

Here, ${\sigma }_{min}^{+}(\cdot )$ denotes the smallest positive singular value of its matrix argument and $\mathit{\alpha }(\mathit{\mu })$ is an easily computable lower bound of ${\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{I}}^{\otimes }\left(\mathit{\mu }\right)\mathit{ℒ}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }\left(\mathit{\mu }\right)$. The bound ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}(\mathit{\mu })$ is the first possible error estimator that provides a rather conservative bound in practice.

To improve the estimator (28), we consider the error bound presented in (Schmidt, Wittwar, and Haasdonk 2020). Assume that $\widehat{\mathfrak{E}}(\mathit{\mu })\approx \mathfrak{E}(\mathit{\mu })$ is an error estimate that is not necessarily an upper bound. An error bound ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}(\mathit{\mu })$ is derived by

(29) $\begin{array}{cc}\parallel \mathfrak{E}(\mathit{\mu }){\parallel }_{\mathrm{F}}& =\parallel \mathfrak{E}(\mathit{\mu })+\hat{\mathfrak{E}}(\mathit{\mu })-\hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}\le \parallel \hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}+\parallel \mathfrak{E}(\mathit{\mu })-\hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}\\ \\ \mathit{\hspace{0.17em}}& \le \parallel \hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}+\frac{1}{\mathit{\alpha }(\mathit{\mu })}\parallel \hat{\mathrm{\Re }}(\mathit{\mu }){\parallel }_{\mathrm{F}}:={\mathrm{\Delta }}_{\mathbf{V}}^{(2)}(\mathit{\mu }),\end{array}$(29)

where the residual $\hat{\mathrm{\Re }}(\mathit{\mu })$ is defined as

$\hat{\mathrm{\Re }}(\mathit{\mu }):=\mathbf{A}(\mathit{\mu })\left(\mathbf{Z}(\mathit{\mu })\mathbf{Z}{(\mathit{\mu })}^{\mathrm{H}}+\hat{\mathfrak{E}}(\mathit{\mu })\right)\mathbf{E}(\mathit{\mu }{)}^{\mathrm{T}}+\mathbf{E}(\mathit{\mu })\left(\mathbf{Z}(\mathit{\mu })\mathbf{Z}{(\mathit{\mu })}^{\mathrm{H}}+\hat{\mathfrak{E}}(\mathit{\mu })\right)\mathbf{A}(\mathit{\mu }{)}^{\mathrm{T}}\text{break}+{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\mathbf{B}(\mathit{\mu })\mathbf{B}(\mathit{\mu }{)}^{\mathrm{T}}{\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu }{)}^{\mathrm{T}}.$

It remains to obtain an error estimate $\hat{\mathfrak{E}}(\mathit{\mu })$ for all $\mathit{\mu }\in {\mathcal{D}}_{\mathrm{Test}}$. To do so, we solve the error Equation (26)(26) $\mathit{ℒ}(\mathit{\mu })\text{ }:=-\mathit{A}(\mathit{\mu })\otimes \mathit{E}(\mathit{\mu })-\mathit{E}(\mathit{\mu })\otimes \mathit{A}(\mathit{\mu }),\mathit{x}(\mathit{\mu })\text{ }:=\text{vec}\left(P_{\text{p}}(\mathit{\mu })\right),\mathit{b}(\mathit{\mu }):=\text{vec}\left({\Pi }_{\text{l}}(\mathit{\mu })\mathit{B}{B}^{\text{T}}{\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}}\right),{\Pi }_{\text{r}}^{\otimes }\left(\mathit{\mu }\right):={\Pi }_{\text{r}}(\mathit{\mu })\otimes {\Pi }_{\text{r}}(\mathit{\mu }),$(26) approximately by modifying the projected ADI iteration from Section 3.2.2 and the reduced basis method in Algorithm 1.

The right-hand side of the error Equation (25)(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25) can be written as the product ${\mathbf{B}}_{\mathrm{l}}(\mathit{\mu }){\mathbf{B}}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{H}}$, where

${\mathbf{B}}_{\mathrm{l}}(\mathit{\mu }):={\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{ccc}\mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })& \mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })& \mathbf{B}(\mathit{\mu })\end{array}\right],$

${\mathbf{B}}_{\mathrm{r}}(\mathit{\mu }):={\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{ccc}\mathbf{E}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })& \mathbf{A}(\mathit{\mu })\mathbf{Z}(\mathit{\mu })& \mathbf{B}(\mathit{\mu })\end{array}\right].$

For simplicity of notation, we consider the modified ADI iteration for the parameter-independent case such that $\mathfrak{E}(\mathit{\mu })\equiv$ $\mathfrak{E},B_{\text{l}}(\mathit{\mu })\equiv B_{\text{l}}$, ${\mathbf{B}}_{\mathrm{r}}(\mathit{\mu })\equiv {\mathbf{B}}_{\mathrm{r}}$, and ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\equiv {\mathbf{\Pi }}_{\mathrm{l}}$. Analogously to the derivation in Section 3.2.2 (and under the same conditions), it can be shown that the iteration

${\mathfrak{E}}_k={\mathbf{Z}}_{\mathrm{l},k}({\mathbf{Z}}_{\mathrm{r},k}{)}^{\mathrm{H}}=\mathbf{S}(p_k)\mathbf{R}(p_k){\mathfrak{E}}_{\mathit{k}-1}\mathbf{R}(p_k{)}^{\mathrm{H}}\mathbf{S}(p_k{)}^{\mathrm{H}}-2\mathrm{Re}(p_k)\mathbf{S}(p_k){\mathbf{B}}_{\mathrm{l}}{\mathbf{B}}_{\mathrm{r}}^{\mathrm{T}}\mathbf{S}(p_k{)}^{\mathrm{H}}=\mathbf{S}(p_k)\mathbf{R}(p_k){\mathbf{Z}}_{\mathrm{l},k-1}({\mathbf{Z}}_{\mathrm{r},k-1}{)}^{\mathrm{H}}\mathbf{R}(p_k{)}^{\mathrm{H}}\mathbf{S}(p_k{)}^{\mathrm{H}}-2\mathrm{Re}(p_k)\mathbf{S}(p_k){\mathbf{B}}_{\mathrm{l}}{\mathbf{B}}_{\mathrm{r}}^{\mathrm{T}}\mathbf{S}(p_k{)}^{\mathrm{H}}$

converges to the solution $\mathfrak{E}$ of (25), where $\mathbf{S}(p_k)$ and $\mathbf{R}(p_k)$ are defined as in (18) for a shift parameter $p_k\in \mathbb{C}$ with $\mathrm{Re}(p_k)<0$. The factors ${\mathbf{Z}}_{\mathrm{l},k}$ and ${\mathbf{Z}}_{\mathrm{r},k}$ can be written as

${\mathbf{Z}}_{\mathrm{l},k}=\left[\begin{array}{cccc}{\mathbf{B}}_{\mathrm{l},0}& {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_{\mathrm{l},0}& \dots & {\mathbf{F}}_1\cdot \dots \cdot {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_{\mathrm{l},0}\end{array}\right],\text{break}{\mathbf{Z}}_{\mathrm{r},k}=\left[\begin{array}{cccc}{\mathbf{B}}_{\mathrm{r},0}& {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_{\mathrm{r},0}& \dots & {\mathbf{F}}_1\cdot \dots \cdot {\mathbf{F}}_{\mathit{k}-1}{\mathbf{B}}_{\mathrm{r},0}\end{array}\right],$

where

${\mathbf{B}}_{\mathrm{l},0}:=\mathit{\kappa }(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}{\mathbf{B}}_{\mathrm{l}},{\mathbf{B}}_{\mathrm{r},0}:=\mathit{\kappa }(p_k)\mathbf{S}(p_k){\mathbf{\Pi }}_{\mathrm{l}}{\mathbf{B}}_{\mathrm{r}}\mathrm{and}{\mathbf{F}}_j:=\frac{\kappa (p_j)}{\mathit{\kappa }(p_{\mathit{j}-1})}\mathbf{S}(p_j)\mathbf{R}(p_{\mathit{j}+1})$

for $\mathit{j}=1,\dots ,\mathit{k}$ and with $\mathit{\kappa }(p_j)=\sqrt{-2\mathrm{Re}(p_j)}$.

We include the modified projected ADI iteration into the reduced basis method presented in Algorithm 1 by solving the error Equation (25)(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25) in Step 2 and 8 to obtain the two factors ${\mathbf{Z}}_{\mathrm{l}}(\mathit{\mu }),{\mathbf{Z}}_{\mathrm{r}}(\mathit{\mu })$ with $Z_{\text{l}}(\mathit{\mu })Z_{\text{r}}{(\mathit{\mu })}^{\text{H}}\approx \mathfrak{E}(\mathit{\mu })$. The orthonormal basis computation in Step 4 and 4 is replaced by ${\mathbf{V}}_{\mathrm{glob}}:=\mathrm{orth}([{\mathbf{V}}_{\mathrm{glob}},{\mathbf{Z}}_{\mathrm{l}}(\mathit{\mu })])$. Note that here, the columns of ${\mathbf{Z}}_{\mathrm{l}}(\mathit{\mu })$ and ${\mathbf{Z}}_{\mathrm{r}}(\mathit{\mu })$ span the same subspace. As error estimator, we use ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}$. After we have determined the orthonormal basis, we solve Equation (25)(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25) on the corresponding subspace to get an approximate error $\widehat{\mathfrak{E}}(\mathit{\mu })$, that we use for the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}(\mathit{\mu })$.

Remark 5. If the factor $\mathbf{Z}(\mathit{\mu })\in {\mathbb{R}}^{\mathit{N},{\mathit{n}}_{\mathbf{Z}}}$ is already of large dimension $n_{\mathbf{Z}}$, the right-hand side of the error Equation (25)(25) $\mathit{A}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{\text{T}}+\mathit{E}(\mathit{\mu })\mathfrak{E}(\mathit{\mu })\mathit{A}{(\mathit{\mu })}^{\text{T}}=-\mathit{\Re }(\mathit{\mu })=-{\Pi }_{\text{l}}(\mathit{\mu })\mathit{\Re }(\mathit{\mu }){\Pi }_{\text{l}}{(\mathit{\mu })}^{\text{T}},\mathfrak{E}(\mathit{\mu })={\Pi }_{\text{r}}(\mathit{\mu })\mathfrak{E}(\mathit{\mu }){\Pi }_{\text{r}}{(\mathit{\mu })}^{\text{T}}.$(25) described by ${\mathbf{B}}_{\mathrm{l}}(\mathit{\mu })$ and ${\mathbf{B}}_{\mathrm{r}}(\mathit{\mu })$ is of even larger rank $2n_{\mathbf{Z}}+\mathit{m}$, and hence, the solution $\mathfrak{E}(\mathit{\mu })$ cannot in general be well-represented by low-rank factors. This means that the error space basis might be large which leads to a time consuming error estimation.

4.2. Strictly dissipative systems

For the model examples from Section 2 we can make use of their structure to derive the bound $\mathit{\alpha }(\mathit{\mu })$ in a computationally advantageous way. Therefore, we will impose an additional assumption on the matrix pencils $\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })$ in (4) and (6) that is formed from submatrices of $\mathbf{E}(\mathit{\mu })$ and $\mathbf{A}(\mathit{\mu })$.

Definition 4.3. The set $\{\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })\in \mathbb{R}[\mathit{s}{]}^{\mathit{n},\mathit{n}}|\mathit{\mu }\in \mathcal{D}\}$ of matrix pencils is called uniformly strictly dissipative, if

$\mathit{E}(\mathit{\mu })=\mathit{E}(\mathit{\mu }{)}^{\mathrm{T}}>0\mathrm{and}{A}^S(\mathit{\mu }):=\frac{1}{2}(\mathit{A}(\mathit{\mu })+\mathit{A}(\mathit{\mu }{)}^{\mathrm{T}})<0\mathrm{\forall }\mathit{\mu }\in \mathcal{D}.$

We will consider this condition for our motivating examples from Section 2 and determine a corresponding lower bound $\mathit{\alpha }(\mathit{\mu })$.

4.2.1. Stokes-like systems

We consider the Stokes-like system (4). Here, the pencil set $\{\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })|\mathit{\mu }\in \mathcal{D}\}$ is naturally uniformly strictly dissipative.

As shown in (Stykel 2006) the projection matrices are given as

(30) ${\text{⁋ia}}_{\mathrm{l}}(\mathit{\mu })={\text{⁋ia}}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}=\left[\begin{array}{cc}\text{⁋ia}(\mathit{\mu })& -\text{⁋ia}(\mathit{\mu })\mathit{A}(\mathit{\mu })\mathit{E}{(\mathit{\mu })}^{-1}\mathit{G}(\mathit{\mu })(\mathit{G}{(\mathit{\mu })}^{\mathrm{T}}\mathit{E}{(\mathit{\mu })}^{-1}\mathit{G}(\mathit{\mu }{))}^{-1}\\ 0& 0\end{array}\right]$(30)

where

(31) $\text{⁋ia}(\mathit{\mu })=I_n-\mathit{G}(\mathit{\mu })(\mathit{G}(\mathit{\mu }{)}^{\mathrm{T}}\mathit{E}(\mathit{\mu }{)}^{-1}\mathit{G}(\mathit{\mu }){)}^{-1}\mathit{G}(\mathit{\mu }{)}^{\mathrm{T}}\mathit{E}(\mathit{\mu }{)}^{-1}.$(31)

Even though for the Stokes-like system in (4) we have $A^S(\mathit{\mu })=\mathit{A}(\mathit{\mu })<0$ for all $\mathit{\mu }\in \mathcal{D}$, in this section we will already treat the more general case of a uniformly strictly dissipative set of matrix pencils $\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })$ with nonsymmetric $\mathit{A}(\mathit{\mu })$ as this case will arise later in Subsection 4.2.2.

We continue finding a value for $\mathit{\alpha }(\mathit{\mu })$.

Theorem 4.4. Let $\mathit{\mu }\in \mathcal{D}$ be given and consider the Stokes-like system (4). Let $A^S(\mathit{\mu })$, $\mathcal{L}(\mathit{\mu })$, ${\mathbf{\Pi }}_{\mathrm{r}}^{\otimes }\left(\mathit{\mu }\right)$, ${\mathbf{\Pi }}_{\mathrm{I}}^{\otimes }\left(\mathit{\mu }\right)$ be defined as in Definition 4.3, (26), (27) with ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })={\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu }{)}^{\mathrm{T}}$ as in (30), respectively. Then the bound

${\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{I}}^{\otimes }(\mathit{\mu })\mathit{ℒ}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })\ge 2{\lambda }_{\mathit{min}}(\mathit{E}(\mathit{\mu })){\lambda }_{\mathit{min}}(-A^S(\mathit{\mu }))$

is satisfied, where ${\lambda }_{min}(\cdot )$ denotes the minimum eigenvalue of its (symmetric) matrix arguments.

Proof. Let $\mathrm{\Pi }(\mathit{\mu })$ be as in (31). Making use of the condition $\mathrm{\Pi }(\mathit{\mu })\mathit{G}(\mathit{\mu })=0$, we can deduce

${\mathbf{\Pi }}_{\mathrm{I}}^{\otimes }(\mathit{\mu })\mathcal{L}(\mathit{\mu }){\mathbf{\Pi }}_{\mathrm{r}}^{\otimes }(\mathit{\mu })=-\left({\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{cc}\mathit{A}(\mathit{\mu })& \mathit{G}(\mathit{\mu })\\ \mathit{G}{(\mathit{\mu })}^{\mathrm{T}}& 0\end{array}\right]{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\right)\otimes \left({\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{cc}\mathit{E}(\mathit{\mu })& 0\\ 0& 0\end{array}\right]{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\right)-\left({\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{cc}\mathit{E}(\mathit{\mu })& 0\\ 0& 0\end{array}\right]{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\right)\otimes \left({\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })\left[\begin{array}{cc}\mathit{A}(\mathit{\mu })& \mathit{G}(\mathit{\mu })\\ \mathit{G}{(\mathit{\mu })}^{\mathrm{T}}& 0\end{array}\right]{\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })\right)=-\left[\begin{array}{cc}\mathrm{\Pi }(\mathit{\mu })\mathit{A}(\mathit{\mu })\mathrm{\Pi }{(\mathit{\mu })}^{\mathrm{T}}& 0\\ 0& 0\end{array}\right]\otimes \left[\begin{array}{cc}\mathrm{\Pi }(\mathit{\mu })\mathit{E}(\mathit{\mu })\mathrm{\Pi }{(\mathit{\mu })}^{\mathrm{T}}& 0\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}\mathrm{\Pi }(\mathit{\mu })\mathit{E}(\mathit{\mu })\mathrm{\Pi }{(\mathit{\mu })}^{\mathrm{T}}& 0\\ 0& 0\end{array}\right]\otimes \left[\begin{array}{cc}\mathrm{\Pi }(\mathit{\mu })\mathit{A}(\mathit{\mu })\mathrm{\Pi }{(\mathit{\mu })}^{\mathrm{T}}& 0\\ 0& 0\end{array}\right]$

Define ${\mathrm{\Pi }}^{\otimes }(\mathit{\mu }):=\mathrm{\Pi }(\mathit{\mu })\otimes \mathrm{\Pi }(\mathit{\mu })$ and the matrix

$\mathit{L}(\mathit{\mu }):=-\mathit{A}(\mathit{\mu })\otimes \mathit{E}(\mathit{\mu })-\mathit{E}(\mathit{\mu })\otimes \mathit{A}(\mathit{\mu }).$

Denoting the smallest positive eigenvalue of a symmetric matrix by ${\lambda }_{min}^{+}(\cdot )$, Corollary 3.1.5 from (Horn and Johnson 1991) and the eigenvalue properties of the Kronecker product, given in (Horn and Johnson 1991; Lancaster and Tismenetsky 1985) (see also (Son and Stykel 2017) as well as a Rayleigh quotient argument lead to the estimate

$\text{}{\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{I}}^{\otimes }(\mathit{\mu })\mathit{ℒ}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })={\sigma }_{{\mathit{min}}^{+}}\left({\Pi }^{\otimes }(\mathit{\mu })\mathit{L}(\mathit{\mu }){\Pi }^{\otimes }{(\mathit{\mu })}^{\text{T}}\right)$

$\ge {\lambda }_{min}^{+}\left(\frac{1}{2}{\mathrm{\Pi }}^{\otimes }(\mathit{\mu })\left(\mathit{L}(\mathit{\mu })+\mathit{L}{(\mathit{\mu })}^{\mathrm{T}}\right){\mathrm{\Pi }}^{\otimes }{(\mathit{\mu })}^{\mathrm{T}}\right)$

$=\underset{\mathit{v}\in \mathrm{im}{\mathrm{\Pi }}^{\otimes }{(\mathit{\mu })}^{\mathrm{T}},\parallel \mathit{v}{\parallel }_2=1}{min}\frac{1}{2}{v}^{\mathrm{T}}{\mathrm{\Pi }}^{\otimes }(\mathit{\mu })\left(\mathit{L}(\mathit{\mu })+\mathit{L}{(\mathit{\mu })}^{\mathrm{T}}\right){\mathrm{\Pi }}^{\otimes }(\mathit{\mu }{)}^{\mathrm{T}}\mathit{v}$

$=\underset{v\in \text{im\hspace{0.05em}}{\Pi }^{\otimes }{(\mathit{\mu })}^{\text{T}},1\mathit{pt}\parallel \mathit{v}{\parallel }_2=1}{\mathit{min}}\frac{1}{2}{v}^{\text{T}}\left(\mathit{L}(\mathit{\mu })+\mathit{L}{(\mathit{\mu })}^{\text{T}}\right)\mathit{v}$

$\ge \underset{\mathit{v}\in {\mathbb{R}}^{n^2},\parallel \mathit{v}{\parallel }_2=1}{min}\frac{1}{2}{v}^{\mathrm{T}}\left(\mathit{L}(\mathit{\mu })+\mathit{L}{(\mathit{\mu })}^{\mathrm{T}}\right)\mathit{v}$

$={\lambda }_{min}\left(\frac{1}{2}\left(\mathit{L}(\mathit{\mu })+\mathit{L}{(\mathit{\mu })}^{\mathrm{T}}\right)\right)=2{\lambda }_{min}(\mathit{E}(\mathit{\mu })){\lambda }_{min}(-A^S(\mathit{\mu })).$

□

Since computing the eigenvalues of $\mathit{E}(\mathit{\mu })$ and $\mathit{A}(\mathit{\mu })$ for all requested parameters $\mathit{\mu }$ can be expensive, there are several ways to find cheap lower bounds of $2{\lambda }_{min}(\mathit{E}(\mathit{\mu })){\lambda }_{min}(-A^S(\mathit{\mu }))$, some of which are derived in (Son and Stykel 2017). Based on these bounds, we can, e.g., estimate

(32) $2{\lambda }_{min}(\mathit{E}(\mathit{\mu })){\lambda }_{min}(-A^S(\mathit{\mu }))\ge 2\underset{k=0,\dots ,n_E}{min}\frac{{\Theta }_k^E(\mathit{\mu })}{{\Theta }_k^E(\bar{\mu })}{\lambda }_{min}(\mathit{E}(\bar{\mu }))\underset{k=0,\dots ,n_A}{min}\frac{{\Theta }_k^A(\mathit{\mu })}{{\Theta }_k^A(\bar{\mu })}{\lambda }_{min}(-A^S(\bar{\mu }))\text{break}=:\mathit{\alpha }(\mathit{\mu }),$(32)

where $\bar{\mu }$ is an arbitrary-fixed parameter in $\mathcal{D}$ and where we make use of the convention that ${\mathrm{\Theta }}_0^E={\mathrm{\Theta }}_0^A\equiv 1$. Here, we compute the eigenvalues of $\mathit{E}(\bar{\mu })$ and $A^S(\bar{\mu })$ only once to find a lower bound $\mathit{\alpha }(\mathit{\mu })$ of ${\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{ℒ}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })$ for every parameter $\mathit{\mu }\in \mathcal{D}$.

4.2.2. Mechanical systems

we turn to the mechanical system from Subsection 2.2. There, the corresponding matrix pencil set

$\left\{\mathit{s}\left[\begin{array}{cc}{I}_{n_x}& 0\\ 0& \mathit{M}(\mathit{\mu })\end{array}\right]-\left[\begin{array}{cc}0& I_{n_x}\\ -\mathit{K}(\mathit{\mu })& -\mathit{D}(\mathit{\mu })\end{array}\right]\in \mathbb{R}{[\mathit{s}]}^{2n_x,2n_x}|\mathit{\mu }\in \mathcal{D}\right\}$

is not uniformly strictly dissipative. However, it can be made uniformly strictly dissipative by appropriate generalized state-space transformations as we will show now.

Here, we apply a transformation presented in (Eid et al. 2011; Panzer, Wolf, and Lohmann 2012). We observe that by adding productive zeros and for an auxiliary function $\mathit{\gamma }:\mathcal{D}\to (0,\mathrm{\infty })$ yet to be chosen, the system (5) is equivalent to the system

$\frac{\mathrm{d}}{\mathrm{d}t}\mathit{K}(\mathit{\mu })\mathit{x}(\mathit{t})+\mathit{\gamma }(\mathit{\mu })\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}\mathit{M}(\mathit{\mu })\mathit{x}(\mathit{t})=-\mathit{\gamma }(\mathit{\mu })\mathit{K}(\mathit{\mu })\mathit{x}(\mathit{t})+\frac{\mathrm{d}}{\mathrm{d}t}\mathit{K}(\mathit{\mu })\mathit{x}(\mathit{t})-\mathit{\gamma }(\mathit{\mu })\frac{\mathrm{d}}{\mathrm{d}t}\mathit{D}(\mathit{\mu })\mathit{x}(\mathit{t})+\mathit{\gamma }(\mathit{\mu })\mathit{G}(\mathit{\mu })\mathit{\lambda }(\mathit{t})+\mathit{\gamma }(\mathit{\mu })B_x\mathit{u}(\mathit{t}),\text{break}\mathit{\gamma }(\mathit{\mu })\frac{\mathrm{d}}{\mathrm{d}t}\mathit{M}(\mathit{\mu })\mathit{x}(\mathit{t})+\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}\mathit{M}(\mathit{\mu })\mathit{x}(\mathit{t})=-\mathit{K}(\mathit{\mu })\mathit{x}(\mathit{t})+\mathit{\gamma }(\mathit{\mu })\frac{\mathrm{d}}{\mathrm{d}t}\mathit{M}(\mathit{\mu })\mathit{x}(\mathit{t})\text{break}-\frac{\mathrm{d}}{\mathrm{d}t}\mathit{D}(\mathit{\mu })\mathit{x}(\mathit{t})+\mathit{G}(\mathit{\mu })\mathit{\lambda }(\mathit{t})+B_x(\mathit{\mu })\mathit{u}(\mathit{t}),\text{break}0=\mathit{\gamma }(\mathit{\mu })\mathit{G}(\mathit{\mu }{)}^{\mathrm{T}}\mathit{x}(\mathit{t}),\text{break}0=\frac{\mathrm{d}}{\mathrm{d}t}\mathit{G}(\mathit{\mu }{)}^{\mathrm{T}}\mathit{x}(\mathit{t}),$

where $\mathit{x}(\mathit{t})$ is equal to $x_1(\mathit{t})$ and $\frac{\mathrm{d}}{\mathrm{d}t}\mathit{x}(\mathit{t})$ equal to $x_2(\mathit{t})$ from system (5). The last equation is a direct consequence from the equation above and poses no extra restrictions. These equations can be written as first-order system. By defining the matrices

(33) $\mathit{E}(\mathit{\mu }):=\left[\begin{array}{cc}\mathit{K}(\mathit{\mu })& \mathit{\gamma }(\mathit{\mu })\mathit{M}(\mathit{\mu })\\ \mathit{\gamma }(\mathit{\mu })\mathit{M}(\mathit{\mu })& \mathit{M}(\mathit{\mu })\end{array}\right],\mathit{A}(\mathit{\mu }):=\left[\begin{array}{cc}-\mathit{\gamma }(\mathit{\mu })\mathit{K}(\mathit{\mu })& \mathit{K}(\mathit{\mu })-\mathit{\gamma }(\mathit{\mu })\mathit{D}(\mathit{\mu })\\ -\mathit{K}(\mathit{\mu })& -\mathit{D}(\mathit{\mu })+\mathit{\gamma }(\mathit{\mu })\mathit{M}(\mathit{\mu })\end{array}\right],$(33)

$\mathit{B}(\mathit{\mu }):=\left[\begin{array}{c}\mathit{\gamma }(\mathit{\mu })B_x(\mathit{\mu })\\ B_x(\mathit{\mu })\end{array}\right],$

and replacing $\mathit{G}(\mathit{\mu })$ by $\mathit{B}(\mathit{\mu }):=\left[\begin{array}{c}\mathit{\gamma }(\mathit{\mu })B_x(\mathit{\mu })\\ B_x(\mathit{\mu })\end{array}\right]$, we generate a system in the form (5) such that we can apply the methods from Subsection 4.2.1 for this kind of system in the following, if $\mathit{\gamma }(\mathit{\mu })$ is chosen properly. This is possible, because even if $\mathit{A}(\mathit{\mu })$ is unsymmetric, the projector $\mathrm{\Pi }(\mathit{\mu })$ in (32) will be the same in both left and right spectral projectors ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })$ and ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$ and hence, the estimate in Theorem 4.4 is still valid. Only an estimation of the form (32) is not possible, since the matrices $\mathit{E}(\mathit{\mu })$ and $A^S(\mathit{\mu })$ cannot be written as affine decompositions in which all summands are positive or negative semidefinite, respectively. So in our numerical experiments we make use of Theorem 4.4 directly.

Theorem 4.5. Consider the mechanical system (5). Then the matrix pencil set $\{\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })|\mathit{\mu }\in \mathcal{D}\}$ with $\mathit{E}(\mathit{\mu })$ and $\mathit{A}(\mathit{\mu })$ defined in (33) is uniformly strictly dissipative, if $0<\mathit{\gamma }(\mathit{\mu })<{\gamma }^{\ast }(\mathit{\mu })$ with

${\gamma }^{\ast }(\mathit{\mu })=\frac{{\mathit{\lambda }}_{min}(\mathit{D}(\mathit{\mu }))}{{\lambda }_{max}(\mathit{M}(\mathit{\mu }))+\frac{1}{4}{\lambda }_{max}{(\mathit{D}(\mathit{\mu }))}^2{\lambda }_{max}(\mathit{K}{(\mathit{\mu })}^{-1})}$

is satisfied for all $\mathit{\mu }\in \mathcal{D}.$Proof. As shown in (Panzer, Wolf, and Lohmann 2012), the pencil set $\{\mathit{sE}(\mathit{\mu })-\mathit{A}(\mathit{\mu })|\mathit{\mu }\in \mathcal{D}\}$ is uniformly strictly dissipative, if we choose $\mathit{\gamma }(\mathit{\mu })>0$ such that

$\mathit{\gamma }(\mathit{\mu })<{\lambda }_{min}\left(\mathit{D}(\mathit{\mu }){\left(M(\mu )+\frac{1}{4}\mathit{D}(\mathit{\mu })\mathit{K}{(\mathit{\mu })}^{-1}\mathit{D}(\mathit{\mu })\right)}^{-1}\right)\mathrm{\forall }\mathit{\mu }\in \mathcal{D}.$

Using the estimate ${\lambda }_{min}(\mathit{XY})\ge {\lambda }_{min}(\mathit{X}){\lambda }_{min}(\mathit{Y})$ for symmetric and positive definite matrices $\mathit{X}$ and $\mathit{Y}$ of conforming dimension (Wang, Zhang 1992, Theorem 4), we arrive at the choice for $\mathit{\gamma }(\mathit{\mu })$ given by

(34) $\mathit{\gamma }(\mathit{\mu }){\lambda }_{max}\left(\mathit{M}(\mathit{\mu })+\frac{1}{4}\mathit{D}(\mathit{\mu })\mathit{K}{(\mathit{\mu })}^{-1}\mathit{D}(\mathit{\mu })\right)<{\lambda }_{min}(\mathit{D}(\mathit{\mu }))\mathrm{\forall }\mathit{\mu }\in \mathcal{D}.$(34)

Because of the symmetry and positive definiteness of $\mathit{D}(\mathit{\mu })\mathit{K}(\mathit{\mu }{)}^{-1}\mathit{D}(\mathit{\mu })$ and the submultiplicativity of the norm we obtain

${\lambda }_{max}\left(\mathit{D}(\mathit{\mu })\mathit{K}{(\mathit{\mu })}^{-1}\mathit{D}(\mathit{\mu })\right)={∥\mathit{D}(\mathit{\mu })\mathit{K}{(\mathit{\mu })}^{-1}\mathit{D}(\mathit{\mu })∥}_2$

$\le \parallel \mathit{D}(\mathit{\mu }){\parallel }_2^2\cdot \parallel \mathit{K}(\mathit{\mu }{)}^{-1}{\parallel }_2={\lambda }_{max}(\mathit{D}(\mathit{\mu }){)}^2\cdot {\lambda }_{max}(\mathit{K}(\mathit{\mu }{)}^{-1}).$

By Weyl’s lemma (Horn and Johnson 1991; Lancaster and Tismenetsky 1985), we also have

${\lambda }_{max}\left(\mathit{M}(\mathit{\mu })+\frac{1}{4}\mathit{D}(\mathit{\mu })\mathit{K}{(\mathit{\mu })}^{-1}\mathit{D}(\mathit{\mu })\right)\le {\lambda }_{max}(\mathit{M}(\mathit{\mu }))+\frac{1}{4}{\lambda }_{max}(\mathit{D}(\mathit{\mu })\mathit{K}(\mathit{\mu }{)}^{-1}\mathit{D}(\mathit{\mu })).$

Therefore, condition (34) is satisfied, if we choose $\mathit{\gamma }(\mathit{\mu })$ such that

$\mathit{\gamma }(\mathit{\mu })\left({\lambda }_{max}(\mathit{M}(\mathit{\mu }))+\frac{1}{4}{\lambda }_{max}{(\mathit{D}(\mathit{\mu }))}^2{\lambda }_{max}(\mathit{K}{(\mathit{\mu })}^{-1})\right)<{\lambda }_{min}(\mathit{D}(\mathit{\mu }))\mathrm{\forall }\mathit{\mu }\in \mathcal{D}.$

This can be achieved by the choice

(35) $\mathit{\gamma }(\mathit{\mu })<\frac{{\mathit{\lambda }}_{min}(\mathit{D}(\mathit{\mu }))}{{\lambda }_{max}(\mathit{M}(\mathit{\mu }))+\frac{1}{4}{\lambda }_{max}{(\mathit{D}(\mathit{\mu }))}^2{\lambda }_{max}(\mathit{K}{(\mathit{\mu })}^{-1})}.$(35)

□

To avoid computing the eigenvalues of $\mathit{M}(\mathit{\mu }),\mathit{D}(\mathit{\mu }),\mathit{K}(\mathit{\mu })$ for every parameter $\mathit{\mu }$, we further use the estimates

${\lambda }_{max}(\mathit{D}(\mathit{\mu }))\le {\lambda }_{max}(D_0)+\sum _{\mathit{k}=1}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{max}(D_k)=\sum _{\mathit{k}=0}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{max}(D_k),$

${\lambda }_{min}(\mathit{D}(\mathit{\mu }))\ge {\lambda }_{min}(D_0)+\sum _{\mathit{k}=1}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{min}(D_k)=\sum _{\mathit{k}=0}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{min}(D_k),$

that are a consequence of Weyl’s lemma. Similar estimates are also valid for $\mathit{M}(\mathit{\mu })$ and $\mathit{K}(\mathit{\mu })$. Then we can further estimate $\mathit{\gamma }(\mathit{\mu })$ as

(36) $\mathit{\gamma }(\mathit{\mu })<{\gamma }^{\ast }(\mathit{\mu })=\frac{{\sum }_{\mathit{k}=0}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{min}(D_k)}{{\sum }_{\mathit{k}=0}^{n_M}{\Theta }_k^M(\mathit{\mu }){\lambda }_{max}(M_k)+\frac{1}{4}{\left({\sum }_{\mathit{k}=0}^{n_D}{\Theta }_k^D(\mathit{\mu }){\lambda }_{max}(D_k)\right)}^2/\left({\sum }_{\mathit{k}=0}^{n_K}{\Theta }_k^K(\mathit{\mu }){\lambda }_{min}(K_k)\right)},$(36)

while ensuring that we can choose $\mathit{\gamma }(\mathit{\mu })>0$ for all $\mathit{\mu }\in \mathcal{D}$.

The benefit of the above estimate is that the extremal eigenvalues of the matrices $M_k$, $D_k$, and $K_k$ have to be computed only once in the beginning and otherwise, its evaluation is cheap, since ${\Theta }_k^M(\cdot )$, ${\Theta }_k^D(\cdot )$, ${\Theta }_k^K(\cdot )$ are scalar-valued functions.

4.3. Treatment of the algebraic parts

The goal of this subsection is to discuss the reduction of the algebraic part of the system which is potentially also parameter-dependent. To do so, we first precompute specific matrices in the offline phase which will later be used to efficiently construct the Markov parameters in the online phase for a particular parameter configuration and then determine a reduced representation of the algebraic subsystem using the construction from the last paragraph in Subsection 3.2.3. In order for this construction to be efficient, we assume that the parameter influence has an additional low-rank structure, i.e., with regard to the affine decomposition (3), we assume that ${\mathbf{E}}_k$, ${\mathbf{A}}_k$, ${\mathbf{B}}_k$, and ${\mathbf{C}}_k$ are low-rank matrices for $\mathit{k}\ge 1$, i.e., they admit low-rank factorizations

(37) $\begin{array}{cc}{\mathbf{E}}_k& ={\mathbf{U}}_{\mathbf{E},k}{\mathbf{V}}_{\mathbf{E},k}^{\mathrm{T}},\mathit{k}=1,\dots ,n_{\mathbf{E}},\\ {\mathbf{A}}_k& ={\mathbf{U}}_{\mathbf{A},k}{\mathbf{V}}_{\mathbf{A},k}^{\mathrm{T}},\mathit{k}=1,\dots ,n_{\mathbf{A}}.\end{array}$(37)

To determine the Markov parameters, for a fixed set of sample points ${\omega }_i\in \mathbb{R}$, $\mathit{i}=1,\dots ,n_{\mathrm{S}}$, we aim to evaluate $\mathcal{G}(\mathit{\mu },\mathrm{i}{\omega }_i)$ efficiently for all needed parameters $\mathit{\mu }\in \mathcal{D}$. In other words, we aim for an efficient procedure to evaluate the transfer function $\mathcal{G}(\mathit{\mu },\mathrm{i}\mathit{\omega })$ for fixed value $\mathit{\omega }\in \mathbb{R}$ and varying $\mathit{\mu }$. Due to (38) and since $n_{\mathbf{E}}$ and $n_{\mathbf{A}}$ are assumed to be very small, also the matrix $\mathrm{i}\mathit{\omega }{\sum }_{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{E}}}{\mathrm{\Theta }}_k^{\mathbf{E}}(\mathit{\mu }){\mathbf{E}}_k-{\sum }_{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{A}}}{\mathrm{\Theta }}_k^{\mathbf{A}}(\mathit{\mu }){\mathbf{A}}_k$ is of low rank and we can factorize it as

$\mathrm{i}\mathit{\omega }\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{E}}}{\mathrm{\Theta }}_k^{\mathbf{E}}(\mathit{\mu }){\mathbf{E}}_k-\sum _{\mathit{k}=1}^{{\mathit{n}}_{\mathbf{A}}}{\mathrm{\Theta }}_k^{\mathbf{A}}(\mathit{\mu }){\mathbf{A}}_k$

$={\left[\begin{array}{c}{\mathbf{U}}_{\mathbf{E},1}^{\mathrm{T}}\\ ⋮\\ {\mathbf{U}}_{\mathbf{E},{\mathit{n}}_{\mathbf{E}}}^{\mathrm{T}}\\ {\mathbf{U}}_{\mathbf{A},1}^{\mathrm{T}}\\ ⋮\\ {\mathbf{U}}_{\mathbf{A},{\mathit{n}}_{\mathbf{A}}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cccccc}\mathrm{i}\mathit{\omega }{\mathrm{\Theta }}_1^{\mathbf{E}}(\mathit{\mu })& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \ddots & \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathrm{i}\mathit{\omega }{\mathrm{\Theta }}_{{\mathit{n}}_{\mathbf{E}}}^{\mathbf{E}}(\mathit{\mu })& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& -{\mathrm{\Theta }}_1^{\mathbf{A}}(\mathit{\mu })& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \ddots & \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& -{\mathrm{\Theta }}_{{\mathit{n}}_{\mathbf{A}}}^{\mathbf{A}}(\mathit{\mu })\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathbf{E},1}^{\mathrm{T}}\\ ⋮\\ {\mathbf{V}}_{\mathbf{E},{\mathit{n}}_{\mathbf{E}}}^{\mathrm{T}}\\ {\mathbf{V}}_{\mathbf{A},1}^{\mathrm{T}}\\ ⋮\\ {\mathbf{V}}_{\mathbf{A},{\mathit{n}}_{\mathbf{A}}}^{\mathrm{T}}\end{array}\right]=:\mathbf{U\Omega }(\mathit{\mu },\mathrm{i}\mathit{\omega }){\mathbf{V}}^{\mathrm{T}}.$

It holds that

$\mathcal{G}(\mathit{\mu },\mathrm{i}\mathit{\omega })=\sum _{\mathit{i}=0}^{{\mathit{n}}_{\mathbf{C}}}\sum _{\mathit{j}=0}^{{\mathit{n}}_{\mathbf{B}}}{\mathrm{\Theta }}_i^{\mathbf{C}}(\mathit{\mu }){\mathrm{\Theta }}_j^{\mathbf{B}}(\mathit{\mu }){\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0+\mathbf{U\Omega }(\mathit{\mu },\mathrm{i}\mathit{\omega }){\mathbf{V}}^{\mathrm{T}}\right)}^{-1}{\mathbf{B}}_j.$

With the help of the Sherman-Morrison-Woodbury identity, we obtain

$\begin{array}{rl}& {\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0+\mathbf{U\Omega }(\mathit{\mu },\mathrm{i}\mathit{\omega }){\mathbf{V}}^{\mathrm{T}}\right)}^{-1}{\mathbf{B}}_j\\ & ={\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}{\mathbf{B}}_j-{\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}\mathbf{U}{\left(\mathbf{\Omega }{(\mathit{\mu },\mathrm{i}\omega )}^{-1}\text{break}+{\mathbf{V}}^{\mathrm{T}}{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}\mathbf{U}\right)}^{-1}{\mathbf{V}}^{\mathrm{T}}{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}{\mathbf{B}}_j.\end{array}$

The matrices ${\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}{\mathbf{B}}_j$, ${\mathbf{C}}_i{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}\mathbf{U}$, ${\mathbf{V}}^{\mathrm{T}}{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}\mathbf{U}$, and ${\mathbf{V}}^{\mathrm{T}}{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}{\mathbf{B}}_j$ are all of small dimension and they can be precomputed in the offline phase. In the online phase, the small matrix $\mathbf{\Omega }(\mathit{\mu },\mathrm{i}\mathit{\omega }{)}^{-1}+{\mathbf{V}}^{\mathrm{T}}{\left(\mathrm{i}\omega {\mathbf{E}}_0-{\mathbf{A}}_0\right)}^{-1}\mathbf{U}$ can be formed (under the assumption that all necessary matrices are invertible) and used for small and dense system solves in order to calculate $\mathcal{G}(\mathit{\mu },\mathrm{i}\mathit{\omega })$ with only little effort.

5. Numerical results

In this section, we present the numerical results for the two differential-algebraic systems from Section 2. The computations have been done on a computer with 2 Intel Xeon Silver 4110 CPUs running at 2.1 GHz and equipped with 192 GB total main memory. The experiments use MATLAB®R2021a (9.3.0.713579) and examples and methods from M-M.E.S.S.-2.2 (Saak, Köhler, and Benner 2022).

5.1. Problem 1: Stokes equation

We consider a system that describes the creeping flow in capillaries or porous media. It has the following structure

(38) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}\mathit{v}(\mathit{\zeta },\mathit{t})& =\mathit{\mu }\Delta \mathit{v}(\mathit{\zeta },\mathit{t})-\mathrm{\nabla p}(\mathit{\zeta },\mathit{t})+\mathit{f}(\mathit{\zeta },\mathit{t}),\\ 0& =\mathrm{div}(\mathit{v}(\mathit{\zeta },\mathit{t})),\end{array}$(38)

with appropriate initial and boundary conditions. The position in the domain $\mathrm{\Omega }\subset {\mathbb{R}}^{\mathrm{\ell }}$ is described by $\mathit{\zeta }\in \mathrm{\Omega }$ and $\mathit{t}\ge 0$ is the time. For simplicity, we use a classical solution concept and assume that the external force $\mathit{f}:\mathrm{\Omega }\times [0,\mathrm{\infty })\to {\mathbb{R}}^{\mathrm{\ell }}$ is continuous and that the velocities $\mathit{v}:\mathrm{\Omega }\times [0,\mathrm{\infty })\to {\mathbb{R}}^{\mathrm{\ell }}$ and pressures $\mathit{p}:\mathrm{\Omega }\times [0,\mathrm{\infty })\to {\mathbb{R}}^{\mathrm{\ell }}$ satisfy the necessary smoothness conditions. The parameter $\mathit{\mu }\in \mathcal{D}$ represents the dynamic viscosity. We discretize system (38) by finite differences as shown in (Mehrmann and Stykel 2005; Stykel 2006) and add an output equation. Then, we obtain a discretized system of the form (4), which is of index 2. The matrix $\mathit{A}(\mathit{\mu })=\mathit{\mu }\cdot \mathit{A}\in {\mathbb{R}}^{\mathit{n},\mathit{n}}$ depends on the viscosity parameter $\mathit{\mu }$.

Figure 1. Error and error estimates of the approximated controllability Gramians for the Stokes system (4) after the first iteration of the reduced basis method.

Figure 2. Results for the reduction of the Stokes system (4).

Our constant example matrices $\mathit{E}=I_n,\mathit{A},\mathit{G},B_1,C_1$ are created by the M-M.E.S.S. function $\mathrm{stokes}\_\mathrm{FVM}$ and we choose the dimensions $\mathit{n}=3280,\mathit{q}=1681$, $\mathit{m},\mathit{p}=1$, and the parameter set $\mathcal{D}=\left[\frac{1}{2},\frac{3}{2}\right]$. The test parameter set ${\mathcal{D}}_{\mathrm{Test}}$ consists of ten equidistant points within $\mathcal{D}$, i. e., ${\mathcal{D}}_{\mathrm{Test}}=\left\{\frac{1}{2}+\frac{1}{9}\cdot \mathit{k}|\mathit{k}=0,\dots ,9\right\}$. The matrices $B_2$ and $C_2$ are zero matrices. The reduced basis method from Section 4 produces bases of dimensions $24$ for both projected continuous-time Lyapunov equations. This basis generation is done in the offline phase, which takes $2.7\cdot {10}^3$ seconds for this example. For the first iteration of the reduced basis method corresponding to the proper controllability Lyapunov equation, we evaluate the errors and their estimates from Equation (28)(28) $\parallel \text{}\mathfrak{E}(\mathit{\mu }){\parallel }_{\text{F}}=\parallel \text{}\mathit{e}(\mathit{\mu })\text{}\parallel 2\le \parallel \text{ }{\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{L}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu }){\parallel }_2\parallel \mathit{r}(\mathit{\mu })\text{}\parallel 2=1{\sigma }_{{\mathit{min}}^{+}}{\Pi }_{\text{l}}^{\otimes }(\mathit{\mu })\mathit{L}(\mathit{\mu }){\Pi }_{\text{r}}^{\otimes }(\mathit{\mu })\parallel \text{}\mathit{\Re }(\mathit{\mu })\text{}\parallel \text{F}\le \frac{1}{\mathit{\alpha }(\mathit{\mu })}\parallel \text{}\mathit{\Re }(\mathit{\mu })\text{}\parallel \text{F}=:{\Delta }_V^{(1)}(\mathit{\mu }).$(28) and (29(29) $\begin{array}{cc}\parallel \mathfrak{E}(\mathit{\mu }){\parallel }_{\mathrm{F}}& =\parallel \mathfrak{E}(\mathit{\mu })+\hat{\mathfrak{E}}(\mathit{\mu })-\hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}\le \parallel \hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}+\parallel \mathfrak{E}(\mathit{\mu })-\hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}\\ \\ \mathit{\hspace{0.17em}}& \le \parallel \hat{\mathfrak{E}}(\mathit{\mu }){\parallel }_{\mathrm{F}}+\frac{1}{\mathit{\alpha }(\mathit{\mu })}\parallel \hat{\mathrm{\Re }}(\mathit{\mu }){\parallel }_{\mathrm{F}}:={\mathrm{\Delta }}_{\mathbf{V}}^{(2)}(\mathit{\mu }),\end{array}$(29) ). They are presented in Figure 1, where the error is evaluated by

$\mathrm{error}\approx \parallel {\mathbf{P}}_{\mathrm{p}}^{\mathrm{a}\mathrm{c}\mathrm{c}}(\mathit{\mu })-\mathbf{Z}(\mathit{\mu })\mathbf{Z}(\mathit{\mu }{)}^{\mathrm{H}}{\parallel }_{\mathrm{F}},$

where ${\mathbf{P}}_{\mathrm{p}}^{\mathrm{a}\mathrm{c}\mathrm{c}}(\mathit{\mu })$ denotes an accurate approximation of the exact solution ${\mathbf{P}}_{\mathrm{p}}(\mathit{\mu })$. For evaluating the error, the Lyapunov equations are also solved with the residual tolerance ${10}^{-15}$ to obtain an accurate estimate of the exact solution. After the first iteration step, we obtain an error that is already smaller than the tolerance $\mathrm{tol}={10}^{-4}$. Additionally, we see that the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ provides an almost sharp bound of the actual error while ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}$ leads to more conservative error bounds.

After we have determined the bases, we approximate the proper controllability and observability Gramians to apply balanced truncation in the online phase for two different parameters $\mathit{\mu }=1.28$ and $\mathit{\mu }=0.65$ where the first parameter is the eighth entry from the test parameter set ${\mathcal{D}}_{\mathrm{Test}}$. The reduced systems for both parameters are of dimension $\mathit{r}=10$ and consist only of differential states. This online phases last for $3.7\cdot {10}^{-2}$ and $1.9\cdot {10}^{-2}$ seconds, respectively, for this example. Figures 2a,b depict the transfer functions of the full and the reduced systems and the corresponding errors. We observe that the error is smaller than ${10}^{-8}$ for all $\mathit{\omega }$ in $[{10}^{-4},{10}^4]$. Additionally, in Figure 3a,b, we show the output $\mathit{y}(\mathit{t})$ of the full system and the output $y_{\mathrm{R}}(\mathit{t})$ of the reduced system for an input $\mathit{u}(\mathit{t})=-2-sin(\mathit{t})$, $\mathit{t}\in [0,10]$, and the initial condition is given as $\mathbf{z}(0)=0$. We also depict the error $\parallel \mathit{y}(\mathit{t})-y_{\mathrm{R}}(\mathit{t}){\parallel }_2$ which is smaller than ${10}^{-8}$ in the entire time range. Hence, we can observe that the system behaviour is well-approximated by the reduced system.

Figure 3. Output and output error of the original and reduced Stokes system (4).

Now, to demonstrate that the methods presented in this paper are also suitable for systems with an improper transfer function, we modify the system in such a way that the matrices $B_2$ and $C_2$ are chosen to be $B_2={\left[\begin{array}{cccc}0& \dots & 0& 1\end{array}\right]}^{\mathrm{T}}$ and $C_2=\left[\begin{array}{ccccc}1& 0& \dots & 0& 1\end{array}\right].$ That way, we obtain more complex restrictions for $G^{\mathrm{T}}\mathit{x}(\mathit{t})$ and also evaluate the pressure described by $\mathit{\lambda }(\mathit{t})$ to generate an output $\mathit{y}(\mathit{t})$. We additionally assume that the external force represented by the input function $B_1(\mathit{\mu })\mathit{u}(\mathit{t})$ varies its magnitude depending on the viscosity of the system so that the input matrix $B_1(\mathit{\mu })=(1+\mathit{\mu })\cdot B_1$ is parameter-dependent, and hence, the polynomial part of the transfer function is parameter-dependent as well.

Again, in the offline phase, we execute the reduced basis method to generate bases that approximate the proper controllability and observability spaces (i.e., the image space of the corresponding proper Gramians) of this system. This phase takes $9.3\cdot {10}^3$ seconds for this example. In Figure 4a,b, the errors and error estimations in the two stages of the reduced basis method corresponding to the proper controllability space are shown. For evaluating the error, the Lyapunov equations are also solved with the residual tolerance ${10}^{-15}$ to obtain an accurate estimate of the exact solution. We can observe that after the first iteration step, we obtain an error that is significantly larger than the tolerance ${10}^{-4}$ as shown in Figure 4a. Additionally, we see that the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ provides an almost sharp bound of the actual error while ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}$ leads to more conservative error bounds. The second iteration of the reduced basis method and the corresponding errors are depicted in Figure 4b. In this case, the error as well as the estimators are smaller than the tolerance $\mathrm{tol}={10}^{-4}$. The final basis projection space dimensions for the projected continuous-time Lyapunov equations are $48$ and $54$.

The next step is to use the bases to generate approximations of the system Gramians and to apply balanced truncation, which leads to a reduced system of dimension $25$ where we obtain $23$ differential states and $2$ algebraic ones for both parameters $\mathit{\mu }=1.28$ and $\mathit{\mu }=0.65$. The online phase lasts for $3.2$ and $3.1$ seconds if we use balanced truncation via the improper Gramians to reduce the improper system components and $2.4$ seconds for both parameters if we extract the polynomial part of the transfer function. However, in this example, the parameter dependency does not have a low-rank structure. Thus, for the determination of the polynomial part, it is necessary to solve one linear system of equations of full order in the online phase for each parameter of interest. Explicit formulas of the improper Gramians are also described in (Stykel 2006). For a fixed parameter, also an explicit formula

(39) ${\mathbf{S}}_{\mathrm{i}}^{\mathrm{T}}\mathbf{A}{\mathbf{R}}_{\mathrm{i}}=\left[\begin{array}{cc}{B}_{11}& C_2{({\mathit{G}}^{\mathrm{T}}\mathit{G})}^{-1}{B}_2\\ C_2{({\mathit{G}}^{\mathrm{T}}\mathit{G})}^{-1}{B}_2& 0\end{array}\right]$(39)

is provided with

$B_{11}:=C_1\mathit{G}(G^{\mathrm{T}}\mathit{G}{)}^{-1}{B}_2+C_2(G^{\mathrm{T}}\mathit{G}{)}^{-1}{G}^{\mathrm{T}}{B}_{12},B_{12}:=B_1-\mathit{AG}(G^{\mathrm{T}}\mathit{G}{)}^{-1}{B}_2,$

compare with (14). This formula replaces the computation of the low-rank factors of the improper Gramians in the online phase. In our setting, (39(38) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}\mathit{v}(\mathit{\zeta },\mathit{t})& =\mathit{\mu }\Delta \mathit{v}(\mathit{\zeta },\mathit{t})-\mathrm{\nabla p}(\mathit{\zeta },\mathit{t})+\mathit{f}(\mathit{\zeta },\mathit{t}),\\ 0& =\mathrm{div}(\mathit{v}(\mathit{\zeta },\mathit{t})),\end{array}$(38) ) can be evaluated efficiently since all involved matrices can be precomputed and then appropriately scaled with $\mathit{\mu }$ in the parametric setting.

Figure 4. Error and error estimates of the approximated gramians for the stokes system (4) with improper parts.

Figure 5. Results for the reduction of the stokes system (4) with improper parts.

We demonstrate the quality of the reduced systems by evaluating the error of the transfer functions. In the following, we denote the reduced system (error system) with an algebraic part generated by balanced truncation as reduced BT (error BT) and the reduced system (error system) with an algebraic part obtained by approximating the polynomial part of the transfer function as reduced TF (error TF). We show the sigma plot of the original and reduced transfer function as well as the error for two parameters. The first parameter, depicted in Figure 5a, is $\mathit{\mu }=1.28\in {\mathcal{D}}_{\mathrm{Test}}$ that is an element from the test parameter set and hence we already know that the error estimators ${\mathrm{\Delta }}_{\mathbf{V}}^{(1/2)}(\mathit{\mu })$ are smaller than the tolerance and the proper controllability space is well-approximated. The second parameter is $\mathit{\mu }=0.65\in \mathcal{D}$, and the corresponding results are shown in Figure 5b. The polynomial part of the transfer function is clearly observable since it does not tend to zero for growing values of $\mathit{\omega }$. We observe that for both parameters, the error is smaller than ${10}^{-5}$ in the entire frequency band $[{10}^{-4},{10}^6]$ and that the parameter chosen from the test parameter set does not lead to significantly better approximations than an arbitrarily chosen parameter from the parameter set $\mathcal{D}$. This allows the conclusion that the basis ${\mathbf{V}}_{\mathrm{glob}}$ approximates the proper controllability space well – for all parameters in $\mathcal{D}$. Additionally, we can observe that both procedures to approximate the improper part of the system lead to satisfactory results.

Finally, to illustrate that the output $\mathit{y}$ is well approximated by the reduced output $y_{\mathrm{R}}$, in Figure 6a,b we evaluate the norms of the output error $\parallel \mathit{y}(\mathit{t})-y_{\mathrm{R}}(\mathit{t}){\parallel }_2$ for $\mathit{t}\in [0,100]$ for the test parameter $\mathit{\mu }=1.28$ and the parameter $\mathit{\mu }=0.65$. We choose the input function $\mathit{u}(\mathit{t})=5-5\cdot \cos (\mathit{t})$ and the initial condition $\mathbf{z}(0)=\left[\begin{array}{c}\mathit{x}(0)\\ \mathit{\lambda }(0)\end{array}\right]=0$, in such a way that the consistency conditions are satisfied. We obtain an output error that is smaller than ${10}^{-5}$ for all $\mathit{t}\in [0,100]$ for both evaluated parameters. Again, we observe that both approximations of the algebraic part of the system lead to results of similar quality.

Figure 6. Output and output error of the original and reduced Stokes system (4) with improper parts.

5.2. Problem 2: Mechanical system

Figure 7. Constrained mass-damper-spring system including three chains of masses that are connected by the mass $m_{3\mathrm{\ell }+1}$.

As an example for system (5), we consider a constrained mass-spring-damper system which is depicted in Figure 7 and which is similar to the one considered in (Ugrica 2020). The matrices are generated using the function triplechain_MSD from the M-M.E.S.S.-2.2. toolbox, see (Saak, Köhler, and Benner 2022) where the mass matrix $\mathit{M}$ is built as

$\mathit{M}=\left[\begin{array}{cccc}{M}_1& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& M_2& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& M_3& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& m_0\end{array}\right],M_i=m_i\cdot I_{\mathrm{\ell }}\in {\mathbb{R}}^{\mathrm{\ell },\mathrm{\ell }},\mathit{i}=1,2,3$

for $m_1=m_2=m_3=m_0=1$. The stiffness matrix is built as

$\mathit{K}=\left[\begin{array}{cccc}{K}_{11}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& {\kappa }_1\\ \mathit{\hspace{0.17em}}& K_{22}& \mathit{\hspace{0.17em}}& {\kappa }_2\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& K_{33}& {\kappa }_3\\ {\kappa }_1^{\mathrm{T}}& {\kappa }_2^{\mathrm{T}}& {\kappa }_3^{\mathrm{T}}& k_1+k_2+k_3+k_0\end{array}\right],K_{\mathit{ii}}=k_i\left[\begin{array}{ccccc}2& -1& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ -1& 2& -1& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \ddots & \ddots & \ddots & \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& -1& 2& -1\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& -1& 2\end{array}\right]\in {\mathbb{R}}^{\mathrm{\ell },\mathrm{\ell }},$

with ${\kappa }_i={\left[\begin{array}{cccc}0& \dots & 0& k_i\end{array}\right]}^{\mathrm{T}}$ and $k_1=k_2=k_3=k_0=1$. We set $\mathrm{\ell }=200$ so that the dimension of $\mathit{M}$ and $\mathit{K}$ is $\mathit{n}=601$ and the matrix $\mathit{G}\in {\mathbb{R}}^{\mathit{n},1}$ is a zero matrix except for the entries $G_{1,1}=1$ and $G_{\mathit{n},1}=-1$. Hence, the first and the last mass are additionally connected by a rigid structure, as shown in Figure 7. Hence, the first-order system is of dimension $2\mathit{n}+2=1204$ after applying the transformation presented in Subsection 4.2.2. The input matrix $B_x\in {\mathbb{R}}^{\mathit{n},1}$ is chosen as a zero matrix with one nonzero entry $(B_x{)}_{450,1}=1.$ The output matrix $C_x\in {\mathbb{R}}^{1,\mathit{n}}$ is given by $C_x=B_x^{\mathrm{T}}$.

These systems occur in the field of damping optimization, see (Benner, Tomljanović, and Truhar 2011; Truhar, Tomljanović, and Puvača 2019; Truhar and Veselić 2009; Ugrica 2020), where the damping matrix $\mathit{D}\in {\mathbb{R}}^{\mathit{n},\mathit{n}}$ consists of an internal damping $D_{\mathrm{int}}$ and some external dampers that need to be optimized. In our example, we choose Rayleigh damping as a model for the internal damping, i.e., $D_{\mathrm{int}}=\mathit{\eta M}+\mathit{\rho K}$, where we set $\mathit{\eta }=0.01$ and $\mathit{\rho }=0.02$. To describe the external damping, we assume that every mass is connected to a grounded damper, where the masses $m_1,\dots ,m_{\mathrm{\ell }}$ are connected by dampers with viscosity ${\mu }_1$, the masses $m_{\mathrm{\ell }+1},\dots ,m_{2\mathrm{\ell }}$ by dampers with viscosity ${\mu }_2$ and the masses $m_{2\mathrm{\ell }+1},\dots ,m_{3\mathrm{\ell }}$ by dampers with viscosity ${\mu }_3$ as depicted in Figure 7. The mass $m_{3\mathrm{\ell }+1}$ is connected to a grounded damper with fixed viscosity $1$. The external dampers can be described by the matrices

$F_1=\left[\begin{array}{cc}{I}_{\mathrm{\ell }}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& 0_{2\mathrm{\ell }+1,2\mathrm{\ell }+1}\end{array}\right],F_2=\left[\begin{array}{ccc}{0}_{\mathrm{\ell },\mathrm{\ell }}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& I_{\mathrm{\ell }}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& 0_{\mathrm{\ell }+1,\mathrm{\ell }+1}\end{array}\right],F_3=\left[\begin{array}{ccc}{0}_{2\mathrm{\ell },2\mathrm{\ell }}& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& I_{\mathrm{\ell }}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& 0\end{array}\right],F_4=\left[\begin{array}{cc}{0}_{3\mathrm{\ell },3\mathrm{\ell }}& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& 1\end{array}\right],$

where $0_{n_1,n_2}$ denotes the zero matrix in ${\mathbb{R}}^{n_1,n_2}$. Therefore, the overall damping matrix can then be constructed as

$\mathit{D}(\mathit{\mu }):=D_{\mathrm{int}}+{\mu }_1{F}_1+{\mu }_2{F}_2+{\mu }_3{F}_3+F_4,$

where $\mathit{\mu }=[{\mu }_1,{\mu }_2,{\mu }_3]\in \mathcal{D}$ are the viscosities of the dampers, which are optimized and thus variable in our model. This model leads to a differential-algebraic system of index 3. We generate the index-2 surrogate system as shown in Subsection 4.2.2 using the value $\mathit{\gamma }(\mathit{\mu })=\frac{1}{2}{\gamma }^{\ast }(\mathit{\mu })$ from (36).

We consider two parameter settings. In the first one, the system is damped more strongly, and in the second one more weakly. The latter setting is used to illustrate that the slightly damped system is more difficult to approximate by projection spaces of small dimensions, which poses a challenge from the numerical point of view. For both parameter settings, the reduced models are obtained by truncating all state variables corresponding to proper Hankel singular values with ${\sigma }_{\mathrm{p},k}/{\sigma }_{\mathrm{p},1}<{10}^{-8}$.

5.2.1. Setting 1: Stronger external damping

We consider the parameter set $\mathcal{D}=[0.1,1{]}^3$ and choose the test parameter set ${\mathcal{D}}_{\mathrm{Test}}$ to consist of 50 randomly selected parameter triples from $\mathcal{D}$. Therefore, we use the $\mathrm{rand}$ operator from MATLAB so that the $\mathit{k}$-th test parameter is equal to $\mathit{\mu }=0.1+0.9\cdot \mathrm{rand}(3,1)$ with fixed seed.

Figure 8. Error and error estimates of the approximated controllability Gramians of the mechanical system (5) with stronger external damping and ${\mu }_1=0.25$, ${\mu }_3=0.35$ after the first iteration of the reduced basis method.

Figure 9. Results for the reduction of the mechanical system (5) with stronger external damping.

In the offline phase, the approximation of the controllability and observability spaces takes $1.4\cdot {10}^4$ seconds for this example. For the proper controllability space, we evaluate the errors and the error estimates that are depicted in Figure 8 for the first iteration of the reduced basis method. For evaluating the error, the Lyapunov equations (9(8) ${\mathbf{P}}_{\mathrm{p}}={\int }_0^{\mathrm{\infty }}{\mathbf{F}}_{\mathbf{J}}(\mathit{t})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{J}}(\mathit{t}{)}^{\mathrm{T}}\mathrm{d}\mathit{t},{\mathbf{P}}_{\mathrm{i}}=\sum _{\mathit{k}=0}^{\mathit{\nu }-1}{\mathbf{F}}_{\mathbf{N}}(\mathit{k})\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{N}}(\mathit{k}{)}^{\mathrm{T}},$(8) ) are also solved by computing the projection matrices ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })$ and ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$ explicitly and solving the corresponding Lyapunov equation. Therefore, we fix the parameters ${\mu }_1=0.25$ and ${\mu }_3=0.35$ and show the errors for the varying second parameter ${\mu }_2\in [0.1,1]$. After the first iteration and for a basis of dimension $219$, the maximal error estimate ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ is equal to $9.48\cdot {10}^{-5}$. Hence, the basis approximates the proper controllability space sufficiently well. Again, we see that the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ approximates the error $\mathfrak{E}(\cdot )$ better than the estimation by ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}$, which massively overestimates the error. Hence, we only use the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ to determine the basis quality.

We need two iterations of the reduced basis method to determine the respective observability subspace of dimension of 304. Hence, we have to solve Lyapunov equations of dimensions 219 and 304 in the online phase to compute the Gramian approximations. We apply the online phase for the parameters $\mathit{\mu }=\left[\begin{array}{c}0.2622,0.1175,0.5169\end{array}\right]\in {\mathcal{D}}_{\mathrm{Test}}\subset \mathcal{D}$ (which is rounded to 4 digits here) and $\mathit{\mu }=\left[\begin{array}{c}0.15,0.15,0.5\end{array}\right]\in \mathcal{D}$. We obtain the two respective reduced systems of dimension 75 and 63, where all states are differential ones. The online phases take $6.4\cdot {10}^{-1}$ and $1.3\cdot {10}^{-1}$ seconds, respectively, for this example.

Figure 10. Output and output error of the original and reduced system (5) with stronger external damping.

In Figure 9, we depict the original, the reduced, and the error transfer function. Additionally, we evaluate the outputs of the original and the reduced system $\mathit{y}$ and $y_{\mathrm{R}}$. The results are shown in Figure 10, where we observe the norm of the output error $\parallel \mathit{y}(\mathit{t})-y_{\mathrm{R}}(\mathit{t}){\parallel }_2$ for $\mathit{t}\in [0,50]$. We choose the input function $\mathit{u}(\mathit{t})=1-1\cdot cos(\mathit{t})$ and the initial condition $\mathbf{z}(0)=\left[\begin{array}{c}\mathit{x}(0)\\ \mathit{\lambda }(0)\end{array}\right]=0$, so that the consistency conditions are satisfied. We observe that the output error is smaller than ${10}^{-2}$ for all $\mathit{t}\in [0,50]$. Note that for damping optimization examples, the proper controllability and observability spaces are often hard to approximate by spaces of small dimensions, see (Benner et al. 2016; Tomljanović, Beattie, and Gugercin 2018). Hence, the tolerances are chosen rather large for this kind of example to achieve satisfying results.

5.2.2. Setting 2: Weaker external damping

Now, we consider a parameter set with viscosities of small magnitudes, that is $\mathcal{D}=[0.01,0.1{]}^3$. The test-parameter set ${\mathcal{D}}_{\mathrm{Test}}$ consists of 50 randomly selected parameters from $\mathcal{D}$. To generate the test parameters, we again use the $\mathrm{rand}$ operator from MATLAB so that the $\mathit{k}$-th parameter is equal to $\mathit{\mu }=0.01+0.09\cdot \mathrm{rand}(3,1)$ with fixed seed.

Figure 11. Error and error estimates of the approximated observability Gramians of the mechanical system (5) with weaker external damping and ${\mu }_1=0.025$, ${\mu }_3=0.035$ after the first iteration of the reduced basis method.

Figure 12. Results for the reduction of the mechanical system (5) with weaker external damping.

First, we apply the offline phase of the reduced basis method that takes $6.9\cdot {10}^3$ seconds. When computing the proper controllability space, two iterations are needed to generate a subspace of dimension 449 that leads to error estimates ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}(\mathit{\mu })$ that are at most $5.6\cdot {10}^{-4}$. For the proper observability space, the errors and the error estimates are illustrated in Figure 11 within the first iteration of the reduced basis method with fixed parameters ${\mu }_1=0.025$ and ${\mu }_3=0.035$ and varying ${\mu }_2\in [0.01,0.1]$. For evaluating the error, the Lyapunov equations (10(9) $\mathbf{E}{\mathbf{P}}_{\mathrm{p}}{\mathbf{A}}^{\mathrm{T}}+\mathbf{A}{\mathbf{P}}_{\mathrm{p}}{\mathbf{E}}^{\mathrm{T}}=-{\mathbf{\Pi }}_{\mathrm{l}}\mathbf{B}{\mathbf{B}}^{\mathrm{T}}{{\mathbf{\Pi }}_{\mathrm{l}}}^{\mathrm{T}},{\mathbf{P}}_{\mathrm{p}}={\mathbf{\Pi }}_{\mathrm{r}}{\mathbf{P}}_{\mathrm{p}}{{\mathbf{\Pi }}_{\mathrm{r}}}^{\mathrm{T}},$(9) ) are also solved by computing the projection matrices ${\mathbf{\Pi }}_{\mathrm{l}}(\mathit{\mu })$ and ${\mathbf{\Pi }}_{\mathrm{r}}(\mathit{\mu })$ explicitly and solving the corresponding Lyapunov equation. After the first iteration, we determine a basis of dimension 223 for which the maximal error estimate ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ is equal to $5.8\cdot {10}^{-4}$, and hence, sufficient accuracy is achieved. As for the first parameter setting, the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ approximates the error $(\cdot )$ better than the estimation by ${\mathrm{\Delta }}_{\mathbf{V}}^{(1)}$. Hence, we use the error estimator ${\mathrm{\Delta }}_{\mathbf{V}}^{(2)}$ within the reduced basis method. We apply the online phase for the parameters $\mathit{\mu }=\left[\begin{array}{c}0.0886,0.0972,0.0882\end{array}\right]\in {\mathcal{D}}_{\mathrm{Test}}\subset \mathcal{D}$ (which is rounded to 4 digits) and $\mathit{\mu }=\left[\begin{array}{c}0.01,0.02,0.01\end{array}\right]\in \mathcal{D}$ to derive a reduced-order model of dimension 78 for the first parameter and the reduced dimension 178 for the second one. The online phases need 2.21 and 2.04 seconds, respectively.

Figure 13. Output and output error of the original and reduced mechanical system (5) with weaker external damping.

Figure 12 depicts the original, the reduced, and the error transfer function. Additionally, we evaluate the norm of the output error $\parallel \mathit{y}(\mathit{t})-y_{\mathrm{R}}(\mathit{t}){\parallel }_2$ for $\mathit{t}\in [0,50]$ in Figure 13, where the input function $\mathit{u}(\mathit{t})=1-1\cdot \cos (\mathit{t})$ and the initial condition $\mathbf{z}(0)=\left[\begin{array}{c}\mathit{x}(0)\\ \mathit{\lambda }(0)\end{array}\right]=0$ are chosen in such a way that the consistency conditions are satisfied. We observe that the output error is smaller than ${10}^{-4}$ for all $\mathit{t}\in [0,50]$.

6. Conclusions

This paper has addressed the model reduction of parameter-dependent differential-algebraic systems by balanced truncation. To apply the balancing procedure, we have utilized specific projectors to eliminate the algebraic constraints. This has enabled us to compute the necessary Gramians efficiently by solving projected Lyapunov equations.

To handle the parameter-dependency, we have applied the reduced basis method, which is split into the offline phase and the online phase. In the offline phase, we have computed the basis of the subspace which approximates the solution space of the parametric Lyapunov equations. To assess the quality of this subspace, we have derived two error estimators. Afterwards, in the online phase, we have solved a reduced Lyapunov equation to obtain an approximation of the Gramian for a parameter value of interest efficiently. Therefore, a balanced truncation for this parameter value can also be carried out fast.

This method has been illustrated by numerical examples of index two and three. In particular, we were often able to reduce the associated Lyapunov equations to very small dimensions. We have evaluated our error estimators, where the second one can often estimate the error almost exactly.

Novelty statement

We propose a new method that reduces parametric differential-algebraic equations by combining projection methods for differential-algebraic equations and the reduced basis methods for Lyapunov equations. To apply the reduced basis method, new error estimators are invented. The performance of our new method is illustrated for several examples.

Code availability

The code and data that has been used to generate the numerical results of this work are freely available under the DOI 10.5281/zenodo.10040792 under the 3-clause BSD license and is authored by Jennifer Przybilla.

CRediT author statement

Jennifer Przybilla: Methodology, Software, Data Curation, Writing – Original Draft, Writing – Review & Editing, Visualization; Matthias Voigt: Conceptualization, Writing – Original Draft, Writing – Review & Editing, Supervision.

Disclosure statement

No potential conflict of interest was reported by the author(s).