# Balancing-based model reduction for switched descriptor systems

Mattia Manucci⋆ and Benjamin Unger⋆ ⋆ Stuttgart Center for Simulation Science
(SC SimTech), University of Stuttgart, Universitätsstr. 32, 70569 Stuttgart,
Germany<EMAIL_ADDRESS>
###### Abstract.

We present a novel certified model order reduction (MOR) algorithm for
switched descriptor systems applicable to large-scale systems. Our algorithm
combines the idea of [Hossain & Trenn, Technical report, 2023] to reformulate
the switched descriptor system as a switched ordinary differential equation
with jumps and an extension of the balanced truncation for switched ODE from
[Pontes Duff et al., IEEE Trans. Automat. Control, 2020]. Besides being the
first MOR method for switched descriptor systems applicable to the large-scale
setting, we give a detailed numerical analysis by incorporating the error in
the computation of the system Gramians in the a-priori error bound for the
output of the reduced system. In more detail, we demonstrate, theoretically
and numerically, that the standard error bound is not applicable, and a
certificate must account for the numerical approximation errors.

Keywords: model order reduction, differential-algebraic equations, switched
systems, balanced truncation, error bound, generalized Lyapunov equation

AMS subject classification: 65F45, 65F55, 65P99, 93A30, 93A15, 93B99

##  1\. Introduction

Consider the system of switched _differential-algebraic equations_ (DAE) of
the form

$\Sigma_{q}\quad\left\\{\quad\begin{aligned}
\bm{E}_{q(t)}\dot{\bm{x}}(t)&=\bm{A}_{q(t)}\bm{x}(t)+\bm{B}_{q(t)}\bm{u}(t),&\bm{x}(t_{0})&=\mathbf{0},\\\
\bm{y}(t)&=\bm{C}_{q(t)}\bm{x}(t),\\\ \end{aligned}\right.$ (1.1)

where $q\colon\mathbb{R}\to\mathcal{J}\vcentcolon=\\{1,\ldots,M\\}$ is the
external switching signal, which we assume to be an element of the set of
allowed switching signals

$\mathcal{S}\vcentcolon=\\{q\colon\mathbb{R}\to\mathcal{J}\mid q\text{ is
right continuous with locally finite number of jumps}\\}.$ (1.2)

The symbols $\bm{x}(t)\in\mathbb{R}^{n}$, $\bm{u}(t)\in\mathbb{R}^{m}$, and
$\bm{y}(t)\in\mathbb{R}^{p}$ denote the _state_ , the controlled _input_ , and
the measured _output_ , respectively. The system matrices
$\bm{E}_{j}\in\mathbb{R}^{n\times n}$, $\bm{A}_{j}\in\mathbb{R}^{n\times n}$,
$\bm{B}_{j}\in\mathbb{R}^{n\times m}$, and $\bm{C}_{j}\in\mathbb{R}^{p\times
n}$ correspond to the _differential-algebraic equation_ (DAE) active in mode
$j\in\mathcal{J}$. We emphasize that $\bm{E}_{j}$ might be singular and we
assume that the finite eigenvalues $\lambda$ of the matrix pair
$(\bm{E}_{j},\bm{A}_{j})$ have negative real part for all $j\in\mathcal{J}$.
We refer to (1.1) as the _full-order model_ (FOM). Sample applications include
robot manipulators, traffic management, automatic gear shifting, and power
systems; see for instance [8] and the references therein.

If (1.1) has to be evaluated repeatedly, for instance in simulation context
for different inputs or switching signals, or if matrix equalities or
inequalities in the context of synthesis have to be solved, then a large
dimension $n$ of the state renders this a computationally expensive task. In
such scenarios, one can rely on _model order reduction_ (MOR) and replace
(1.1) by the _reduced-order model_ (ROM)

$\tilde{\Sigma}_{q}\quad\left\\{\quad\begin{aligned}
\tilde{\bm{E}}_{q(t)}\dot{\tilde{\bm{x}}}(t)&=\tilde{\bm{A}}_{q(t)}\tilde{\bm{x}}(t)+\tilde{\bm{B}}_{q(t)}\bm{u}(t),&\tilde{\bm{x}}(t_{0})&=\mathbf{0},\\\
\tilde{\bm{y}}(t)&=\tilde{\bm{C}}_{q(t)}\tilde{\bm{x}}(t),\\\
\end{aligned}\right.$ (1.3)

with $\tilde{\bm{E}}_{j},\tilde{\bm{A}}_{j}\in\mathbb{R}^{r\times r}$,
$\tilde{\bm{B}}_{j}\in\mathbb{R}^{r\times m}$, and
$\tilde{\bm{C}}_{j}\in\mathbb{R}^{p\times r}$, and $r\ll n$. In many cases,
see for instance [1], the reduced system matrices are obtained via
Petrov–Galerkin projection, i.e., one constructs matrices
$\bm{V},\bm{W}\in\mathbb{R}^{n\times r}$ and then defines

$\displaystyle\tilde{\bm{E}}_{j}$
$\displaystyle\vcentcolon=\bm{W}^{\mathsf{T}}\bm{E}_{j}\bm{V},$
$\displaystyle\tilde{\bm{A}}_{j}$
$\displaystyle\vcentcolon=\bm{W}^{\mathsf{T}}\bm{A}_{j}\bm{V},$
$\displaystyle\tilde{\bm{B}}_{j}$
$\displaystyle\vcentcolon=\bm{W}^{\mathsf{T}}\bm{B}_{j},$
$\displaystyle\tilde{\bm{C}}_{j}$ $\displaystyle\vcentcolon=\bm{C}_{j}\bm{V}.$
(1.4)

The goal of MOR is thus to derive in a computationally efficient and robust
way the matrices $\bm{W},\bm{V}$ such that the error $\bm{y}-\tilde{\bm{y}}$
is small in some given norm.

###  1.1. Main contributions

We combine a reformulation of the switched DAE (1.1) taken from [20, 21] with
a generalization of the MOR method from [39] to obtain a novel model reduction
scheme for switched DAEs in Section 3. This is the first MOR method for
switched DAEs that can be applied in a large-scale setting (cf. Section 4) and
provides an a-priori error bound Theorem 4.14. To enable our MOR method in a
large-scale setting, we apply a stationary algorithm in combination with a
subspace projection framework to solve generalized Lyapunov equations. We
provide an in-depth error analysis in Section 4 to certify the numerical
algorithms. In particular, we demonstrate that the classical balanced
truncation error bound is corrupted by the numerical approximation of the
Gramians (cf. Proposition 4.9) such that the error bound has to be modified to
account for the numerical errors; see Theorem 4.14. The modified error bound
employs bounds on the approximation error for the generalized Lyapunov error,
which we provide in Theorem 4.8. The results are verified in Section 5 on a
constrained mass-spring-damper system (cf. Section 5.1) and a discretized
Stokes equation (cf. Section 5.2).

###  1.2. Literature review and state-of-the-art

Our method is mainly inspired by [39] and [21], with [21] based on the thesis
[20], appears to be the only reference for MOR for switched descriptor
systems. Nevertheless, [21] is restricted to switched DAEs where the switching
sequence is known a priori and no error certification is provided. Moreover,
the authors of [21] explicitly mention that in the current version, their
algorithm is limited to small or medium-sized problems and thus not readily
available for the large-scale setting. Let us also mention [41], where a
dimensionality reduction of a switched descriptor system is realized via one
of the Wong sequences defined in (2.3), which is subsequently used to obtain
stability conditions [42].

The main idea of [20] is to reformulate the switched DAE (1.1) as a switched
ODE with jumps and impulses and subsequently apply the midpoint-based balanced
truncation [22] on this system. The literature for MOR for switched ODEs
involves several different approaches. In particular, the authors of [34, 35]
propose constructing the ROMs for each mode independent of the other modes,
completely ignoring the transition from one mode to another. Moreover, if a
state-dependent switching signal is allowed, any approximation of the switched
DAE may be arbitrarily bad. Both phenomena are detailed with examples in [46].
In [61], an approach based on a set of coupled linear matrix inequalities is
proposed, which becomes infeasible in a large-scale context. Nevertheless, the
matrix inequalities can be used to guarantee quadratic stability of the
reduced system and to derive an error bound [38]. Instead of matrix
inequalities, [48, 16] propose solving a set of coupled Lyapunov equations to
compute Gramians, which can then be used for a balancing-based model
reduction. There is no guarantee that this approach has a solution [28], and
in a large-scale setting, the computational complexity may be very demanding.
If the Gramians for each mode can be simultaneously diagonalized, then
classical balanced truncation methods can be adapted as discussed in [33]. In
contrast, the methods reported in [46, 39] rely on reformulating the switched
system as a non-switched system by suitably interpreting the switching signal
as a control input. In [46], the switched system is recast as a linear system
such that standard methods can be applied, while in [39], the system is recast
as a bilinear system, and a balanced truncation approach for such systems is
employed. In this paper, we closely follow the second strategy. We emphasize
that both methods provide an a-priori error certification. If prior
information about the switching sequence is available, then the method in [39]
can be further specialized as reported in [17]. Although we do not pursue this
approach, our method can be adapted similarly. Interpolation-based techniques
are discussed in [43] for hybrid and in [2] for switched systems. For a data-
driven approach to obtain reduced models, we refer to [15]. Let us emphasize
that model reduction is closely related to realization theory. For switched
systems, the associated connections are illustrated in [37] and the references
therein. We conclude by mentioning [36], where the authors discretize the
control variable to obtain a switched system of autonomous PDEs that they
approximate with ROMs in a model predictive control framework.

As one of our main contributions relates to the numerical solution of a GLE,
we also provide a short overview of numerical methods for GLEs. The existence
and uniqueness of solutions of GLEs are established, for instance, in [12,
Thm. 3.6.1], which is essentially based on operator splitting ideas developed
already in [45]. We also mention [60, Lem. 4.2], and for the more general case
of generalized Silvester equations, a similar characterization of solvability
in [23, Thm. 2.1]. The interpretation of the solution as Gramians of bilinear
and stochastic linear control systems and their relation to energy functionals
is discussed in [5]. Conditions for a fast singular value decay of the
solution matrix, essential for a low-rank approximation of the solution, have
been developed in [4, 23]. In terms of numerical schemes for the solution of
GLEs, we mention the ADI-preconditioned Krylov subspace method [11], the
bilinear ADI method [4], the alternative linear scheme [24], which can be
interpreted as finding $\mathcal{H}_{2}$ optimal search directions [3, 7], and
the stationary iteration from [49], which we employ here.

###  1.3. Organization of the manuscript

After this introduction, we recall important concepts from the literature
about DAEs and switched descriptor systems in Sections 2.1 and 2.2,
respectively. Our novel MOR algorithm and its subsequent analysis is presented
in Section 3. To enable our algorithm in a large-scale setting, we discuss
several numerical aspects in Section 4 and present two examples in Section 5.

###  1.4. Notation

The symbols $\bm{I}$ and $\mathrm{GL}_{n}$ denote the identity matrix of
appropriate size and the set of $n\times n$ real nonsingular matrices,
respectively. Let $\bm{A}\in\mathbb{R}^{n\times n}$. Then we call $A$
asymptotically stable if all eigenvalues of $\bm{A}$ are contained in the open
left-half complex plane. Moreover, we write $\bm{A}\succeq\mathbf{0}$ or
$\bm{A}\preceq\mathbf{0}$ if $\bm{A}$ is positive or negative semidefinite,
respectively. The preimage of $A$ with respect to a linear subspace
$\mathcal{N}\subseteq\mathbb{R}^{n}$ is denoted with

$\bm{M}^{-1}(\mathcal{N})\coloneqq\\{\bm{x}\in\mathbb{R}^{n}\mid\bm{M}\bm{x}\in\mathcal{N}\\}.$
(1.5)

The smallest $\bm{A}$-invariant subspace containing $\mathcal{N}$ and the
largest $\bm{A}$-invariant subspace contained in $\mathcal{N}$ are

$\displaystyle\langle\bm{A}\mid\mathcal{N}\rangle$
$\displaystyle\coloneqq\mathcal{N}+\bm{A}\mathcal{N}+\ldots+\bm{A}^{n-1}\mathcal{N},$
$\displaystyle\langle\mathcal{N}\mid\bm{A}\rangle$
$\displaystyle\coloneqq\mathcal{N}\cap\bm{A}^{-1}\mathcal{N}\cap\ldots\cap\bm{A}^{-(n-1)}\mathcal{N}.$

##  2\. Preliminaries

###  2.1. Differential-algebraic equations

Since for constant $q$ the switched DAE (1.1) reduces to a DAE, we also define
the non-switched descriptor system

$\Sigma_{\textsf{DAE}}\quad\left\\{\quad\begin{aligned}
\bm{E}\dot{\bm{x}}(t)&=\bm{A}\bm{x}(t)+\bm{B}\bm{u}(t),\quad\bm{x}(t_{0})=\bm{x}_{0},\\\
\bm{y}(t)&=\bm{C}\bm{x}(t),\\\ \end{aligned}\right.$ (2.1)

with $\bm{E}$, $\bm{A}\in\mathbb{R}^{n\times n}$,
$\bm{B}\in\mathbb{R}^{n\times m}$, $\bm{C}\in\mathbb{R}^{p\times n}$, and
initial value $\bm{x}_{0}\in\mathbb{R}^{n}$. To ensure existence and
uniqueness of solutions of the DAE (2.1), the matrix pair $(\bm{E},\bm{A})$
has to satisfy certain properties; see for instance [25, Cha. 2]. In more
detail, we assume that the matrix pair $(\bm{E},\bm{A})$ is _regular_ , i.e.,
$\det(s\bm{E}-\bm{A})\in\mathbb{C}[s]\setminus\\{0\\}$. In this case, one can
show that in the space of piecewise-smooth distributions [55] the initial
trajectory problem associated with the DAE (2.1) has a unique solution for any
initial value and any right-hand side. Regularity can be characterized by the
Weierstraß form [14] or the slightly simplified _quasi-Weierstraß form_ (QWF)
[6].

######  Theorem 2.1 (Quasi-Weierstrass Form, [6]).

A matrix pair $(\bm{E},\bm{A})\in\mathbb{R}^{n\times
n}\times\mathbb{R}^{n\times n}$ is regular if and only if there exists
matrices $\bm{S},\bm{T}\in\mathrm{GL}_{n}$ such that

$\left(\bm{S}\bm{E}\bm{T},\bm{S}\bm{A}\bm{T}\right)=\Bigg{(}\begin{bmatrix}\bm{I}&\mathbf{0}\\\
\mathbf{0}&\bm{N}\end{bmatrix},\begin{bmatrix}\bm{J}&\mathbf{0}\\\
\mathbf{0}&\bm{I}\end{bmatrix}\Bigg{)},$ (2.2)

where $\bm{N}\in\mathbb{R}^{n_{\bm{N}}\times n_{\bm{N}}}$ is nilpotent with
nilpotency index $\nu$ and $\bm{J}\in\mathbb{R}^{n_{\bm{J}}\times
n_{\bm{J}}}$, with $n_{\bm{J}}=n-n_{\bm{N}}$.

######  Remark 2.2.

The assumption that the finite eigenvalues of the matrix pair
$(\bm{E},\bm{A})$ have negative real part implies that the matrix $\bm{J}$
appearing in (2.2) is asymptotically stable.

The matrices $\bm{S},\bm{T}\in\mathrm{GL}_{n}$ can be constructed via the
_Wong sequences_ [59], which are defined as

$\displaystyle\mathcal{V}^{0}$ $\displaystyle\coloneqq\mathbb{R}^{n},$
$\displaystyle\mathcal{V}^{i+1}$
$\displaystyle\coloneqq\bm{A}^{-1}(\bm{E}\mathcal{V}^{i}),$ $\displaystyle
i\in\mathbb{N},$ (2.3a) $\displaystyle\mathcal{W}^{0}$
$\displaystyle\coloneqq\\{0\\},$ $\displaystyle\mathcal{W}^{j+1}$
$\displaystyle\coloneqq\bm{E}^{-1}(\bm{A}\mathcal{W}^{j}),$ $\displaystyle
j\in\mathbb{N},$ (2.3b)

where we use the notation for the preimage as in (1.5). After finitely many
steps the sequences in (2.3) converge and the limits are given by

$\mathcal{V}^{\star}\coloneqq\bigcap_{i\in\mathbb{N}}\mathcal{V}^{i}\qquad\text{and}\qquad\mathcal{W}^{\star}\coloneqq\bigcup_{i\in\mathbb{N}}\mathcal{W}^{i}.$
(2.4)

######  Theorem 2.3 (QWF via Wong sequences, [6, Thm. 2.6]).

Consider a regular matrix pair $(\bm{E},\bm{A})$ with corresponding Wong
limits $\mathcal{V}^{\star}$ and $\mathcal{W}^{\star}$. For any full rank
matrices $\hat{\mathbf{V}}$ and $\hat{\mathbf{W}}$ such that
$\operatorname{img}(\hat{\mathbf{V}})=\mathcal{V}^{\star}$ and
$\operatorname{img}(\hat{\mathbf{W}})=\mathcal{W}^{\star}$, the matrices

$\bm{T}=[\hat{\mathbf{V}},\hat{\mathbf{W}}],\quad\bm{S}=[\bm{E}\hat{\mathbf{V}},\bm{A}\hat{\mathbf{W}}]^{-1}$
(2.5)

are invertible and transform $(\bm{E},\bm{A})$ into QWF (2.2).

With these preparations, we can aim for a geometric description of the
solution of (2.1) independent on the particular choice of coordinates for
$\mathcal{V}^{\star}$ and $\mathcal{W}^{\star}$. Let us decompose the solution
$\bm{x}$ of (2.1) as
$\bm{x}=\bm{x}^{\mathrm{diff}}\oplus\bm{x}^{\mathrm{imp}}$ with
$\bm{x}^{\mathrm{diff}}(t)\in\mathcal{V}^{\star}$ and
$\bm{x}^{\mathrm{imp}}(t)\in\mathcal{W}^{\star}$ for all $t\in\mathbb{R}$ and
define the matrices

$\displaystyle\mathbf{\Pi}_{(\bm{E},\bm{A})}$
$\displaystyle\coloneqq\bm{T}\begin{bmatrix}\bm{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\bm{T}^{-1},$
$\displaystyle\mathbf{\Pi}^{\mathrm{diff}}_{(\bm{E},\bm{A})}$
$\displaystyle\coloneqq\bm{T}\begin{bmatrix}\bm{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\bm{S},$
$\displaystyle\mathbf{\Pi}^{\mathrm{imp}}_{(\bm{E},\bm{A})}$
$\displaystyle\coloneqq\bm{T}\begin{bmatrix}\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\bm{I}\end{bmatrix}\bm{S},$ (2.6)

which are called the _consistency projector_ , the _differential selector_ ,
and the _impulse selector_ , respectively. Note that the projectors do not
depend on the specific choice of $\bm{S},\bm{T}\in\mathrm{GL}_{n}$ (see [54,
Sec. 4.2.2]) and hence the matrices

$\displaystyle\bm{A}^{\mathrm{diff}}$
$\displaystyle\coloneqq\mathbf{\Pi}^{\mathrm{diff}}_{(\bm{E},\bm{A})}\bm{A},$
$\displaystyle\bm{B}^{\mathrm{diff}}$
$\displaystyle\coloneqq\mathbf{\Pi}^{\mathrm{diff}}_{(\bm{E},\bm{A})}\bm{B},$
$\displaystyle\bm{C}^{\mathrm{diff}}$
$\displaystyle\coloneqq\bm{C}\mathbf{\Pi}_{(\bm{E},\bm{A})},$ (2.7a)
$\displaystyle\bm{E}^{\mathrm{imp}}$
$\displaystyle\coloneqq\mathbf{\Pi}^{\mathrm{imp}}_{(\bm{E},\bm{A})}\bm{E},$
$\displaystyle\bm{B}^{\mathrm{imp}}$
$\displaystyle\coloneqq\mathbf{\Pi}^{\mathrm{imp}}_{(\bm{E},\bm{A})}\bm{B},$
$\displaystyle\bm{C}^{\mathrm{imp}}$
$\displaystyle\coloneqq\bm{C}\left(\bm{I}-\mathbf{\Pi}_{(\bm{E},\bm{A})}\right),$
(2.7b)

are also independent of $\bm{S}$ and $\bm{T}$.

For the switched DAE (1.3) we define accordingly
$\mathbf{\Pi}_{k}\vcentcolon=\mathbf{\Pi}_{(\bm{E}_{k},\bm{A}_{k})}$,
$\mathbf{\Pi}_{k}^{\mathrm{diff}}\vcentcolon=\mathbf{\Pi}^{\mathrm{diff}}_{(\bm{E}_{k},\bm{A}_{k})}$,
$\mathbf{\Pi}_{k}^{\mathrm{imp}}\vcentcolon=\mathbf{\Pi}^{\mathrm{imp}}_{(\bm{E}_{k},\bm{A}_{k})}$,
and the corresponding matrices $\bm{A}_{k}^{\mathrm{diff}}$,
$\bm{B}_{k}^{\mathrm{diff}}$, $\bm{C}_{k}^{\mathrm{diff}}$,
$\bm{E}_{k}^{\mathrm{imp}}$, $\bm{B}_{k}^{\mathrm{imp}}$, and
$\bm{C}_{k}^{\mathrm{imp}}$ for $k\in\mathcal{J}$.

######  Remark 2.4.

Note that the matrices (2.7) and (2.5) should not be formed explicitly. In
particular, in a large-scale context, the sparsity of the system matrices of
the switched system may be destroyed by applying the operators in (2.6). For
this reason, assuming a sparse representation of (2.5) (see the forthcoming
discussion in Section 4.1), the matrices (2.7) are only used implicitly
through suitable matrix-vector products and sparse linear system solves.

######  Remark 2.5.

For the existence of a continuously differentiable solution the input $\bm{u}$
needs to be sufficiently smooth and the initial value
$\bm{x}_{0}\in\mathbb{R}^{n}$ has to satisfy the _consistency condition_

$(\bm{I}-\mathbf{\Pi}_{(\bm{E},\bm{A})})\bm{x}_{0}=\sum_{i=0}^{\nu-1}\left(\bm{E}^{\mathrm{imp}}\right)^{i}\bm{B}^{\mathrm{imp}}\bm{u}^{(i)}(t_{0}).$
(2.8)

###  2.2. Switched descriptor systems

For a switched DAE (1.1), there is no guarantee that the initial value after a
switch satisfies the consistency condition (2.8). Hence, we cannot expect to
assume the existence of classical solutions. Instead, we rely on _piecewise-
smooth distributions_ , first presented in [55] and applied to switched DAEs
in [54, 56]. In particular, regularity of each mode of the switched DAE and
assuming that the switching times do not accumulate is sufficient for the
existence of a unique solution of the initial trajectory problem for (1.1);
see [56, Cor. 5.2]. For an overview of other distributional solution concepts
in the context of DAEs, we refer to [57]. Asymptotic stability of switched
DAEs and related reduced order model is beyond the proposes of this work
although fundamental for practical applications. We refer to [29] for
sufficient conditions for asymptotic stability of switched DAEs under our
class of switching signals defined in (1.2) or under sufficiently slow average
dwell-time switching.

To formulate our forthcoming MOR method, we make use of an equivalent
characterization of the switched descriptor system as a switched ODE with
jumps and impulses, which was first derived in [20, 21]. In more detail, let
$q\in\mathcal{S}$ and consider the switched ODEs with jumps and Dirac impulses
given by

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\bm{z}}(t)&=\bm{A}^{\mathrm{diff}}_{q_{k}}\bm{z}(t)+\bm{B}_{q_{k}}^{\mathrm{diff}}\bm{u}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{z}(t_{k}^{+})&=\mathbf{\Pi}_{q_{k}}\bm{z}(t_{k}^{-})+\mathbf{\Pi}_{q_{k}}\bm{U}_{q_{k}}^{-}(t_{k}^{-}),&\bm{z}(t_{0}^{-})=0,\\\
\bm{y}(t)&=\bm{C}^{\mathrm{diff}}_{q_{k}}\bm{z}(t)+\bm{D}_{q_{k}}\bm{U}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{y}[t_{k}]&=-\bm{C}^{\mathrm{imp}}_{q_{k}}\sum_{i=1}^{\nu_{k}-1}(\bm{E}_{q_{k}}^{\mathrm{imp}})^{i}\bm{z}(t_{k}^{-})\delta_{t_{k}}^{(i)}+\bm{U}_{q_{k}}^{\mathrm{imp}}(t_{k}),\end{aligned}\right.$
(2.9)

where $q_{k}\vcentcolon=q(t_{k})$ is constant along $(t_{k},t_{k+1})$,
$\nu_{q_{k}}$ is the nilpotency index of matrix pencil
$(E_{q_{k}},A_{q_{k}})$, and

$\displaystyle\bm{U}_{q_{k}}^{-}(t_{k}^{-})$
$\displaystyle\vcentcolon=\sum_{i=0}^{\nu_{k}-1}\left(\bm{E}_{q_{k-1}}^{\mathrm{imp}}\right)^{i}\bm{B}_{q_{k-1}}^{\mathrm{imp}}\bm{u}^{(j)}(t_{k}^{-}),$
(2.10a) $\displaystyle\bm{U}_{q_{k}}^{\mathrm{imp}}(t_{k})$
$\displaystyle\vcentcolon=\bm{C}^{\mathrm{imp}}_{q_{k}}\sum_{i=1}^{\nu_{k}-1}\left(\bm{E}_{q_{k}}^{\mathrm{imp}}\right)^{i}\left(\bm{U}_{q_{k+1}}^{-}(t_{k}^{+})-\bm{U}_{q_{k}}^{-}(t^{-}_{k})\right)\delta_{t_{k}}^{(i)},$
(2.10b) $\displaystyle\bm{U}(t)$
$\displaystyle\vcentcolon=\left[\bm{u}^{\mathsf{T}}(t),\dot{\bm{u}}^{\mathsf{T}}(t),\ldots,\bm{u}^{{(n-1)}^{\mathsf{T}}}(t)\right]^{\mathsf{T}},$
(2.10c) $\displaystyle\bm{D}_{q_{k}}$
$\displaystyle\vcentcolon=-\bm{C}^{\mathrm{imp}}_{q_{k}}[\bm{E}_{q_{k}}^{\mathrm{imp}}\bm{B}_{q_{k}}^{\mathrm{imp}},\ldots,(\bm{E}_{q_{k}}^{\mathrm{imp}})^{\nu_{k}-1}\bm{B}_{q_{k}}^{\mathrm{imp}}].$
(2.10d)

Note that we do not recall the precise definition of left-and right-sided
evaluation at $t$, denoted with $t^{-}$ and $t^{+}$, the impulsive component
denoted with $[t]$ here, and the $i$th derivative of the Dirac impulse
$\delta_{t}^{(i)}$, since these aspects are not relevant in the remainder of
the manuscript. Instead, we refer to [55].

######  Theorem 2.6 ([20, Thm. 7.11]).

Assume that the switched DAE (1.1) is regular. Then for every switching signal
$q\in\mathcal{S}$, the systems (1.1) and (2.9) have the same input-output
behavior in the space of piecewise-smooth distributions.

Let us emphasize that in contrast to a standard switched linear system, the
switched system (2.9) has the following additional features:

  1. (i)

The state transition in (2.9) from one mode to another may depend on
derivatives of the input such that additional input-dependent jumps may occur
in the solution.

  2. (ii)

As a result of the impulsive component of the original switched descriptor
system, derivatives of the input may appear as a direct feedthrough term in
the output of (2.9).

  3. (iii)

Inconsistent states before a switch may result in additional Dirac impulses in
the output of (2.9).

For the moment, let us consider the following assumptions; see [20, Thm.
7.12].

######  Assumption 2.7.

Consider the switched system (2.9).

  1. (i)

The system matrices satisfy
$\mathbf{\Pi}_{q_{k}}(\bm{E}_{q_{k-1}}^{\mathrm{imp}})^{i}\bm{B}_{q_{k-1}}^{\mathrm{imp}}=\mathbf{0}$
for all $i\geq 0$ and for all $k\in\mathcal{J}$.

  2. (ii)

The system matrices satisfy
$\bm{C}^{\mathrm{imp}}_{q_{k}}(\bm{E}_{q_{k}}^{\mathrm{imp}})^{i}\bm{\Pi}_{q_{k-1}}=\mathbf{0}$
for all $i\geq 0$ and for all $k\in\mathcal{J}$.

Under 2.7 the switched ODE with jumps (2.9) (and thus (1.1)) simplifies to

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\bm{z}}(t)&=\bm{A}_{q_{k}}^{\mathrm{diff}}\bm{z}(t)+\bm{B}_{q_{k}}^{\mathrm{diff}}\bm{u}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{z}(t_{k}^{+})&=\mathbf{\Pi}_{q_{k}}\bm{z}(t_{k}^{-}),&\bm{z}(t_{0}^{-})=\mathbf{0},\\\
\bm{y}(t)&=\bm{C}_{q_{k}}^{\mathrm{diff}}\bm{z}(t)+\bm{D}_{q_{k}}\bm{U}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{y}[t_{k}]&=\bm{U}_{q_{k}}^{\mathrm{imp}}(t_{k});\end{aligned}\right.$
(2.11)

see [21, Cor. 5].

###  2.3. Reachability and Observability

Let $\bm{\phi}(t,t_{0},\bm{x}_{0},\bm{u},q)$ denote the state trajectory at
time $t$ of the switched system (2.11) starting from
$\bm{x}(t_{0})=\bm{x}_{0}$ with input $\bm{u}$ and switching path
$q\in\mathcal{S}$ with $\mathcal{S}$ given in (1.2).

######  Definition 2.8.

Let $q\in\mathcal{S}$ be a given switching path. A state
$\bm{x}\in\mathbb{R}^{n}$ is called

  1. (i)

_reachable via $q$_ if there exist a time instant $t_{\mathrm{f}}>t_{0}$, and
an input $\bm{u}\colon[t_{0},t_{f}]\rightarrow\mathbb{R}^{m}$, such that
$\bm{\phi}(t_{f},t_{0},0,\bm{u},q)=\bm{x}$;

  2. (ii)

_unobservable via $q$_, if there exists an input $\bm{u}$, such

$\bm{C}^{\mathrm{diff}}_{q}\bm{\phi}(t,t_{0},\bm{x},\bm{u},q)=\bm{C}^{\mathrm{diff}}_{q}\bm{\phi}(t,t_{0},\mathbf{0},\bm{u},q)\quad\text{for
all }t\geq t_{0}.$

The _reachable_ and _unobservable set_ via $q$, denoted by $\mathcal{R}_{q}$
and $\mathcal{U}\mathcal{O}_{q}$, respectively, is the set of states which are
reachable and unobservable via $q$, respectively. The _observable set_ via $q$
of (2.11), denoted by $\mathcal{O}_{q}$, is defined as
$\mathcal{O}_{q}\vcentcolon=(\mathcal{U}\mathcal{O}_{q})^{\perp}$. The set of
reachable states $\mathcal{R}$ and the set of observable states $\mathcal{O}$
of (2.11) can be defined as

$\displaystyle\mathcal{R}\vcentcolon=\bigcup_{q\in\mathcal{S}}\mathcal{R}_{q}\qquad\text{and}\qquad\mathcal{O}\vcentcolon=\bigcup_{q\in\mathcal{S}}\mathcal{O}_{q}.$
(2.12)

######  Lemma 2.9 (See [26, 20]).

For a given switching signal $q\in\mathcal{S}$ with $\kappa\in\mathbb{N}$
switches at time instants $t_{k}$ with $k=0,\ldots,\kappa$, define
$\tau_{k}\coloneqq t_{k+1}-t_{k}$, and consider the following recursive
relation

$\displaystyle\mathcal{M}_{0}$ $\displaystyle\coloneqq\mathcal{R}_{q_{0}},$
$\displaystyle\mathcal{M}_{k}$
$\displaystyle\coloneqq\mathcal{R}_{q_{k}}+\operatorname{e}^{\bm{A}_{q_{k}}^{\mathrm{diff}}\tau_{k}}\mathbf{\Pi}_{q_{k}}\mathcal{M}_{k-1},$
$\displaystyle k$ $\displaystyle=1,\ldots,\kappa,$
$\displaystyle\mathcal{N}_{\kappa}$
$\displaystyle\coloneqq\mathcal{U}\mathcal{O}_{q_{\kappa}},$
$\displaystyle\mathcal{N}_{k}$
$\displaystyle\coloneqq\mathcal{U}\mathcal{O}_{q_{k}}\cap\left(\operatorname{e}^{-\bm{A}_{q_{k}}^{\mathrm{diff}}\tau_{k}}\mathbf{\Pi}^{-1}_{q_{k+1}}\mathcal{N}_{q_{k+1}}\right),$
$\displaystyle k$ $\displaystyle=\kappa-1,\ldots,0,$

where
$\mathcal{R}_{q_{k}}\vcentcolon=\big{\langle}\bm{A}^{\mathrm{diff}}_{q_{k}}\mid\operatorname{img}(\bm{B}^{\mathrm{diff}}_{q_{k}})\big{\rangle}$
and
$\mathcal{U}\mathcal{O}_{q_{k}}\vcentcolon=\big{\langle}\operatorname{ker}(\bm{C}^{\mathrm{diff}}_{q_{k}})\mid\bm{A}^{\mathrm{diff}}_{q_{k}}\big{\rangle}$.
Then the reachable and observable set via $q$ satisfy

$\mathcal{R}_{q}=\mathcal{M}_{\kappa}\qquad\text{and}\qquad\mathcal{O}_{q}=\left(\mathcal{U}\mathcal{O}_{q}\right)^{\perp}\;=\;\mathcal{N}_{0}^{\perp},$

respectively.

##  3\. Model Reduction

In this section, we present a MOR method for the switched DAE (1.1) starting
from the simplified switched ODE with state jumps (2.11). Note that, the
additional feedthrough term $\bm{D}_{k}\bm{U}(t)$ does not need to be
approximated, since it does not depend on the high-dimensional state, and can
be exactly transferred to the reduced system. The same holds for the
additional impulsive output $\bm{U}_{k}^{\mathrm{imp}}$. Finally, we discuss
how the presented MOR framework can be adapted to deal with the more general
class of switched DAEs represented by (2.9).

###  3.1. Model reduction algorithm

Our idea consists of applying the balancing MOR method proposed in [39] for
switched ODEs, where the symmetric positive semi-definite Gramians are
computed from the _generalized Lyapunov equations_ (GLEs)

$\displaystyle\mathbf{\mathscr{A}}\mathbf{\mathcal{P}}+\mathbf{\mathcal{P}}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\mathbf{\mathcal{P}}\mathbf{\mathscr{F}}_{j}^{\mathsf{T}}+\bm{B}^{\mathrm{diff}}_{j}\left(\bm{B}^{\mathrm{diff}}_{j}\right)^{\mathsf{T}}\right)$
$\displaystyle=\mathbf{0},$ (3.1a)
$\displaystyle\mathbf{\mathscr{A}}^{\mathsf{T}}\mathbf{\mathcal{Q}}+\mathbf{\mathcal{Q}}\mathbf{\mathscr{A}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}^{\mathsf{T}}_{j}\mathbf{\mathcal{Q}}\mathbf{\mathscr{F}}_{j}+\left(\bm{C}^{\mathrm{diff}}_{j}\right)^{\mathsf{T}}\bm{C}^{\mathrm{diff}}_{j}\right)$
$\displaystyle=\mathbf{0},$ (3.1b)

where $\mathbf{\mathscr{A}}\vcentcolon=\bm{A}^{\mathrm{diff}}_{1}$ and
$\mathbf{\mathscr{F}}_{j}\vcentcolon=\bm{A}^{\mathrm{diff}}_{j}-\bm{A}^{\mathrm{diff}}_{1}$
for $j=1,\ldots,M$. Let $\mathbf{\mathcal{P}}=\bm{S}\bm{S}^{\mathsf{T}}$ and
$\mathbf{\mathcal{Q}}=\bm{R}\bm{R}^{\mathsf{T}}$ denote the Cholesky
decomposition of the Gramians and construct the projection matrices $\bm{V}$
and $\bm{W}$ for the ROM (1.4) using the _singular value decomposition_ (SVD)
of the product of the Cholesky factors, i.e.,

$\bm{V}=\bm{S}\bm{U}_{1}\mathbf{\Sigma}^{-1/2}_{1}\quad\text{and}\quad\bm{W}=\bm{R}\bm{V}_{1}\mathbf{\Sigma}^{-1/2}_{1}$
(3.2)

with

$\bm{\mathcal{H}}\coloneqq\bm{S}^{\mathsf{T}}\bm{R}=\left[\bm{U}_{1},\bm{U}_{2}\right]\begin{bmatrix}\mathbf{\Sigma}_{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{\Sigma}_{2}\end{bmatrix}[\bm{V}_{1},\bm{V}_{2}]^{\mathsf{T}}.$
(3.3)

where $\bm{\mathcal{H}}$ is the Hankel matrix. This procedure is denoted as
square-root balanced truncation (see [1, Sec. 7.3]). Under certain conditions,
see [39, Ass. 1, Prop. 1], the ROM (1.3) obtained through the Petrov-Galerkin
projection by (3.2) is quadratically stable [38]. If the Gramians computed in
(3.1) satisfy a set of _linear matrix inequalities_ (LMIs) of [38, Def. 10],
then one can show, see [39, Thm. 4] and [38, Thm. 6], that the error between
the output of the FOM and the one of the ROM is bounded by

$\|\bm{y}-\hat{\bm{y}}\|_{L_{2}}\leq
2\\!\\!\sum_{i=r+1}^{n}\sigma_{i}\|\bm{u}\|_{L_{2}},$ (3.4)

for every switching signal $q$, where $\sigma_{i}$ are the neglected singular
values of the Hankel matrix; see [58] for a relation to the Kolmogorov
$n$-widths.

Two issues prevent us from straightforwardly using this method. First, it
requires a switched system with asymptotically stable matrices, while, by
(2.7a), it is clear that $\bm{A}^{\mathrm{diff}}_{j}$ may be singular. Second,
the fact that $\mathcal{P}$ and $\mathcal{Q}$ are suitable in representing,
respectively, the reachable and observable states was proven for a switched
system without state jumps at the switching times; see [39, Thm. 3].
Nevertheless, since the state transition matrices are projections induced by
the active DAE mode, we can show that the transition matrices can be neglected
in the computation of the reachable and observable sets in Lemma 2.9.

######  Theorem 3.1.

Consider the system (2.11) and let $\mathcal{R}$ and $\mathcal{O}$ denote the
reachable and observable states. Moreover, let $\tilde{\mathcal{R}}$ and
$\tilde{\mathcal{O}}$ denote the reachable and observable states for (2.11)
with artificially modified $\mathbf{\Pi}_{k}=\bm{I}$ for all
$k\in\mathcal{J}$. Then

$\displaystyle\mathcal{R}=\tilde{\mathcal{R}}\qquad\text{and}\qquad\mathcal{O}=\tilde{\mathcal{O}}.$
(3.5)

###### Proof.

Direct calculations using the definitions (2.6) and (2.7) yield
$\big{(}\bm{A}_{k}^{\mathrm{diff}}\big{)}^{j}\mathbf{\Pi}_{k}=\big{(}\bm{A}_{k}^{\mathrm{diff}}\big{)}^{j}$
for any $j\in\mathbb{N}_{0}$ and $k\in\mathcal{J}$. Thus, the first identity
follows immediately from Lemma 2.9 and (2.12). For the second identity, we
show inductively
$\mathbf{\Pi}_{q_{k+1}}^{-1}\mathcal{N}_{q_{k+1}}=\mathcal{N}_{q_{k+1}}$ by
exploiting
$\bm{A}_{k}^{\mathrm{diff}}\mathbf{\Pi}_{k}=\bm{A}_{k}^{\mathrm{diff}}$, which
implies

$\bm{C}_{k}^{\mathrm{diff}}A_{k}^{\mathrm{diff}}\bm{x}=0\quad\iff\quad\bm{C}_{k}^{\mathrm{diff}}A_{k}^{\mathrm{diff}}\mathbf{\Pi}_{k}\bm{x}=0$

for any $\bm{x}\in\mathbb{R}^{n}$, which together with Lemma 2.9 and (2.12)
completes the proof. ∎

It remains to address the required asymptotic stability, which we do by a
suitable reformulation that we propose in the next subsection, which is then
used to set up the GLEs from which we recover the Gramians for the switched
system (2.11).

###  3.2. Reformulation of the switched system

Let us observe that $\bm{A}_{q_{k}}^{\mathrm{diff}}$ is singular because the
switched system (2.11) has the state variable expressed in the full dimension
$n$ but it describes the state dynamic using only the differential part of
each DAE subsystem. Here we provide a reformulation of (2.11) that preserves
input-output mapping and thus does not alter reachability and observability
sets, and, at the same time, ensures to have asymptotically stable matrices
for each subsystem by making use of the appropriate number of degrees of
freedom for the state variable of each system mode.

We start by writing down the following switched system with switching state
coordinates

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\hat{\bm{z}}}_{q_{k}}(t)&=\begin{bmatrix}\bm{J}_{q_{k}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\hat{\bm{z}}_{q_{k}}(t)+\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k}}}}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\bm{S}_{q_{k}}\bm{B}_{q_{k}}\bm{u}(t),&t\in(t_{k},t_{k+1}),\\\
\hat{\bm{z}}_{q_{k}}(t_{k}^{+})&=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{q_{k}}}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\bm{T}_{q_{k}}^{-1}\bm{T}_{q_{k-1}}\hat{\bm{z}}_{{q_{k-1}}}(t_{k}^{-}),&\bm{z}(t_{0}^{-})=\mathbf{0},\\\
\bm{z}(t)&=\bm{T}_{q_{k}}\hat{\bm{z}}_{q_{k}}(t),\\\
\bm{y}(t)&=\bm{C}_{q_{k}}\bm{T}_{q_{k}}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k}}}}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\hat{\bm{z}}_{q_{k}}(t)+\bm{D}_{q_{k}}\bm{U}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{y}[t_{k}]&=\bm{U}_{q_{k}}^{\mathrm{imp}}(t_{k}).\end{aligned}\right.$
(3.6)

Since (3.6) is obtained via a state-space transformation, the input-output
mapping is not altered. We immediately observe that the last
$n-n_{\bm{J}_{q_{k}}}$ entries of $\hat{\bm{z}}_{q_{k}}(t)$ are always zero,
therefore (3.6) can be reduced to

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\tilde{\bm{z}}}_{q_{k}}(t)&=\bm{J}_{q_{k}}\tilde{\bm{z}}_{q_{k}}(t)+\begin{bmatrix}\bm{I}_{n_{\bm{J}_{q_{k}}}}&\mathbf{0}\end{bmatrix}\bm{S}_{q_{k}}\bm{B}_{q_{k}}\bm{u}(t),&t\in(t_{k},t_{k+1}),\\\
\tilde{\bm{z}}_{q_{k}}(t_{k}^{+})&=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{q_{k}}}}&\mathbf{0}\end{bmatrix}\bm{T}_{q_{k}}^{-1}\bm{T}_{q_{k-1}}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{q_{k-1}}}}&\mathbf{0}\end{bmatrix}^{\mathsf{T}}\tilde{\bm{z}}_{q_{k-1}}(t_{k}^{-}),&\bm{z}(t_{0}^{-})=\mathbf{0},\\\
\tilde{\bm{z}}_{q_{k}}(t)&=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k}}}}}&\mathbf{0}\end{bmatrix}\hat{\bm{z}}_{q_{k}}(t),\\\
\bm{y}(t)&=\bm{C}_{q_{k}}\bm{T}_{q_{k}}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k}}}}}&\mathbf{0}\end{bmatrix}^{\mathsf{T}}\tilde{\bm{z}}_{q_{k}}(t)+\bm{D}_{q_{k}}\bm{U}(t),&t\in(t_{k},t_{k+1}),\\\
\bm{y}[t_{k}]&=\bm{U}_{q_{k}}^{\mathrm{imp}}(t_{k})\end{aligned}\right.$ (3.7)

with rectangular state transition matrices

$\bm{P}_{{q_{k}},{q_{k-1}}}\vcentcolon=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k}}}}}&\mathbf{0}\end{bmatrix}\bm{T}_{{q_{k}}}^{-1}\bm{T}_{{q_{k-1}}}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{{q_{k-1}}}}}&\mathbf{0}\end{bmatrix}^{\mathsf{T}}\;\in\;\bm{R}^{n_{\bm{J}_{{q_{k}}}}\times
n_{\bm{J}_{q_{k-1}}}}.$

With these preparations, we can now rewrite system (3.7) as a bilinear-system
following the strategy presented in [39, Section II.$A$]. We omit the details
and define the resulting matrices for the GLEs (3.1)

$\displaystyle\begin{aligned}
\mathbf{\mathscr{A}}&\vcentcolon=\bm{J}_{1},&\mathbf{\mathscr{F}}_{j}&\vcentcolon=\bm{P}_{1,j}\bm{J}_{j}\bm{P}_{j,1}-\bm{J}_{1},\\\
\mathbf{\mathscr{B}}_{j}&\vcentcolon=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{j}}}&\mathbf{0}_{n_{\bm{J}_{1}}-n_{\bm{J}_{j}}}\end{bmatrix}^{\mathsf{T}}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{j}}}&\mathbf{0}\end{bmatrix}\bm{S}_{j}\bm{B}_{j},&\mathbf{\mathscr{C}}_{j}&\vcentcolon=\bm{C}_{j}\bm{T}_{j}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{j}}}&\mathbf{0}\end{bmatrix}^{\mathsf{T}}\bm{P}_{j,1}.\end{aligned}$
(3.8)

for $j=1,\ldots,M$ and with $\bm{J}_{1}$ associated to the system mode with
the largest differential dimension $n_{\bm{J}_{k}}$.

The construction of the matrices $\mathbf{\mathcal{P}}$ and
$\mathbf{\mathcal{Q}}$ can now be justified as follows: Provided 2.7 is
satisfied, Theorem 2.6 and the modifications in this subsection guarantee that
(1.1) and (3.7) have the same input-output mapping. For system (3.7), we have
shown with Theorem 3.1 that the reachable and observable states can be
computed by ignoring the consistency projectors, which motivates us to solve
the GLEs, whose solutions encode the reachable and observable states [39, Thm.
$3.2$].

###  3.3. Handling the general class of switched DAEs

In this subsection we discuss how to relax 2.7. First, we 2.7 (i). We observe
that for each input $\bm{u}$, the vectors $\bm{U}_{q_{k}}^{-}(t_{k}^{-})$
defined in (2.10) satisfy $\bm{U}_{q_{k}}^{-}\in\mathcal{R}^{\mathrm{imp}}$
with

$\displaystyle\mathcal{R}^{\mathrm{imp}}\vcentcolon=\bigcup_{j=1}^{M}\left(\bigcup_{i=0}^{\nu_{j}-1}\operatorname{span}\left(\left(\bm{E}_{q_{j}}^{\mathrm{imp}}\right)^{i}\bm{B}_{q_{j}}^{\mathrm{imp}}\right)\right).$
(3.9)

Note that the jumps induced by $\bm{U}_{q_{k}}^{-}$ are applied to the state
via the projection matrix $\mathbf{\Pi}_{q_{k}}$ such that we immediately
obtain
$\mathbf{\Pi}_{q_{k}}\mathcal{R}^{\mathrm{imp}}\subseteq\mathcal{R}^{\mathrm{imp}}$.
Clearly, the jumps influence the reachability set of the switched system, and
one possible way to include them in the reachability set of the switched
system is to introduce an additional term in the GLE associated with the
reachability Gramian, i.e.,

$\displaystyle\mathbf{\mathscr{A}}\mathbf{\mathcal{P}}+\mathbf{\mathcal{P}}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\mathbf{\mathcal{P}}\mathbf{\mathscr{F}}_{j}^{\mathsf{T}}+\mathbf{\mathcal{B}}_{j}\mathbf{\mathcal{B}}_{j}^{\mathsf{T}}+\mathbf{\mathcal{B}}^{\mathrm{imp}}_{j}(\mathbf{\mathcal{B}}_{j}^{\mathrm{imp}})^{\mathsf{T}}\right)=\mathbf{0},$
(3.10)
$\displaystyle\mathbf{\mathcal{B}}^{\mathrm{imp}}_{j}\vcentcolon=\begin{bmatrix}\bm{I}_{n_{\bm{J}_{1}}}&\mathbf{0}\end{bmatrix}\bm{T}_{1}^{-1}\begin{bmatrix}(\bm{E}_{j}^{\mathrm{imp}})^{0},\ldots,(\bm{E}_{j}^{\mathrm{imp}})^{\nu_{j}-1}\end{bmatrix}\bm{B}_{j}^{\mathrm{imp}}.$
(3.11)

The rigorous proof for this is beyond the purpose of our work but ideally
could be formulated combining the proof of [39, Thm. 3.2], the precise
characterization of the input-dependent jumps (3.9), and the definition and
characterization of reachable states for a given switched signal (see Lemma
2.9).

Now, we consider the case of 2.7 (ii) not satisfied. We notice that, even in
this case, one may still not be interested in considering the state-dependent
output at the switching times. Indeed, due to the presence of the Dirac delta,
the output at the switching times is meaningful when evaluated inside an
integral operator. Otherwise, the Dirac delta impulse can only be interpreted
as instantaneous infinite energy in response to the state jumps. However, if
interested in the approximation of this term, we can adopt the idea for the
input-dependent jumps and replace the GLE for the observability Gramian with
the modified GLE

$\displaystyle\mathbf{\mathscr{A}}^{\mathsf{T}}\mathbf{\mathcal{Q}}+\mathbf{\mathcal{Q}}\mathbf{\mathscr{A}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}^{\mathsf{T}}_{j}\mathbf{\mathcal{Q}}\mathbf{\mathscr{F}}_{j}+\mathbf{\mathcal{C}}_{j}^{\mathsf{T}}\mathbf{\mathcal{C}}_{j}+\left(\mathbf{\mathcal{C}}^{\mathrm{imp}}_{j}\right)^{\mathsf{T}}\mathbf{\mathcal{C}}^{\mathrm{imp}}_{j}\right)=\mathbf{0},$
$\displaystyle\mathbf{\mathcal{C}}^{\mathrm{imp}}_{j}\vcentcolon=\bm{C}^{\mathrm{imp}}_{j}\begin{bmatrix}\bm{E}_{j}^{\mathrm{imp}},\ldots,(\bm{E}_{j}^{\mathrm{imp}})^{\nu_{j}-1}\end{bmatrix}^{\mathsf{T}}\bm{T}_{j}\begin{bmatrix}\bm{I}_{n_{\bm{J}_{j}}}&\mathbf{0}\end{bmatrix}^{\mathsf{T}}\bm{P}_{j,1}.$

##  4\. Efficient and certified numerics for large scale systems

For the large-scale setting it is crucial to exploit sparsity and state-of-
the-art methods in every step of the MOR algorithm presented in Section 3.
Besides computational efficiency, we further aim for an error certification
and thus have to balance the errors of the different methods. In more detail,
we discuss the efficient computation of the decoupling matrices for the QWF
(2.2) in Section 4.1. In Section 4.2 we present the algorithm used to
approximate the solutions of (3.1) and we derive efficiently computable error
estimates to certify the quality of the approximation. Finally, in Section 4.3
we show how to combine the error estimates derived in the previous section
with the error estimate for the MOR of switched linear system through balance
truncation.

###  4.1. Efficient computation of the Wong sequences space

The first numerical issue concerns the determination of the two subspaces
$\mathcal{V}^{\star}$ and $\mathcal{W}^{\star}$ from (2.4) and two
representative matrices, $\hat{\mathbf{V}}$ and $\hat{\mathbf{W}}$ such that
$\operatorname{img}{(\hat{\mathbf{V}})}=\mathcal{V}^{\star}$ and
$\operatorname{img}{(\hat{\mathbf{W}})}=\mathcal{W}^{\star}$, necessary to
transform the DAE (2.1) in the QWF; see Theorem 2.1. In essence, an algorithm
determining $\mathcal{V}^{\star}$ and $\mathcal{W}^{\star}$ for large systems,
should retain the following two features:

  1. (i)

if the matrix pair $(\bm{E},\bm{A})$ is sparse, then forming dense matrices of
size $n\times n$ to determine $\hat{\mathbf{V}}$, $\hat{\mathbf{W}}$, and the
matrices in (2.7) should be avoided;

  2. (ii)

the computation of $\hat{\mathbf{V}}$ and $\hat{\mathbf{W}}$ should be
efficient, i.e., involving at the most $\mathcal{O}(n)$ floating point
operations for sparse matrices and $\mathcal{O}(n^{2})$ for full matrices.

The starting point is (2.3) and the observation that the sequence converges
after $\nu$ steps, see [6, Prop. 2.10], where $\nu$ is the nilpotent index of
the DAE. The kernel of $\bm{E}$ has a crucial role, indeed, from (2.3), one
can easily check that

$\mathcal{W}^{1}=\ker(\bm{E}),\quad\mathcal{V}^{1}=\bm{A}^{-1}\bm{E}\ker(\bm{E})^{\perp}.$
(4.1)

We emphasize that in many applications, see for instance [32, Sec. 5], $\nu$
is typically known and is less or equal three, i.e., the limiting spaces
$\mathcal{V}^{\star}$ and $\mathcal{W}^{\star}$ in (2.4) are determined after
few iterations. Moreover, the systems arise in a highly structured form that
either simplify the computation or can be used to directly read of the desired
matrices $\hat{\mathbf{V}}$ and $\hat{\mathbf{W}}$ in a sparse manner.
Therefore, it is often the case that a sparse matrix representation with
orthogonal columns of $\ker(\bm{E})$, say $\bm{K}$, can be straightforwardly
derived. For instance, this is the case when

$\bm{E}=\begin{bmatrix}\bm{M}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}$ (4.2)

with $\bm{M}$ a square invertible matrix. Then, we can distinguish between two
cases. First, when $\dim(\ker(\bm{E}))$ is independent from $n$, like in the
mass-spring system of Section 5.1, and second, when $\dim(\ker(\bm{E}))$
increases proportionally with respect to $n$, like in the Stokes control
system of Section 5.2. Without going into the details, it is possible to show
that, whenever a sparse orthogonal representation of $\ker(\bm{E})$ and of its
orthogonal space is available, then requirement (i) can be always satisfied
avoiding forming explicitly the pair $(\hat{\mathbf{V}},\hat{\mathbf{W}})$
while (ii) is generally only satisfied when the dimension of $\ker(\bm{E})$ is
independent from $n$.

###  4.2. Error certification for numerical approximation of GLEs

Once the matrices $\hat{\mathbf{V}}$ and $\hat{\mathbf{W}}$ are available,
either implicitly or explicitly, they can be used to transform each mode of
the switched DAE (1.1) to QWF according to Theorem 2.3 and then identifying
the matrices defined in (2.7). Our next task is thus to solve the GLEs (3.1)
associated with the switched linear system with jumps (3.7).

Algorithm 1 Stationary iterations for the GLEs (4.4)

Input: Matrices $\mathbf{\mathscr{A}}$, $\mathbf{\mathscr{F}}_{j}$,
$\mathbf{\mathscr{B}}_{j}^{\mathrm{diff}}$ for $j=1,\ldots,M$

Output: $\tilde{\bm{Z}}$ such that
$\bm{X}\approx\tilde{\bm{Z}}\tilde{\bm{Z}}^{\mathsf{T}}$ is an approximation
to the solution of (4.4)

1:Set
$\tilde{\bm{B}}_{1}\;\vcentcolon=\;[\mathbf{\mathscr{B}}_{1},\ldots,\mathbf{\mathscr{B}}_{M}]$

2:Approximately solve
$\mathbf{\mathscr{A}}\bm{X}+\bm{X}\mathbf{\mathscr{A}}^{\mathsf{T}}+\tilde{\bm{B}}_{1}\tilde{\bm{B}}_{1}^{\mathsf{T}}=\mathbf{0}$
for $\bm{X}_{1}=\bm{Z}_{1}\bm{Z}_{1}^{\top}$

3:for $k=2,3,\ldots$ do

4: Set
$\tilde{\bm{B}}_{k}\vcentcolon=[\mathbf{\mathscr{F}}_{1}\bm{Z}_{k-1},\ldots,\mathbf{\mathscr{F}}_{M}\bm{Z}_{k-1},\tilde{\bm{B}}_{1}]$

5: Approximately solve
$\mathbf{\mathscr{A}}\bm{X}+\bm{X}\mathbf{\mathscr{A}}^{\mathsf{T}}+\tilde{\bm{B}}_{k}\tilde{\bm{B}}^{\mathsf{T}}_{k}=\mathbf{0}$
(4.3) forfor $\bm{X}_{k}=\bm{Z}_{k}\bm{Z}_{k}^{\mathsf{T}}$

6: if sufficiently accurate then stop

7:end for

For simplicity we focus on the approximation of the reachability Gramian
$\mathbf{\mathcal{P}}$ through the GLE

$\mathbf{\mathscr{A}}\bm{X}+\bm{X}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\bm{X}\mathbf{\mathscr{F}}_{j}^{\mathsf{T}}\right)+\sum_{j=1}^{M}\left(\mathbf{\mathscr{B}}_{j}\mathbf{\mathscr{B}}_{j}^{\mathsf{T}}\right)=\mathbf{0},\\\
$ (4.4)

where $\mathbf{\mathscr{A}}$, $\mathbf{\mathscr{F}}_{j}$, and
$\mathbf{\mathscr{B}}_{j}$ are as defined in (3.8). The approximate solution
of the controllability Gramian $\mathbf{\mathcal{Q}}$ follows analogously by
replacing $\mathbf{\mathscr{A}}$, $\mathbf{\mathscr{F}}_{j}$, and
$\mathbf{\mathscr{B}}_{j}$ with $\mathbf{\mathscr{A}}^{\mathsf{T}}$,
$\mathbf{\mathscr{F}}^{\mathsf{T}}_{j}$, and $\mathbf{\mathscr{C}}_{j}$,
respectively.

Existence and uniqueness of a solution of (4.4) can be studied via Kronecker
algebra, i.e., by introducing the matrices

$\displaystyle\mathcal{L}\vcentcolon=\bm{I}_{n}\otimes\mathbf{\mathscr{A}}+\mathbf{\mathscr{A}}\otimes\bm{I}_{n}\in\mathbb{R}^{n^{2}\times
n^{2}}\qquad\text{and}\qquad\Pi\vcentcolon=\sum_{j=1}^{M}\mathbf{\mathscr{F}}_{j}\otimes\mathbf{\mathscr{F}}_{j}\in\mathbb{R}^{n^{2}\times
n^{2}},$ (4.5)

such that the vectorized form of (4.4) is given by
$(\mathcal{L}+\Pi)\operatorname{vec}2(\bm{X})=-\operatorname{vec}2(\bm{B}\bm{B}^{\mathsf{T}})$.
Hence, the GLE (4.4) is uniquely solvable if and only if $\mathcal{L}+\Pi$ is
nonsingular.

######  Assumption 4.1.

The GLE (4.4) is uniquely solvable.

We rely on the stationary iteration presented in Algorithm 1 taken from [49,
Alg. 2.1] to solve the GLE (4.4) numerically. It consists of a fixed-point
iteration where in each iteration we need to solve a classical Lyapunov
equation. To ensure that Algorithm 1 is feasible and convergent, we make the
following assumption; see for instance [23].

######  Assumption 4.2.

The matrix $\mathbf{\mathscr{A}}\in\mathbb{R}^{n\times n}$ in (4.4) is
asymptotically stable, i.e., all its eigenvalues are contained in the open
left-half complex plane, and $\|\mathcal{L}^{-1}\Pi\|_{2}<1$.

######  Remark 4.3.

On a first glance, the assumption $\|\mathcal{L}^{-1}\Pi\|_{2}<1$ sounds
rather restrictive. Nevertheless, one can always proceed by defining the
rescaled matrix

$\tilde{\Pi}=\frac{1}{\|\mathcal{L}^{-1}\Pi\|_{2}+\delta}\Pi,$

for some $\delta>0$ and observing that rescaling $\Pi$ does not affect the
space spanned by the columns of the solution matrix $\bm{X}$; see [10, 5].

An efficient implementation of Algorithm 1 passes through the solution of the
following numerical issues:

  1. (i)

At each iteration the Lyapunov equation requires an efficient solver.

  2. (ii)

An opportune stopping criteria is required.

In the following subsections, we describe how to address these points.

####  4.2.1. Approximate solution of Lyapunov equations involving large sparse
matrices

Consider (4.3) with large and sparse $\mathbf{\mathscr{A}}$. The literature
for the case of large and sparse matrices is vast; see, for example, the
survey [51] about methods for matrix equations and references therein. Almost
all methods in the large-scale setting have in common that they construct a
sequence of approximations converging to the true solution. Thus, a suitable
stopping criterion is required to balance the approximation error in the
computation of the solution of the Lyapunov equation and the solution of the
stationary iteration from Algorithm 1. To this end, assume that
$\bm{X}_{\ell}\in\mathcal{S}^{n}_{\succeq}$ is the approximation of the
solution $\bm{X}\in\mathcal{S}^{n}_{\succeq}$ of the Lyapunov equation (4.3)
in the $\ell$th step (and $k$th iteration of Algorithm 1) and let

$\bm{R}_{\ell}\vcentcolon=\mathbf{\mathscr{A}}\bm{X}_{\ell}+\bm{X}_{\ell}\mathbf{\mathscr{A}}^{\mathsf{T}}+\tilde{\bm{B}}_{k}\tilde{\bm{B}}_{k}^{\mathsf{T}}$
(4.6)

denote the corresponding residual. We obtain the following error-residual
relation.

######  Proposition 4.4 (Error bound for approximate solution of Lyapunov
equations).

Consider the Lyapunov equation (4.3) with asymptotically stable
$\mathbf{\mathscr{A}}$ and assume that
$\bm{X}_{\ell}\in\mathcal{S}^{n}_{\succeq}$ is an approximation of the unique
solution $\bm{X}\in\mathcal{S}^{n}_{\succeq}$ of (4.3). Let
$\bm{E}_{\ell}\vcentcolon=\bm{X}_{\ell}-\bm{X}$. Then

$\|\bm{E}_{\ell}\|_{2}\leq\frac{\|\bm{R}_{\ell}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})},$
(4.7)

where $\sigma_{\min}(\mathbf{\mathscr{A}})$ denotes the smallest singular
value of the matrix $\mathbf{\mathscr{A}}$ and $\bm{R}_{\ell}$ is as defined
in (4.6).

###### Proof.

Substituting (4.6) into the residual equation (4.6) yields the error-residual
relation

$\displaystyle\bm{R}_{\ell}=\mathbf{\mathscr{A}}\bm{E}_{\ell}+\bm{E}_{\ell}\mathbf{\mathscr{A}}^{\mathsf{T}},$
(4.8)

which together with the symmetry of $\bm{R}_{\ell}$ yields that
$\bm{E}_{\ell}$ is symmetric; see, e.g., [27, Cha. 12.3, Thm. 3]. Using
Kronecker algebra and the notation introduced in (4.5) we can reformulate
(4.8) in vectorial from as
$\operatorname{vec}2{(\bm{R}_{\ell})}=\mathcal{L}\operatorname{vec}2{(\bm{E}_{\ell})}$.
Using $\|\bm{E}\|_{2}\leq\|\bm{E}\|_{\mathrm{F}}$ and
$\sigma_{\min}(\mathcal{L})=2\sigma_{\min}(\mathbf{\mathscr{A}})$ yields

$\displaystyle\|\bm{E}_{\ell}\|_{2}$
$\displaystyle\leq\|\bm{E}_{\ell}\|_{\mathrm{F}}=\|\operatorname{vec}2(\bm{E}_{\ell})\|_{2}=\|\mathcal{L}^{-1}\operatorname{vec}2{(\bm{R}_{\ell})}\|_{2}\leq\|\mathcal{L}^{-1}\|_{2}\|\operatorname{vec}2{(\bm{R}_{\ell})}\|_{2}$
$\displaystyle=\frac{\|\bm{R}_{\ell}\|_{\mathrm{F}}}{\sigma_{\min}(\mathcal{L})}=\frac{\|\bm{R}_{\ell}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})}.\qed$

We observe that the smallest singular value
$\sigma_{\min}(\mathbf{\mathscr{A}})$ needs to be computed only once. The
computational cost is $\mathcal{O}(n)$ floating point operations if
$\mathbf{\mathscr{A}}$ is sparse. In the large-scale case, on can for instance
rely one the method presented in [47]. The questions that remain to be
answered are the following: ($i$) the choice of the iterative method to
determine $\bm{X}_{\ell}$ and ($ii$) how to efficiently compute the residual
norm $\|\bm{R}_{\ell}\|_{\mathrm{F}}$. Since these questions are closely tied
together, we first briefly describe the method for computing $\bm{X}_{\ell}$.
As proposed in [51] we project the Lyapunov equation (4.6) onto a smaller
subspace and then use standard solvers for the small equation111See [51, Sec.
5.1] for a review on the solution of small-scale Lyapunov equations, thus
obtaining the coordinates of the approximation in the subspace.

More precisely, consider the matrix $\bm{V}_{\ell}\in\mathbb{R}^{n\times
n_{\ell}}$ whose $n_{\ell}$ columns form an orthonormal basis of the
appropriate subspace $\mathcal{S}_{\ell}$ to be used for the projection. Let
${\bm{X}}_{\ell}=\bm{V}_{\ell}\bm{Y}_{\ell}\bm{V}_{\ell}^{\mathsf{T}}$ be a
low-rank approximation to the solution of (4.6) with
$\bm{Y}_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell}}$. Imposing the Galerkin
condition $\bm{V}_{\ell}^{\mathsf{T}}\bm{R}_{\ell}\bm{V}_{\ell}=\mathbf{0}$ on
the residual, we obtain the lower-order projected Lyapunov equation

$\mathbf{\mathscr{A}}_{\ell}\bm{Y}_{\ell}+\bm{Y}_{\ell}\mathbf{\mathscr{A}}_{\ell}^{\mathsf{T}}+\tilde{\bm{B}}_{\ell}\tilde{\bm{B}}_{\ell}^{\mathsf{T}}=\mathbf{0},$
(4.9)

with
$\mathbf{\mathscr{A}}_{\ell}\vcentcolon=\bm{V}_{\ell}^{\mathsf{T}}\mathbf{\mathscr{A}}\bm{V}_{\ell}$
and
$\tilde{\bm{B}}_{\ell}\vcentcolon=\bm{V}^{\mathsf{T}}_{\ell}\tilde{\bm{B}}_{k}$.
For $\mathcal{S}_{\ell}$ one can chose, for instance, the proposal from [40]
and use the block Krylov subspace

$\mathcal{K}_{\ell}(\mathbf{\mathscr{A}},\tilde{\bm{B}}_{k})=\operatorname{span}\\{\tilde{\bm{B}}_{k},\mathbf{\mathscr{A}}\tilde{\bm{B}}_{k},\mathbf{\mathscr{A}}^{2}\tilde{\bm{B}}_{k},\ldots,\mathbf{\mathscr{A}}^{\ell-1}\tilde{\bm{B}}_{k}\\}.$
(4.10)

For more sophisticated techniques such as rational Krylov subspaces we refer
to [50, 51].

Besides computational efficiency, the low-rank approximation can also be
exploited for an efficient evaluation of the Frobenius norm of the residual.
Indeed, using the Krylov subspace (4.10), one can show
$\|\bm{R}_{\ell}\|=\|\bm{V}_{\ell+1}^{\mathsf{T}}\bm{R}_{\ell}\bm{V}_{\ell+1}\|$
for the spectral and the Frobenius norm, see [50, Prop. 3.3]. This result is
particularly convenient since, for
$\bm{V}_{\ell+1}\vcentcolon=[\bm{V}_{\ell},\;\tilde{\bm{V}}]$, we have

$\displaystyle\bm{V}^{\mathsf{T}}_{\ell+1}\bm{R}_{\ell}\bm{V}_{\ell+1}=\begin{bmatrix}\bm{V}_{\ell}^{\mathsf{T}}\bm{R}_{\ell}\bm{V}_{\ell}&\bm{V}_{\ell}^{\mathsf{T}}\bm{R}_{\ell}\tilde{\bm{V}}\\\
\tilde{\bm{V}}^{\mathsf{T}}\bm{R}_{\ell}\bm{V}_{\ell}&\tilde{\bm{V}}^{\mathsf{T}}\bm{R}_{\ell}\tilde{\bm{V}}\end{bmatrix}=\begin{bmatrix}\mathbf{0}&\bm{Y}_{\ell}\bm{V}_{\ell}^{\mathsf{T}}\mathbf{\mathscr{A}}^{\mathsf{T}}\tilde{\bm{V}}\\\
\tilde{\bm{V}}^{\mathsf{T}}\mathbf{\mathscr{A}}\bm{V}_{\ell}\bm{Y}_{\ell}&\mathbf{0}\end{bmatrix}$
(4.11)

where we used expression (4.6), the fact that
${\bm{X}}_{\ell}=\bm{V}_{\ell}\bm{Y}_{\ell}\bm{V}_{\ell}^{\mathsf{T}}$, and
that $\tilde{\bm{V}}$ is by construction orthogonal to $\bm{V}_{\ell}$ and
$\tilde{\bm{B}}_{k}$. From (4.11) it directly follow that

$\|\bm{V}^{\mathsf{T}}_{\ell+1}\bm{R}_{\ell}\bm{V}_{\ell+1}\|_{\mathrm{F}}=\;2\|\tilde{\bm{V}}^{\mathsf{T}}\mathbf{\mathscr{A}}\bm{V}_{\ell}\bm{Y}_{\ell}\|_{\mathrm{F}}$

and the right hand side can be used to evaluate $\|\bm{R}\|_{\mathrm{F}}$ in a
computational efficient way.

####  4.2.2. Stopping criteria for Algorithm 1

Consider the approximate solution $\bm{X}_{k}$ for the GLE (4.4) provided by
Algorithm 1 at step $k$ and let us assume for the moment that the Lyapunov
equation in 5 is solved exactly. To understand if such a solution is accurate
enough for our requirements we need to estimate the norm of
$\bm{E}_{k}\vcentcolon=\bm{X}-\bm{X}_{k}$.

######  Proposition 4.5 (Stopping criterion for Algorithm 1).

Suppose 4.1 and 4.2 are satisfied and consider the error
$\bm{E}_{k}\vcentcolon=\bm{X}-\bm{X}_{k}$, where $\bm{X}$ is the unique
solution of the GLE (4.4) and $\bm{X}_{k}=\bm{Z}_{k}\bm{Z}_{k}^{\mathsf{T}}$
is the unique solution of the Lyapunov equation (4.3) in iteration $k$ of
Algorithm 1. Then Algorithm 1 converges and

$\|\bm{E}_{k}\|_{2}\leq\gamma\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}},$ (4.12)

with

$\gamma\vcentcolon=\frac{\|\mathcal{L}^{-1}\Pi\|}{1-\|\mathcal{L}^{-1}\Pi\|}.$
(4.13)

###### Proof.

Due to Algorithm 1 we obtain

$\mathbf{\mathscr{A}}\bm{X}_{k}+\bm{X}_{k}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\bm{X}_{k-1}\mathbf{\mathscr{F}}_{j}^{\mathsf{T}}\right)+\tilde{\bm{B}}_{1}\tilde{\bm{B}}_{1}^{\mathsf{T}}=\mathbf{0}$
(4.14)

for all $k\geq 2$. Then, from (4.4) we have that
$\tilde{\bm{B}}_{1}\tilde{\bm{B}}_{1}^{\mathsf{T}}=-\mathbf{\mathscr{A}}\bm{X}-\bm{X}\mathbf{\mathscr{A}}^{T}-\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\bm{X}\mathbf{\mathscr{F}}_{j}^{\mathsf{T}}\right)$,
and substituting into (4.14) we get

$\mathbf{\mathscr{A}}\bm{E}_{k}+\bm{E}_{k}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\bm{E}_{k-1}\mathbf{\mathscr{F}}_{k}^{\mathsf{T}}\right)=\mathbf{0}.$
(4.15)

For convenience we write (4.15) in its vectorial form
$\mathcal{L}\operatorname{vec}2(\bm{E}_{k})=-\Pi\operatorname{vec}2(\bm{E}_{k-1})$
with $\mathcal{L}$ and $\Pi$ defined in (4.5), thus we have

$\operatorname{vec}2(\bm{E}_{k})=-\mathcal{L}^{-1}\Pi\operatorname{vec}2(\bm{E}_{k-1}).$
(4.16)

From relation (4.16) directly follows that, numerically, Algorithm 1 converges
if and only if $\rho(\mathcal{L}^{-1}\Pi)<1$ where $\rho(\mathcal{L}^{-1}\Pi)$
denotes the spectral radius of $\mathcal{L}^{-1}\Pi$, see also [11, 49].
Observe that $\bm{E}_{k-1}=\bm{E}_{k}+(\bm{X}_{k}-\bm{X}_{k-1})$, therefore,
substituting in (4.16) and considering the norm on both sides, we get

$\displaystyle\begin{aligned}
\|\bm{E}_{k}\|_{\mathrm{F}}\;=\;\|\operatorname{vec}2(\bm{E}_{k})\|_{2}\;=&\;\|\mathcal{L}^{-1}\Pi\operatorname{vec}2\left(\bm{E}_{k}+\bm{X}_{k}-\bm{X}_{k-1}\right)\|_{2}\\\
\leq\;&\;\|\mathcal{L}^{-1}\Pi\|_{2}{\Big{(}\|\bm{E}_{k}\|_{\mathrm{F}}+\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}\Big{)}}.\end{aligned}$

By 4.2 we have $\|\mathcal{L}^{-1}\Pi\|<1$ and hence

$\|\bm{E}_{k}\|_{2}\leq\|\bm{E}_{k}\|_{\mathrm{F}}\leq\frac{\|\mathcal{L}^{-1}\Pi\|_{2}}{1-\|\mathcal{L}^{-1}\Pi\|_{2}}\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}.\qed$

To use (4.12) as a stopping criteria for Algorithm 1 requires the following
clarification. For the computation or approximation of
$\|\mathcal{L}^{-1}\Pi\|_{2}$ we observe

$\displaystyle\|\mathcal{L}^{-1}\Pi\|_{2}$
$\displaystyle\leq\|\mathcal{L}^{-1}\|_{2}\|\Pi\|_{2}=\frac{1}{2\sigma_{\mathrm{min}}(\mathbf{\mathscr{A}})}\Bigg{\|}\sum_{j=1}^{M}\mathbf{\mathscr{F}}_{j}\otimes\mathbf{\mathscr{F}}_{j}\Bigg{\|}_{2}$
$\displaystyle\leq\frac{1}{2\sigma_{\text{min}}(\mathbf{\mathscr{A}})}\sum_{j=1}^{M}\Big{\|}\mathbf{\mathscr{F}}_{j}\otimes\mathbf{\mathscr{F}}_{j}\Big{\|}_{2}=\frac{1}{2\sigma_{\mathrm{min}}(\mathbf{\mathscr{A}})}\sum_{j=1}^{M}\sigma_{1}(\mathbf{\mathscr{F}}_{j})^{2},$
(4.17)

and the right-hand side of (4.2.2) is computable in $\mathcal{O}(n)$ floating
points operations if the involved matrices are sparse. This motivates for any
$\delta>0$ to perform the scaling

$\tilde{\mathbf{\mathscr{F}}}_{j}\vcentcolon=\frac{\mathbf{\mathscr{F}}_{j}}{\sqrt{\beta+\delta}},\quad\text{with}\quad\beta\coloneqq\frac{1}{2\sigma_{\text{min}}(\mathbf{\mathscr{A}})}\sum_{j=1}^{M}\sigma_{1}(\mathbf{\mathscr{F}}_{j})^{2},\quad
j=1,\ldots M.$ (4.18)

Hence, it is not restrictive to consider $\gamma\approx 1$. Indeed, it is
sufficient to scale (4.18) with $\delta=\beta$ to get $\gamma<1$ and thus we
can get rid of $\gamma$ in (4.12).

######  Remark 4.6.

It is also convenient to scale the known term, i.e., to run Algorithm 1 with
$\tilde{\mathbf{\mathscr{B}}}_{j}\;\vcentcolon=\frac{\mathbf{\mathscr{B}}_{j}}{\sqrt{\beta+\delta}}$.
Then, after convergence, the output $\tilde{\bm{Z}}$ is multiplied by the same
factor.

Last, we have to compute $\|\bm{X}_{j}-\bm{X}_{j-1}\|_{\mathrm{F}}$, which
would require $\mathcal{O}(n^{2})$ floating point operations. To overcome this
issue we employ the following result.

######  Proposition 4.7.

Let $\bm{X}_{k}=\bm{Z}_{k}\bm{Z}^{\mathsf{T}}_{k}$ with
$\bm{Z}_{k}\in\mathbb{R}^{n\times n_{k}}$ and
$\bm{X}_{k-1}=\bm{Z}_{k-1}\bm{Z}^{\mathsf{T}}_{k-1}$ with
$\bm{Z}_{k}\in\mathbb{R}^{n\times n_{k-1}}$ be the solutions at iteration $k$
and $k-1$ of Algorithm 1. Then

$\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}\leq\left(\|\bm{Z}_{k}\|_{\mathrm{F}}+\|\bm{Z}_{k-1}\|_{\mathrm{F}}\right)\|\bm{Z}_{k}-\bm{Z}_{k-1}\|_{\mathrm{F}}.$
(4.19)

###### Proof.

Without loss of generality assume $n_{k}=n_{k-1}$. By making use of the
triangular inequality it directly follows

$\displaystyle\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}$
$\displaystyle=\|\bm{Z}_{k}\bm{Z}_{k}^{\mathsf{T}}-\bm{Z}_{k}\bm{Z}_{j-1}^{\mathsf{T}}+\bm{Z}_{k}\bm{Z}_{k-1}^{\mathsf{T}}-\bm{Z}_{k-1}\bm{Z}_{k-1}^{\mathsf{T}}\|_{\mathrm{F}}$
$\displaystyle\leq\|\bm{Z}_{k}\|_{\mathrm{F}}\|\bm{Z}_{k}^{\mathsf{T}}-\bm{Z}_{k-1}^{\mathsf{T}}\|_{\mathrm{F}}+\|\bm{Z}_{k}-\bm{Z}_{k-1}\|_{\mathrm{F}}\|\bm{Z}_{k-1}^{\mathsf{T}}\|_{\mathrm{F}}$
$\displaystyle=\left(\|\bm{Z}_{k}\|_{\mathrm{F}}+\|\bm{Z}_{k-1}\|_{\mathrm{F}}\right)\|\bm{Z}_{k}-\bm{Z}_{k-1}\|_{\mathrm{F}}.\qed$

Replacing $\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}$ with the right hand side
of (4.19) has the advantage that we only need $\mathcal{O}(n)$ floating point
operations for its evaluation (assuming $n_{k},n_{k-1}\ll n$). With these
preparations we are now ready to formulate our main result for the approximate
solution of the GLE (4.4).

######  Theorem 4.8 (Error bound for the approximate solution of GLE).

Suppose 4.1 and 4.2 are satisfied. Let $\bm{X}=\bm{Z}\bm{Z}^{\mathsf{T}}$ with
$\bm{Z}\in\mathbb{R}^{n\times n}$ be the unique solution of the GLE (4.4) and
$\tilde{\bm{X}}_{k}=\tilde{\bm{Z}}_{k}\tilde{\bm{Z}}_{k}^{\mathsf{T}}$ with
$\bm{Z}\in\mathbb{R}^{n\times n_{k}}$ be the approximate solution of the GLE
computed at iteration $k$ of Algorithm 1 applying the subspace projection
method in 5. Then

$\displaystyle\|\bm{X}-\tilde{\bm{X}}_{k}\|_{2}\leq\gamma\left(\|\tilde{\bm{Z}}_{k}\|_{\mathrm{F}}+\|\tilde{\bm{Z}}_{k-1}\|_{\mathrm{F}}\right)\|\tilde{\bm{Z}}_{k}-\tilde{\bm{Z}}_{k-1}\|_{\mathrm{F}}+\frac{(1+\gamma)\|\bm{R}_{k}\|_{\mathrm{F}}+\gamma\|\bm{R}_{k-1}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})},$
(4.20)

where $\bm{R}_{k}$ and $\bm{R}_{k-1}$ are the residuals at the $k$th and
$(k-1)$th iterations of Algorithm 1 due to the subspace approximation,
respectively, and $\gamma$ is as defined in (4.13).

###### Proof.

Let $\bm{X}_{k}$ be the exact solution of the Lyapunov equation (4.3) computed
at step $k$ of Algorithm 1. Using Propositions 4.4 and 4.5 we obtain

$\displaystyle\|\bm{X}-\tilde{\bm{X}}_{k}\|_{2}$
$\displaystyle\leq\|\bm{X}-\bm{X}_{k}\|_{2}+\|\bm{X}_{k}-\tilde{\bm{X}}_{k}\|_{2}\leq\gamma\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}+\frac{\|\bm{R}_{k}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})}.$
(4.21)

Furthermore, the triangle inequality and Propositions 4.4 and 4.7 imply

$\displaystyle\|\bm{X}_{k}-\bm{X}_{k-1}\|_{\mathrm{F}}$
$\displaystyle\leq\|\bm{X}_{k}-\tilde{\bm{X}}_{k}\|_{\mathrm{F}}+\|\tilde{\bm{X}}_{k}-\tilde{\bm{X}}_{k-1}\|_{\mathrm{F}}+\|\tilde{\bm{X}}_{k-1}-\bm{X}_{k-1}\|_{\mathrm{F}}$
$\displaystyle\leq\frac{\|\bm{R}_{k}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})}+\left(\|\tilde{\bm{Z}}_{k}\|_{\mathrm{F}}+\|\tilde{\bm{Z}}^{\mathsf{T}}_{k-1}\|_{F}\right)\|\tilde{\bm{Z}}_{k}-\tilde{\bm{Z}}_{k-1}\|_{\mathrm{F}}+\frac{\|\bm{R}_{k-1}\|_{\mathrm{F}}}{2\sigma_{\min}(\mathbf{\mathscr{A}})},$

yielding (4.20). ∎

###  4.3. Combining reduction error with offline numerical approximations

In this section, we show how to connect the MOR error for switched systems
with the numerical error introduced to efficiently approximate the system
Gramians. We will provide an error bound that, under some general numerical
assumptions, can be efficiently computed once the tolerance for the resolution
of the Gramians is fixed.

####  4.3.1. Relation between the exact Hankel singular values and their
approximation

The singular values of the Hankel matrix (3.3) are involved in estimating the
error between the output of the full and reduced systems; see (3.4). However,
one should observe that we never compute the exact solution $\bm{X}$ of the
GLE (4.4) and thus the exact Hankel matrix but only its approximation
$\bm{X}_{k}$ at the $k$th iteration of Algorithm 1 and then balance the
switched system using the approximated Gramians. If $\bm{X}_{k}\succeq\bm{X}$
for all $k$, then using the concept of generalized Gramians (see [13, Sec.
$4.7$]) introduced for standard systems without switching, one could replace
the singular values of $\bm{\mathcal{H}}$ in (3.4) with the ones of the
approximated Hankel matrix to keep a reliable error estimate. Unfortunately,
this is not generally the case as stated by the following result.

######  Proposition 4.9.

Suppose 4.1 and 4.2 are satisfied and let $\bm{X}$ and $\bm{X}_{k}$ denote the
unique solution of the GLE (4.4) and its approximation at the $k$th iteration
computed through Algorithm 1, respectively, where the Lyapunov equation (4.3)
is solved exactly. Then

$\bm{X}\succeq\bm{X}_{k}\quad\text{for all }k\in\mathbb{N}.$ (4.22)

###### Proof.

Let $\bm{E}_{k}=\bm{X}-\bm{X}_{k}$. It is sufficient to show that
$\bm{E}_{k}\succeq 0$ holds for all $k\in\mathbb{N}$. First, observe that
$\bm{E}_{0}\succeq 0$ since Algorithm 1 is (implicitly) initialized with
$\bm{X}_{0}=\mathbf{0}$ and $\bm{X}$ is symmetric positive definite. Next, let
us assume $\bm{E}_{k-1}\succeq 0$ for some $k\in\mathbb{N}$. By taking (4.4)
and subtracting (4.14), we obtain

$\mathbf{\mathscr{A}}\bm{E}_{k}+\bm{E}_{k}\mathbf{\mathscr{A}}^{\mathsf{T}}+\sum_{j=1}^{M}\left(\mathbf{\mathscr{F}}_{j}\bm{E}_{k-1}\mathbf{\mathscr{F}}_{k}^{\mathsf{T}}\right)=\mathbf{0},$

and since $\bm{E}_{k-1}\succeq 0$ by assumption it directly follows
$\bm{E}_{k}\succeq 0$ by the fact that $\mathbf{\mathscr{A}}$ is
asymptotically stable. The result follows via induction. ∎

Proposition 4.9 states that we cannot directly replace the singular values in
(3.4) with the ones arising by our balancing with the approximated Gramians.
Even when considering an approximation error $\tilde{\bm{E}}_{k}$ in the
solution of (4.3), we get $\hat{\bm{E}}_{k}={\bm{E}}_{k}+\tilde{\bm{E}}_{k}$.
By Proposition 4.9 we have ${\bm{E}}_{k}\succeq 0$ while $\tilde{\bm{E}}_{k}$
can be made arbitrarily smaller for each $k$, therefore it is reasonable to
assume that either $\hat{\bm{E}}_{k}\succeq 0$ or indefinite.

The following proposition relates the singular values of the exact and
approximated Hankel matrix in a computable way.

######  Proposition 4.10.

Consider matrices $\bm{\mathcal{H}}=\bm{S}^{\mathsf{T}}\bm{Z}$ and
$\tilde{\bm{\mathcal{H}}}=\tilde{\bm{S}}^{\mathsf{T}}\tilde{\bm{Z}}$ with
$\bm{S},\tilde{\bm{S}},\bm{Z},\tilde{\bm{Z}}\in\mathbb{R}^{k\times n}$. Then

$\sigma_{i}(\bm{\mathcal{H}})\leq\sigma_{i}(\tilde{\bm{\mathcal{H}}})+\sigma_{1}(\bm{\mathcal{H}}-\tilde{\bm{\mathcal{H}}})\qquad\text{for
all $i=1,\ldots,n$}$

and

$\sigma_{1}(\bm{\mathcal{H}}-\tilde{\bm{\mathcal{H}}})\leq\|\bm{S}-\tilde{\bm{S}}\|_{2}\|\bm{Z}-\tilde{\bm{Z}}\|_{2}+\|\tilde{\bm{S}}\|_{2}\|\bm{Z}-\tilde{\bm{Z}}\|_{2}+\|\bm{S}-\tilde{\bm{S}}\|_{2}\|\tilde{\bm{Z}}\|_{2}.$
(4.23)

###### Proof.

For two generic square matrices $\bm{M}_{1}$ and $\bm{M}_{2}$ of size $n\times
n$ we have

$\sigma_{i+j-1}(\bm{M}_{1}+\bm{M}_{2})\leq\sigma_{i}(\bm{M}_{1})+\sigma_{j}(\bm{M}_{2}),\qquad
1\leq i,j\leq n\quad\text{and}\quad i+j\leq n+1;$

see [19, Thm. 3.3.16]. The first result is thus a simple consequence setting
$j=1$. For (4.23) we have

$\displaystyle\sigma_{1}(\bm{\mathcal{H}}-\tilde{\bm{\mathcal{H}}})$
$\displaystyle=\|\bm{\mathcal{H}}-\tilde{\bm{\mathcal{H}}}\|_{2}=\|\bm{S}^{\mathsf{T}}\bm{Z}-\tilde{\bm{S}}^{\mathsf{T}}\tilde{\bm{Z}}\|_{2}=\|\bm{S}^{\mathsf{T}}\bm{Z}-\bm{S}^{\mathsf{T}}\tilde{\bm{Z}}+\bm{S}^{\mathsf{T}}\tilde{\bm{Z}}-\tilde{\bm{S}}^{\mathsf{T}}\tilde{\bm{Z}}\|_{2}$
$\displaystyle\leq\|\bm{S}^{\mathsf{T}}-\tilde{\bm{S}}^{\mathsf{T}}+\tilde{\bm{S}}^{\mathsf{T}}\|_{2}\|\bm{Z}-\tilde{\bm{Z}}\|_{2}+\|\bm{S}^{\mathsf{T}}-\tilde{\bm{S}}^{\mathsf{T}}\|_{2}\|\tilde{\bm{Z}}\|_{2}$
$\displaystyle\leq\|\bm{S}^{\mathsf{T}}-\tilde{\bm{S}}^{\mathsf{T}}\|_{2}\|\bm{Z}-\tilde{\bm{Z}}\|_{2}+\|\tilde{\bm{S}}^{\mathsf{T}}\|_{2}\|\bm{Z}-\tilde{\bm{Z}}\|_{2}+\|\bm{S}^{\mathsf{T}}-\tilde{\bm{S}}^{\mathsf{T}}\|_{2}\|\tilde{\bm{Z}}\|_{2}.\qed$

######  Remark 4.11.

The assumptions that the matrices $\bm{S}$ and $\tilde{\bm{S}}$ in Proposition
4.10 need to have the same dimension can be relaxed by adding zero rows until
the dimensions match.

To exploit Proposition 4.10, we need a computable upper bound for
$\|\bm{Z}-\tilde{\bm{Z}}\|_{2}$, which we derive in the following using
perturbation theory for eigenvectors and eigenvalues.

######  Theorem 4.12.

Let $\bm{X},\tilde{\bm{X}}\in\mathbb{R}^{n\times n}$ be symmetric and positive
definite and semi-definite, respectively, and let
$\bm{X}=\bm{U}_{\bm{X}}\bm{\Sigma}_{\bm{X}}\bm{U}_{\bm{X}}^{\top}$ and
$\tilde{\bm{X}}=\bm{U}_{\tilde{\bm{X}}}\bm{\Sigma}_{\tilde{\bm{X}}}\bm{U}_{\tilde{\bm{X}}}^{\mathsf{T}}$
denote their eigendecomposition with eigenvalues sorted in descending order.
Define $\bm{Z}\vcentcolon=\bm{U}_{\bm{X}}\bm{\Sigma}_{\bm{X}}^{1/2}$ and
$\tilde{\bm{Z}}\vcentcolon=\bm{U}_{\tilde{\bm{X}}}\bm{\Sigma}_{\tilde{\bm{X}}}^{1/2}$.
Then

$\displaystyle\|\bm{Z}-\tilde{\bm{Z}}\|_{2}$
$\displaystyle\leq\left(\|\tilde{\bm{Z}}^{\dagger}\|_{2}+C(\tilde{\bm{X}})\right)\|\bm{X}-\tilde{\bm{X}}\|_{2}+\sigma_{n_{\tilde{{\bm{X}}}}+1}(\bm{Z})+\mathcal{O}(\|\bm{X}-\tilde{\bm{X}}\|^{2}_{2}),$
(4.24)

where $\tilde{\bm{Z}}^{\dagger}$ denotes the pseudoinverse of
$\tilde{\bm{Z}}$,
$n_{\tilde{\bm{X}}}\vcentcolon=\operatorname{rank}(\tilde{\bm{X}})$,
$C(\tilde{\bm{X}})\vcentcolon=\sum_{i=1}^{n_{\tilde{{\bm{X}}}}}{\delta_{i}(\tilde{\bm{X}})}$,
and

$\delta_{i}(\tilde{\bm{X}})\vcentcolon=\lambda^{1/2}_{i}(\tilde{\bm{X}})\min_{\begin{subarray}{c}j=1,\ldots,n_{\tilde{\bm{X}}}+1\\\
\lambda_{i}\neq\lambda_{j}\end{subarray}}|\lambda_{i}(\tilde{\bm{X}})-\lambda_{j}(\tilde{\bm{X}})|,\quad
i=1,\ldots,n_{\tilde{\bm{X}}}.$

###### Proof.

Partition the matrices $\bm{U}_{\bm{X}}=[\bm{U}_{1},\bm{U}_{2}]$,
$\bm{\Sigma}_{\bm{X}}=\operatorname{diag}(\bm{\Sigma}_{1},\bm{\Sigma}_{2})$,
$\bm{U}_{\tilde{\bm{X}}}=[\tilde{\bm{U}}_{1},\tilde{\bm{U}}_{2}]$,
$\bm{\Sigma}_{\tilde{\bm{X}}}=\operatorname{diag}(\tilde{\bm{\Sigma}}_{1},0)$
with $\bm{U}_{1},\tilde{\bm{U}}_{1}\in\mathbb{R}^{n\times
n_{\tilde{\bm{X}}}}$,
$\bm{\Sigma}_{1},\tilde{\bm{\Sigma}}_{1}\in\mathbb{R}^{n_{\tilde{\bm{X}}}\times
n_{\tilde{\bm{X}}}}$. Then

$\displaystyle\begin{aligned}
\|\bm{Z}-\tilde{\bm{Z}}\|_{2}&\leq\|\bm{U}_{1}(\bm{\Sigma}_{1}^{1/2}-\tilde{\bm{\Sigma}}_{1}^{1/2})\|_{2}+\|(\bm{U}_{1}-\tilde{\bm{U}}_{1})\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{2}+\|\bm{U}_{2}\bm{\Sigma}_{2}^{1/2}\|_{2}\\\
&\leq\|\bm{\Sigma}_{1}^{1/2}-\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{2}+\|(\bm{U}_{1}-\tilde{\bm{U}}_{1})\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{2}+\sigma_{n_{\tilde{{\bm{X}}}}+1}(\bm{Z}).\end{aligned}$

For the first term, using Weyl’s theorem [18, Thm. 4.3.1], we have

$\displaystyle\begin{aligned}
\|\bm{\Sigma}_{1}^{1/2}-\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{2}&=\|(\bm{\Sigma}_{1}-\tilde{\bm{\Sigma}}_{1})(\bm{\Sigma}_{1}^{1/2}+\tilde{\bm{\Sigma}}_{1}^{1/2})^{-1}\|_{2}\leq\frac{\|\bm{X}-\tilde{\bm{X}}\|_{2}}{\lambda^{1/2}_{n_{\tilde{{\bm{X}}}}}(\bm{X})+\lambda^{1/2}_{n_{\tilde{{\bm{X}}}}}(\tilde{\bm{X}})}\\\
&\leq\frac{\|\bm{X}-\tilde{\bm{X}}\|_{2}}{\lambda^{1/2}_{n_{\tilde{{\bm{X}}}}}(\tilde{\bm{X}})}=\|\tilde{\bm{Z}}^{\dagger}\|_{2}\|\bm{X}-\tilde{\bm{X}}\|_{2}.\end{aligned}$

For the second term, we can write

$\displaystyle\|(\bm{U}_{1}-\tilde{\bm{U}}_{1})\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{2}$
$\displaystyle\leq\|(\bm{U}_{1}-\tilde{\bm{U}}_{1})\tilde{\bm{\Sigma}}_{1}^{1/2}\|_{\mathrm{F}}=\|\operatorname{vec}2((\bm{U}_{1}-\tilde{\bm{U}}_{1})\tilde{\bm{\Sigma}}_{1}^{1/2})\|_{2}$
$\displaystyle\leq\sum_{i=1}^{n_{\tilde{{\bm{X}}}}}\lambda^{1/2}_{i}(\tilde{\bm{X}})\|\bm{u}_{i}-\tilde{\bm{u}}_{i}\|_{2},$

where $\bm{u}_{i}$ and $\tilde{\bm{u}}_{i}$ is the $i$th column of
$\bm{U}_{1}$ and $\tilde{\bm{U}}_{1}$, respectively. and using eigenvectors
perturbation theory, see [52, Cha. V, Sec. 2.2.3], we get

$\|\bm{u}_{i}-\tilde{\bm{u}}_{i}\|_{2}\leq\max_{\begin{subarray}{c}j=1,\ldots,n_{\tilde{\bm{X}}}+1\\\
\lambda_{i}\neq\lambda_{j}\end{subarray}}\frac{\|\bm{X}-\tilde{\bm{X}}\|_{2}}{|\lambda_{i}(\tilde{{\bm{X}}})-\lambda_{j}(\tilde{\bm{X}})|}+\mathcal{O}(\|\bm{X}-\tilde{\bm{X}}\|^{2}_{2}),$
(4.25)

which completes the proof (noticing that (4.25) takes into account non-simple
eigenvalues). ∎

######  Remark 4.13.

The constant $C(\tilde{\bm{X}})$ can be explicitly computed once
$\tilde{\bm{Z}}$ is available. Moreover, the spectral gaps appearing in its
computation are not as relevant as they may appear. What matters is only the
last term
$\delta_{n_{\tilde{\bm{X}}}}(\tilde{\bm{X}})=\lambda^{-1/2}_{n_{\tilde{\bm{X}}}}(\tilde{{\bm{X}}})$.
Indeed, since non-simple eigenvalues are not problematic because it is
sufficient to consider the spectral gap concerning the set of eigenvalues
strictly different from the one considered, also small spectral gaps, say, for
instance, spectral gaps such that
$\delta_{i}(\tilde{\bm{X}})\leq\delta_{n_{\tilde{\bm{X}}}}(\tilde{\bm{X}})$,
could be ignored treating the two eigenvalues as a couple of non-simple
eigenvalues and looking to the spectral gaps with respects to a more distant
set of eigenvalues. This is possible because, in our MOR context, the basis
chosen to represent the subspace $\operatorname{span}(\tilde{\bm{X}})$ is not
relevant. What is relevant is that $\operatorname{span}(\tilde{\bm{X}})$ is a
good approximation
$\operatorname{span}(\bm{U}_{1}\bm{\Sigma}_{1}\bm{U}_{1}^{\mathsf{T}})$ and
the perturbation theory of invariant subspaces guarantees this, see [52, Cha.
V, Sec. 2].

####  4.3.2. A new certified error bound for the reduction error

We are now able to present our main result, which establishes a connection
between the balance truncation error bound for switched systems and the
approximated solutions achieved by using the projection subspace technique for
solving (4.3) within Algorithm 1.

######  Theorem 4.14.

Assume that the switched DAE (1.1) satisfies 2.7, the GLEs (3.1) with the
matrices defined in (3.8), satisfy 4.1, and let $\mathbf{\mathcal{P}}$ and
$\mathbf{\mathcal{Q}}$ denote the exact solutions. Assume that
$\mathbf{\mathcal{P}}$ and $\mathbf{\mathcal{Q}}$ satisfy the LMIs

$\displaystyle(\mathbf{\mathscr{F}}_{j}+\bm{J}_{1})\mathbf{\mathcal{P}}+\mathbf{\mathcal{P}}(\mathbf{\mathscr{F}}_{j}+\bm{J}_{1})^{\mathsf{T}}+\mathbf{\mathscr{B}}_{j}\mathbf{\mathscr{B}}_{j}^{\mathsf{T}}$
$\displaystyle\preceq\mathbf{0},$ (4.26a)
$\displaystyle(\mathbf{\mathscr{F}}_{j}+\bm{J}_{1})^{\mathsf{T}}\mathbf{\mathcal{Q}}+\mathbf{\mathcal{Q}}(\mathbf{\mathscr{F}}_{j}+\bm{J}_{1})+\mathbf{\mathscr{C}}_{j}^{\mathsf{T}}\mathbf{\mathscr{C}}_{j}$
$\displaystyle\preceq\mathbf{0},$ (4.26b)

for every $j\in\mathcal{J}$. Consider approximations
$\tilde{\mathbf{\mathcal{P}}}=\tilde{\bm{Z}}\tilde{\bm{Z}}^{\mathsf{T}}$ and
$\tilde{\mathbf{\mathcal{Q}}}=\tilde{\bm{S}}\tilde{\bm{S}}^{\mathsf{T}}$ with
$\tilde{\bm{Z}}\in\mathbb{R}^{n\times n_{\mathbf{\mathcal{P}}}}$ and
$\tilde{\bm{S}}\in\mathbb{R}^{n\times n_{\mathbf{\mathcal{Q}}}}$, the
associated ROM (1.3), and define
$\tilde{n}\vcentcolon=\min{(n_{\mathbf{\mathcal{P}}},n_{\mathbf{\mathcal{Q}}})}$.
Then, for every piecewise right continuous switching signal
$q\colon\mathbb{R}\to\mathcal{J}$ with locally finite number of jumps, the
relative output error between the FOM and the ROM is bounded by

$\displaystyle\begin{multlined}\frac{\|\bm{y}-\tilde{\bm{y}}\|_{L^{2}}}{\|\bm{u}\|_{L^{2}}}\leq
2(\tilde{n}-r)\left(c_{1}\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|_{2}+c_{2}\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|_{2}+\sigma_{n_{\mathbf{\mathcal{P}}}+1}({{\bm{Z}}})\|\tilde{\bm{S}}\|_{2}+\sigma_{n_{\mathbf{\mathcal{Q}}}+1}({{\bm{S}}})\|\tilde{\bm{Z}}\|_{2}\right)\\\
+2\\!\\!\sum_{i=r+1}^{\tilde{n}}\sigma_{i}(\tilde{\bm{\mathcal{H}}})+2\\!\\!\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}})+\mathcal{O}(\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|^{2}_{2})+\mathcal{O}(\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|^{2}_{2})+\mathcal{O}(\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|_{2}\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|_{2})\end{multlined}\frac{\|\bm{y}-\tilde{\bm{y}}\|_{L^{2}}}{\|\bm{u}\|_{L^{2}}}\leq
2(\tilde{n}-r)\left(c_{1}\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|_{2}+c_{2}\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|_{2}+\sigma_{n_{\mathbf{\mathcal{P}}}+1}({{\bm{Z}}})\|\tilde{\bm{S}}\|_{2}+\sigma_{n_{\mathbf{\mathcal{Q}}}+1}({{\bm{S}}})\|\tilde{\bm{Z}}\|_{2}\right)\\\
+2\\!\\!\sum_{i=r+1}^{\tilde{n}}\sigma_{i}(\tilde{\bm{\mathcal{H}}})+2\\!\\!\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}})+\mathcal{O}(\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|^{2}_{2})+\mathcal{O}(\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|^{2}_{2})+\mathcal{O}(\|\mathbf{\mathcal{P}}-\tilde{\mathbf{\mathcal{P}}}\|_{2}\|\mathbf{\mathcal{Q}}-\tilde{\mathbf{\mathcal{Q}}}\|_{2})$
(4.29)

with
$c_{1}\coloneqq\left(\|\tilde{\bm{Z}}^{\dagger}\|_{2}+C(\tilde{\mathbf{\mathcal{P}}})\right)\|\tilde{\bm{S}}\|_{2}$,
$c_{2}\coloneqq\left(\|\tilde{\bm{S}}^{\dagger}\|_{2}+C(\tilde{\mathbf{\mathcal{Q}}})\right)\|\tilde{\bm{Z}}\|_{2}$.

###### Proof.

2.7 and Theorem 2.6 allow us to study the reformulation (2.11), which itself
is equivalent to (3.7). By assumptions, the Gramians for this system satisfy
the LMIs (4.26) such that we can apply [38, Thm. 6] to obtain the error bound
(3.4). Using Proposition 4.10 to estimate the Hankel singular values, we
obtain

$\displaystyle\frac{\|\bm{y}-\tilde{\bm{y}}\|_{L^{2}}}{\|\bm{u}\|_{L^{2}}}$
$\displaystyle\leq
2\\!\\!\sum_{i=r+1}^{\tilde{n}}\sigma_{i}(\bm{\mathcal{H}})+2\\!\\!\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}})$
$\displaystyle\leq
2\\!\\!\sum_{i=r+1}^{\tilde{n}}\left(\sigma_{i}({\tilde{\bm{\mathcal{H}}}})+\sigma_{1}({\bm{\mathcal{H}}}-\tilde{\bm{\mathcal{H}}})\right)+2\\!\\!\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}})$
$\displaystyle=2(\tilde{n}-r)\sigma_{1}({\bm{\mathcal{H}}}-\tilde{\bm{\mathcal{H}}})+2\\!\\!\sum_{i=r+1}^{\tilde{n}}\sigma_{i}({\tilde{\bm{\mathcal{H}}}})+2\\!\\!\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}}).$

Then, using Proposition 4.10 in combination with Remark 4.11 and Theorem 4.12
we directly get the statement of Theorem 4.14. ∎

Several remarks are in order for the practical use of Theorem 4.8 and Theorem
4.14.

  1. (i)

The norms $\|\mathcal{P}-\tilde{\mathcal{P}}\|_{2}$ and
$\|\mathcal{Q}-\tilde{\mathcal{Q}}\|_{2}$ in (4.29) can be efficiently
estimated using Theorem 4.8, which in turn can be used to fix the exit
tolerance for the subspace projection employed in 5 and the exit tolerance for
Algorithm 1. In particular, considering $\gamma=1$, we can we require

$\displaystyle\|\bm{R}_{k}\|_{\mathrm{F}}\leq\tfrac{\sigma_{\min}(\mathbf{\mathscr{A}})}{3}\texttt{tol}\qquad\text{and}\qquad\left(\|\tilde{\bm{Z}}_{k}\|_{\mathrm{F}}+\|\tilde{\bm{Z}}_{k-1}\|_{\mathrm{F}}\right)\|\tilde{\bm{Z}}_{k}-\tilde{\bm{Z}}_{k-1}\|_{\mathrm{F}}$
$\displaystyle\leq\tfrac{1}{2}\texttt{tol}$

during Algorithm 1 to obtain $\|\bm{X}-\bm{X}_{k}\|_{2}\leq\texttt{tol}$.

  2. (ii)

The term $2\sum_{i=\tilde{n}+1}^{n}\sigma_{i}(\bm{\mathcal{H}})$ can be
reasonably assumed to be negligible with respect to the contribution of the
other terms in (4.29), thus we neglect it in practice.

  3. (iii)

The solution $\bm{X}$ of (4.4) is never formed explicitly, indeed only the
possible low-rank approximate Cholesky factor $\tilde{\bm{Z}}$ is determined.
In view of the proposed Krylov subspace method discussed in Section 4.2.1, the
Cholesky factor is obtained by computing the singular value decomposition of
the symmetric solution $\bm{Y}$ of (4.9), i.e.,

$\displaystyle\bm{Y}$
$\displaystyle=\bm{U}^{\bm{Y}}\bm{\Sigma}^{\bm{Y}}\left(\bm{U}^{\bm{Y}}\right)^{\mathsf{T}},$
$\displaystyle\bm{U}^{\bm{Y}}$
$\displaystyle\vcentcolon=\begin{bmatrix}\bm{U}^{\bm{Y}}_{1}&\bm{U}^{\bm{Y}}_{2}\end{bmatrix},$
$\displaystyle\bm{\Sigma}^{\bm{Y}}$
$\displaystyle\vcentcolon=\begin{bmatrix}\bm{\Sigma}^{\bm{Y}}_{1}&\mathbf{0}\\\
\mathbf{0}&\bm{\Sigma}^{\bm{Y}}_{2}\end{bmatrix},$

and then truncating using the approximation
$\tilde{\bm{Y}}\vcentcolon=\bm{U}^{\bm{Y}}_{1}\bm{\Sigma}^{\bm{Y}}_{1}(\bm{U}^{\bm{Y}}_{1})^{\mathsf{T}}$,
which implies
$\|\bm{Y}-\tilde{\bm{Y}}\|_{2}=\sigma_{1}(\bm{\Sigma}^{\bm{Y}}_{2})$. The
Cholesky factor is formed as
$\tilde{\bm{Z}}=\bm{V}\bm{U}^{\bm{Y}}(\bm{\Sigma}_{1}^{\bm{Y}})^{1/2}$, where
the columns of $\bm{V}$ form the basis of the Krylov subspace. In particular,
$\tilde{\bm{Z}}$ is in the form of Theorem 4.12. It is reasonable to choose
$\sigma_{1}(\bm{\Sigma}^{\bm{Y}}_{2})$ as the first singular value of $\bm{Y}$
smaller than tol. Then one can easily verify that
$\|\bm{Z}^{\dagger}\|_{2}\approx\tfrac{1}{\sqrt{\texttt{tol}}}$.

  4. (iv)

Using the previous observations, Remark 4.13, and the rapid decay of the
singular value of the Gramians, $C(\tilde{\mathbf{\mathcal{P}}})$
($C(\tilde{\mathbf{\mathcal{Q}}})$) $\approx\tfrac{1}{\sqrt{\texttt{tol}}}$.
Moreover, $\sigma_{n_{\mathbf{\mathcal{P}}}+1}(\tilde{{\bm{Z}}})$
($\sigma_{n_{\mathbf{\mathcal{Q}}}+1}(\tilde{{\bm{S}}})$) can be estimated as
${\sqrt{\texttt{tol}}}$.

  5. (v)

Normalizing the Gramians during Algorithm 1 gives
$\|\tilde{\bm{Z}}\|_{2},\|\tilde{\bm{S}}\|_{2}\approx 1$. Then, considering
the previous observations, the relative error bound (4.29) simplifies to

$\displaystyle\frac{\|\bm{y}-\tilde{\bm{y}}\|_{L^{2}}}{\|\bm{u}\|_{L^{2}}}$
$\displaystyle\lesssim
2(\tilde{n}-r)\left(6\sqrt{\texttt{tol}}\right)+2\\!\\!\sum_{i=r+1}^{\tilde{n}}\sigma_{i}(\tilde{\bm{\mathcal{H}}}),$
(4.30)

where we can reasonably neglect the second order terms. Relation (4.30) can be
used as a practical error estimate.

##  5\. Numerical Experiments

In this section we show numerical results to validate our method. We first
test our approach on the so called constrained mass-spring-damper system and
then we apply it to the instationary Stokes control system. With these test
problems we demonstrate that the MOR method developed, and the certified
numerics for the resolution of GLE, are effective.

All the computations were performed using MATLAB 2023a on a laptop with a
Apple M2 Pro processor. The codes for the test problems considered are
available in [30].

###  5.1. Constrained mass-spring-damper system

We consider the holonomically constrained mass-spring-damper system presented
in [31, Sec. 4]. The vibration of this system is described by the descriptor
system

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\bm{p}}(t)&=\bm{v}(t),&\bm{M}\dot{\bm{v}}(t)&=\bm{K}\bm{p}(t)+\bm{D}\bm{v}(t)-\bm{G}^{\mathsf{T}}\bm{\lambda}(t)+\bm{B}_{2}\bm{u}(t),\\\
\mathbf{0}&=\bm{G}\bm{p}(t),\\\
\bm{y}(t)&=\bm{C}_{1}\bm{p}(t),\end{aligned}\right.$ (5.1)

$0$$20$$40$$60$$-1$$0$$1$$t$$\bm{y}(t)$ and $\tilde{\bm{y}}(t)$FOM
($n=10001$ROM ($r=20$) (a) $\bm{u}(t)\coloneqq\sin(t)$

$0$$2$$4$$6$$8$$10$$12$$-0.4$$-0.2$$0$$0.2$$0.4$$t$$\bm{y}(t)$ and
$\tilde{\bm{y}}(t)$FOM ($n=10001$ROM ($r=20$) (b)
$\bm{u}(t)\coloneqq\sin(2\pi(\operatorname{e}^{\frac{t}{8}}))$

Figure 1. Constrained mass-spring-damper system with $M=5$ modes. Qualitative
comparison of full and reduced outputs (specifically second output) for two
different switching paths and input signals.

where $\bm{p}(t)\in\mathbb{R}^{g}$ is the position vector,
$\bm{v}(t)\in\mathbb{R}^{g}$ is the velocity vector,
$\bm{\lambda}(t)\in\mathbb{R}$ is the Lagrange multiplier,
$\bm{M}=\text{diag}(m_{1},\ldots,m_{g})$ is the mass matrix, $\bm{D}$ and
$\bm{K}$ are the tridiagonal damping and stiffness matrices,
$\bm{G}=[1,0,\ldots,0,-1]\in\mathbb{R}^{g}$ is the constraint matrix,
$\bm{B}_{2}=\bm{e}_{1}$, and
$\bm{C}_{1}=[\bm{e}_{1},\bm{e}_{2},\bm{e}_{g-1}]^{\mathsf{T}}$ where
$\bm{e}_{j}$ denotes the $j$th column of the identity matrix $\bm{I}_{g}$. The
descriptor system arising from (5.1) is of index $3$ and its associated
matrices read

$\displaystyle\bm{E}$
$\displaystyle\vcentcolon=\begin{bmatrix}\bm{I}_{g}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\bm{M}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}\end{bmatrix},$ $\displaystyle\bm{A}$
$\displaystyle\vcentcolon=\begin{bmatrix}\mathbf{0}&\bm{I}_{g}&\mathbf{0}\\\
\bm{K}&\bm{D}&-\bm{G}^{\mathsf{T}}\\\ \bm{G}&\mathbf{0}&\mathbf{0}\\\
\end{bmatrix},$ $\displaystyle\bm{B}$
$\displaystyle\vcentcolon=\begin{bmatrix}\mathbf{0}\\\ \bm{B}_{2}\\\
\mathbf{0}\\\ \end{bmatrix},$
$\displaystyle\bm{C}\vcentcolon=\begin{bmatrix}\bm{C}_{1}^{\mathsf{T}}\\\
\mathbf{0}\\\ \mathbf{0}\\\ \end{bmatrix}^{\mathsf{T}}.$ (5.2)

Considering $g$ masses, the state dimension of the system is $n=2g+1$, while
input and output dimensions are $m=1$ and $p=3$. The specific setting of
parameters is taken from [31, Sec. 4]. Note that the matrices in (5.2) are
sparse, and the kernel of $\bm{E}$ is of dimension one. To obtain a switched
system, we sample the matrices

$\displaystyle\bm{E}_{j}$
$\displaystyle\vcentcolon=\bm{E}+\mathcal{U}_{[0,1]}\bm{E},$
$\displaystyle\bm{A}_{j}$
$\displaystyle\vcentcolon=\bm{A}+\begin{bmatrix}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathcal{U}_{[0,0.35]}\bm{I}_{g}&\mathbf{0}\\\
0.5{\bm{e}_{j+1}^{\mathsf{T}}}&\mathbf{0}&\mathbf{0}\\\ \end{bmatrix},$
$\displaystyle\bm{B}_{j}$
$\displaystyle\vcentcolon=\bm{B}+\begin{bmatrix}\bm{e}_{j}\\\ \bm{e}_{j}\\\
\mathbf{0}\end{bmatrix}$ (5.3)

for $j=2,\ldots,M$, with $j=1$ given by (5.2), and where $\mathcal{U}_{[a,b]}$
denotes the uniform distribution in the interval $[a,b]$. With these choices,
all finite eigenvalues of $(\bm{E}_{j},\bm{A}_{j})$ have negative real part
and 2.7 (i) holds, i.e., there are no jumps in the state due to the input.
Moreover, the algebraic constraints are different for each mode and the input
matrices of each mode are linearly independent. For our numerical experiments
we choose $M=5$ modes and for clarity we only report the second output.

In our first experiment, we consider a FOM with $g=5000$ masses and construct
a ROM of size $r=20$ and compare the ROM and the FOM with two different input
functions and switching signals. The results are presented in Figure 1,
detailing that we can accurately represent the system dynamics for different
inputs and different switching signals with the reduced model. In our second
experiment, we construct ROMs for varying sizes of the reduced dimension. The
decay of the approximated Hankel singular values (red), the modified error
bound (4.30) derived from Theorem 4.14 (blue), and the relative reduction
error (green) with respect to the size of the reduced problem $r$ are given in
Figure 2(b). The GLEs to compute the approximated Gramians is solved with an
accuracy of $\texttt{tol}=10^{-12}$ such that (4.30) indicates that we can at
most predict an error around roughly $10^{-6}$. This is indeed confirmed by
the plot. The sum of the neglected singular values is not enough to guarantee
the upper bound of the reduction error when the reduced dimension is large
enough, while the error estimate derived in Theorem 4.14 is a reliable upper
bound. In our third experiment for the constrained mass-spring-damper system,
reported in Figure 2(b), we present the computational time to solve the GLEs
with respect to the size of the full problem $n$. The expected linear increase
confirms that our implementation for sparse matrices is effective.

$0$$10$$20$$30$$40$$10^{-11}$$10^{-7}$$10^{-3}$$10^{1}$$r$error bound
(4.29)$2\sum_{r+1}^{n}\sigma_{i}(\tilde{\mathbf{\mathcal{H}}})$$\|\bm{y}-\tilde{\bm{y}}\|_{L_{2}}/\|\bm{u}\|_{L_{2}}$
(a) Decay of the scaled reduction error and error estimate with respect to the
reduced dimension $r$.

$10^{2}$$10^{3}$$10^{4}$$10^{0}$$10^{1}$$10^{2}$$n$Computational time
(s)$\mathbf{\mathcal{P}}$$\mathbf{\mathcal{Q}}$$\mathcal{O}(n)$ (b)
Computational time to solve the GLEs with respect to the size of the full
problem $n$.

Figure 2. Constrained mass-spring-damper system. The input signal used to
produce Figure 2(a) is $\bm{u}(t)\coloneqq\sin(t)$.

Finally, to evaluate the impact of 2.7 (i), we consider a switched DAE with
two modes. The first mode $(\bm{E}_{1},\bm{A}_{1},\bm{B}_{1},\bm{C}_{1})$ is
equivalent to the reference system (5.2) and the second mode
$(\bm{E}_{2},\bm{A}_{2},\bm{B}_{2},\bm{C}_{2})$ is such that $\bm{E}_{2}$ is
equal to $\bm{E}_{1}$ but with the last non-zero row set to zero, $\bm{A}_{2}$
modifies $\bm{A}_{1}$ in such a way that the index of the DAE is still three
and that the finite eigenvalue of $(\bm{E}_{2},\bm{A}_{2})$ have negative real
part. The matrices $(\bm{B}_{2},\bm{C}_{2})$ are equal to
$(\bm{B}_{1},\bm{C}_{1})$. The resulting switched system does not satisfy 2.7
(i). Therefore, we compute the reachability Gramian by solving (3.10) instead
of (3.1a). Figure 3(a) compares the third output of the full and reduced
system. The ROM matches the output of the FOM well, and we can also observe
the jumps in the output variable, which, for this test example, are entirely
caused by the state jumps at the switching times. In Figure 3(b), we compare
the reduction error with respect to the size of the reduced problem for
reduced systems computed using $\mathbf{\mathcal{P}}$ from (3.10) (blue line)
and $\mathbf{\mathcal{P}}$ from the GLE using the matrices in (3.8) (red).
Moreover, to have a worst-case scenario, we artificially compute
$\mathbf{\mathcal{P}}$ from (3.8) but enforce orthogonality with respect to
the set $\mathcal{R}^{\mathrm{imp}}$ defined in (3.9) (green line). For the
constructed third case, we observe no decay in the reduction error since the
reduced system cannot approximate the jumps at the switching times. Moreover,
the blue line is always below the red line. In general, when input state jumps
are possible, one should always solve (3.10) to obtain $\mathbf{\mathcal{P}}$
since, without including the vectors originating the input jumps, there is no
guarantee that $\mathbf{\mathcal{P}}$ solving GLE without the terms
$\mathbf{\mathcal{B}}^{\mathrm{imp}}_{j}$ in (3.11) encloses the information
to well represents the input state jumps, therefore the reduction error may
not decay at all, as the green line shows.

$0$$2$$4$$6$$-20$$0$$20$$t$$\bm{y}(t)$ and $\tilde{\bm{y}}(t)$FOM $n=10001$ROM
$r=20$ (a) Qualitative comparison of the third output of the full and reduced
system

$0$$10$$20$$30$$10^{-10}$$10^{-6}$$10^{-2}$$10^{2}$$r$$\|\bm{y}-\tilde{\bm{y}}\|_{L_{2}}/\|\bm{u}\|_{L_{2}}$$\mathbf{\mathcal{P}}$
from (3.10)$\mathbf{\mathcal{P}}$ from GLE with
(3.8)$\mathbf{\mathcal{P}}\perp\mathcal{R}^{\mathrm{imp}}$ in (3.9) (b) Decay
of the scaled reduction error with respect to the reduced dimension $r$

Figure 3. Constrained mass-spring-damper system with $M=2$ allowing input
dependent state-jumps. The input used is $\bm{u}(t)\coloneqq\sin(t^{2}+t)$.

###  5.2. Instationary Stokes equations

The Stokes equations describe the flow of fluids at very low velocities
without convection and coincide with the linearization of the Navier-Stokes
equations around the zero-state. After a semi-discretization in space (see [9]
which is based on [44]), we obtain the descriptor system

$\displaystyle\left\\{\quad\begin{aligned}
\dot{\bm{v}}(t)&=\nu\bm{A}_{11}\bm{v}(t)+\bm{A}_{12}\bm{\rho}(t)+\bm{B}_{1}\bm{u}(t),\\\
\mathbf{0}&=\bm{A}^{\mathsf{T}}_{12}\bm{v}(t)+\bm{B}_{2}\bm{u}(t),\\\
\bm{y}(t)&=\bm{C}_{1}\bm{v}(t)+\bm{C}_{2}\bm{\rho}(t),\end{aligned}\right.$
(5.4)

where $\bm{v}(t)\in\mathbb{R}^{n_{\bm{v}}}$ and
$\bm{\rho}(t)\in\mathbb{R}^{n_{\bm{\rho}}}$ are the semidiscretized vectors of
velocities and pressures, respectively. The DAE (5.4) has index two and the
order $n=n_{\bm{v}}+n_{\bm{\rho}}$ of system (5.4) depends on the fineness of
the discretization and is usually very large, whereas the number $m$ of inputs
and the number $p$ of outputs are typically small. The representation of
system (5.4) in the standard DAE form reads as

$\displaystyle\bm{E}$ $\displaystyle=\begin{bmatrix}\bm{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix},$ $\displaystyle\bm{A}$
$\displaystyle=\begin{bmatrix}\nu\bm{A}_{11}&\bm{A}_{12}\\\
\bm{A}_{12}^{\mathsf{T}}&\mathbf{0}\end{bmatrix},$ $\displaystyle\bm{B}$
$\displaystyle=\begin{bmatrix}\bm{B}_{1}\\\ \bm{B}_{2}\end{bmatrix},$
$\displaystyle\bm{C}$
$\displaystyle=\begin{bmatrix}\bm{C}_{1}&\bm{C}_{2}\end{bmatrix}.$ (5.5)

To form a switched DAE we assume that system (5.4) is switching the diffusion
$\nu$ in the interval $[0.65,1.35]$. In particular, based on the number of
chosen modes, we select the values of $\nu$ uniformly spaced in this interval.
We consider $m=3$ and $p=3$ for the input and output dimensions, the matrices
$(\bm{E}_{j},\bm{B}_{j},\bm{C}_{j})$ are chosen as

$\displaystyle\bm{E}_{j}$
$\displaystyle\vcentcolon=\bm{E}+\mathcal{U}_{[0,1]}\bm{E},$
$\displaystyle\bm{B}_{j}$
$\displaystyle\vcentcolon=\bm{B}+\begin{bmatrix}\bm{e}_{j}&\bm{e}_{j+1}&\bm{e}_{n_{\bm{v}}+j}\end{bmatrix},$
$\displaystyle\bm{C}_{j}$
$\displaystyle\vcentcolon=\bm{C}+\begin{bmatrix}\bm{e}_{j}&\bm{e}_{j+1}&\bm{e}_{n_{\bm{v}}+j}\end{bmatrix}^{\mathsf{T}}$

for $j=2,\ldots,M$, with $j=1$ given by (5.5) and where $\bm{e}_{j}$ is the
$j$th column of the identity matrix $\bm{I}_{n}$. This choice allows verifying
2.7 (i) and to have a sparse orthonormal representation for the kernel of
$\bm{E}_{j}$ that has dimension $n_{\bm{\rho}}$, i.e., it increases linearly
with respect to the total dimension $n$ of the system.

$0$$1$$2$$3$$4$$5$$-0.2$$0$$0.2$$0.4$$t$$\bm{y}(t)$ and $\tilde{\bm{y}}(t)$FOM
($n=4719$)ROM ($r=18$) (a) Qualitative comparison of full and reduced outputs
(specifically third output)

$0$$20$$40$$60$$10^{-9}$$10^{-5}$$10^{-1}$$10^{3}$$r$error bound
(4.29)$2\sum_{r+1}^{n}\sigma_{i}(\tilde{\mathbf{\mathcal{H}}})$$\|\bm{y}-\tilde{\bm{y}}\|_{L_{2}}/\|\bm{u}\|_{L_{2}}$
(b) Decay of the scaled reduction error with respect to the reduced dimension
$r$

Figure 4. Instationary Stokes problem with $M=5$. The input used is
$\bm{u}(t)\coloneqq\sin(t^{2}+t)$

In Figure 4(a), we compare the FOM of dimension $n=4719$ with the ROM of
dimension $r=18$ for the input $\bm{u}(t)\coloneqq\sin(t^{2}+t)$, observing a
good match of the outputs. We emphasize that the output jumps are entirely due
to the output matrix $\bm{C}_{j}$ switching. In Figure 4(b), we demonstrate
that for the exact switching signal and same input as in Figure 4(a), the
classical error bound is indeed bounding the actual error without the
additional need for the modified error bound. The GLE for this test problem is
solved with an accuracy of $\texttt{tol}=10^{-12}$.

##  6\. Conclusions and perspectives

We introduced a novel MOR method for SDAEs. This is the first contribution of
its kind, specifically designed to handle large-scale systems with a generic
external switching signal. Importantly, our method provides an a priori bound
for the error between the output of the full and reduced problem, a feature
not previously available for MOR for switched descriptor systems. Our twofold
contribution significantly advances the field of model reduction for SDAEs:

  1. (i)

We show that a reduced order system can be constructed with a balancing
reduction method. First, this is shown for a particular class of switched
DAEs; then, we discuss how the method can handle also the generic class of
switched DAEs.

  2. (ii)

We show how to solve the GLE with prescribed accuracy efficiently. Then, we
use the obtained results to derive, assuming certain LMIs are verified, a
rigorous a priori error bound for the reduction error. We motivate
theoretically and practically why the classical balancing error bound cannot
be reliable for our approximation method.

The numerical results confirm our theoretical findings. The a priori error
bound relies on the fact that the exact solution of the stationary iteration
algorithm for the solution of the GLE satisfies the set of LMIs. Due to
rescaling operations, this assumption may not hold. Therefore, we highlight
that enforcing the verification of the LMIs for the computed Gramians is
crucial. A possible way to tackle the problem is to formulate it as a matrix
nearness problem where the computed Gramians are modified through a certain
minium-norm low-rank perturbation in order to satisfy the set of LMIs.

## Acknowledgment

MM acknowledges funding by the BMBF (grant no. 05M22VSA). BU is funded by
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC 2075 - 390740016. Both authors acknowledge
the support by the Stuttgart Center for Simulation Science (SimTech).

## References

  * [1] A. C. Antoulas.  Approximation of large-scale dynamical systems.  Adv. Des. Control. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2005.
  * [2] M. Baştuğ, M. Petreczky, R. Wisniewski, and J. Leth.  Model reduction by nice selections for linear switched systems.  IEEE Trans. Automat. Control, 61(11):3422–3437, 2016.
  * [3] P. Benner and T. Breiten.  Interpolation-based $\mathcal{H}_{2}$-model reduction of bilinear control systems.  SIAM J. Matrix Anal. Appl., 33(3):859–885, 2012.
  * [4] P. Benner and T. Breiten.  Low rank methods for a class of generalized Lyapunov equations and related issues.  Numer. Math., 124:441–470, 2013.
  * [5] P. Benner and T. Damm.  Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic systems.  SIAM J. Cont. Optim., 49(2):686–711, 2011.
  * [6] T. Berger, A. Ilchmann, and S. Trenn.  The quasi-Weierstraß form for regular matrix pencils.  Linear Algebra Appl., 436(10):4052–4069, 2012.
  * [7] T. Breiten and E. Ringh.  Residual-based iterations for the generalized Lyapunov equation.  BIT Numer. Math., 59:823–852, 2019.
  * [8] D. Cheng.  Stabilization of planar switched systems.  Systems Control Lett., 51:79–88, 2004.
  * [9] T. M. Community.  Stokes equation.  MORwiki – Model Order Reduction Wiki, 2018.
  * [10] M. Condon and R. Ivanov.  Nonlinear systems – algebraic gramians and model reduction.  Compel, 24(1):202–219, 2005.
  * [11] T. Damm.  Direct methods and ADI-preconditioned Krylov subspace methods for generalized Lyapunov equations.  Numer. Lin. Alg. Appl., 15(9):853–871, 2008.
  * [12] T. Damm.  Rational Matrix Equations in Stochastic Control.  Lecture Notes in Control and Information Sciences. Springer-Verlag, 2004\.
  * [13] G. Dullerud and F. Paganini.  A Course in Robust Control Theory a convex approach.  01 2000.
  * [14] F. R. Gantmacher.  The Theory of Matrices, volume 2.  Chelsea, New York, 1959.
  * [15] I. V. Gosea, M. Petreczky, and A. C. Anntoulas.  Data-driven model order reduction of linear switched systems in the Loewner framework.  SIAM J. Sci. Comput., 40(2):B572–B610, 2018.
  * [16] I. V. Gosea, M. Petreczky, A. C. Anntoulas, and C. Fiter.  Balanced truncation for linear switched systems.  Adv. Comput. Math., 44(6):1845–1886, 2018.
  * [17] I. V. Gosea, I. Pontes Duff, P. Benner, and A. C. Antoulas.  Model order reduction of switched linear systems with constrained switching.  In J. Fehr and B. Haasdonk, editors, IUTAM Symposium on Model Order Reduction of Coupled Systems, Stuttgart, Germany, May 22–25, 2018, volume 36 of IUTAM Bookseries, pages 41–53. Springer-Verlag, 2020\.
  * [18] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Matrix Analysis. Cambridge University Press, 1985.
  * [19] R. A. Horn and C. R. Johnson.  Topics in Matrix Analysis.  Cambridge University Press, 1991.
  * [20] M. S. Hossain.  Reduced realizations and model reduction for switched linear systems: a time-varying approach.  PhD thesis, University of Groningen, 2022.
  * [21] M. S. Hossain and S. Trenn.  Model reduction for switched differential-algebraic equations with known switching signal.  Technical report, University of Groningen, 2023.
  * [22] M. S. Hossain and S. Trenn.  Midpoint-based balanced truncation for switched linear systems with known switching signal.  IEEE Trans. Automat. Control, 69(1):535–542, 2024.
  * [23] E. Jarlebring, G. Mele, D. Palitta, and E. Ringh.  Krylov methods for low-rank commuting generalized Sylvester equations.  Numer. Lin. Alg. Appl., 28(6):e2176, 2018.
  * [24] D. Kressner and P. Sirković.  Truncated low-rank methods for solving general linear matrix equations.  Numer. Lin. Alg. Appl., 22(3):564–583, 2015.
  * [25] P. Kunkel and V. Mehrmann.  Differential-Algebraic Equations. Analysis and Numerical Solution.  European Mathematical Society, Zürich, 2006.
  * [26] F. Küsters and S. Trenn.  Duality of switched DAEs.  Math. Control Signals Systems, 28(3):25, 2016.
  * [27] P. Lancaster and M. Tismenetsky.  The Theory of Matrices: With Applications.  Computer Science and Scientific Computing. Elsevier Science, 2 edition, 1985.
  * [28] D. Liberzon.  Switching in Systems and Control.  Springer Science & Business Media, New York, NY, USA, 2003.
  * [29] D. Liberzon and S. Trenn.  Switched nonlinear differential algebraic equations: Solution theory, lyapunov functions, and stability.  Automatica, 48(5):954–963, 2012.
  * [30] M. Manucci.  Mattiamanucci/sdae-mor, April 2024\.
  * [31] V. Mehrmann and T. Stykel.  Balanced truncation model reduction for large-scale systems in descriptor form.  In P. Benner, D. C. Sorensen, and V. Mehrmann, editors, Dimension Reduction of Large-Scale Systems, volume 45 of Lecture Notes in Computational Science and Engineering, pages 83–115. Springer-Verlag, 2005\.
  * [32] V. Mehrmann and B. Unger.  Control of port-Hamiltonian differential-algebraic systems and applications.  Acta Numer., 32:395–515, 2023.
  * [33] N. Monshizadeh, H. L. Trentelman, and M. K. Camlibel.  A simultaneous balanced truncation approach to model reduction of switched linear systems.  IEEE Trans. Automat. Control, 47(12):3118–3131, 2012.
  * [34] A. V. Papadopoulos and M. Prandini.  Model reduction of switched affine systems: a method based on balanced truncation and randomized optimization.  In Proc. 17th Int. Conf. Hybrid Syst. Comput. Control, pages 113–122, Berlin, Germany, 2014.
  * [35] A. V. Papadopoulos and M. Prandini.  Model reduction of switched affine systems.  Automatica J. IFAC, 70:57–65, 2016.
  * [36] S. Peitz and S. Klus.  Koopman operator-based model reduction for switched-system control of PDEs.  Automatica J. IFAC, 106:184–101, 2019.
  * [37] M. Petreczky and I. V. Gosea.  Model reduction and realization theory of linear switched systems.  In C. Beattie, P. Benner, M. Embree, S. Gugercin, and S. Lefteriu, editors, Realization and Model Reduction of Dynamical Systems, pages 197–212. Springer-Verlag, 2022.
  * [38] M. Petreczky, R. Wisniewski, and J. Leth.  Balanced truncation for linear switched systems.  Nonlinear Anal. Hybrid Syst, 10:4–20, 2013.
  * [39] I. Pontes Duff, S. Grundel, and P. Benner.  New gramians for switched linear systems: Reachability, observability, and model reduction.  IEEE Trans. Automat. Control, 65(6):2526–2535, 2020.
  * [40] Y. Saad.  Numerical solution of large Lyapunov equations.  In Proceedings of the international symposium MTNS-89 signal processing, scattering, operator theory, and numerical methods., 1990.
  * [41] S. Sajja, M. Corless, E. Zeheb, and R. Shorten.  On dimensionality reduction and the stability of a class of switched descriptor systems.  Automatica J. IFAC, 49:1855–1860, 2013.
  * [42] S. Sajja, M. Corless, E. Zeheb, and R. Shorten.  Some stability tests for switched descriptor systems.  Automatica J. IFAC, 106:257–265, 2019.
  * [43] G. Scarciotti and A. Astolfi.  Model reduction for hybrid systems with state-dependent jumps.  IFAC-PapersOnLine, 49(18):850–855, 2016.
  * [44] M. Schmidt.  Systematic discretization of input/output maps and other contributions to the control of distributed parameter systems.  PhD thesis, TU Berlin, 2007.
  * [45] H. Schneider.  Positive operators and an inertia theorem.  Numer. Math., 7:11–17, 1965.
  * [46] P. Schulze and B. Unger.  Model reduction for linear systems with low-rank switching.  SIAM J. Cont. Optim., 56(6):4365–4384, 2018.
  * [47] H. Schwetlick and U. Schnabel.  Iterative computation of the smallest singular value and the corresponding singular vectors of a matrix.  Linear Algebra Appl., 371:1–30, 2003.
  * [48] H. R. Shaker and R. Wisniewski.  Model reduction of switched systems based on switching generalized Gramians.  Int. J. Innov. Comput. I., 8(7(B)):5025–5044, 2012.
  * [49] S. D. Shank, V. Simoncini, and D. B. Szyld.  Efficient low-rank solution of generalized Lyapunov equations.  Numer. Math., 134(2):327–342, 2016.
  * [50] V. Simoncini.  A new iterative method for solving large-scale lyapunov matrix equations.  SIAM Journal on Scientific Computing, 29(3):1268–1288, 2007.
  * [51] V. Simoncini.  Computational methods for linear matrix equations.  SIAM Rev., 58(3):377–441, 2016.
  * [52] G. W. Stewart and J. Sun.  Matrix perturbation theory.  Academic Press, Boston, 1990.
  * [53] Z. Sun and S. Ge.  Switched Linear Systems: Control and Design.  Springer-Verlag, 2005.
  * [54] S. Trenn.  Distributional differential algebraic equations.  Dissertation, Institut für Mathematik, Technische Universität Ilmenau, Universitätsverlag Ilmenau, Germany, 2009.
  * [55] S. Trenn.  Regularity of distributional differential algebraic equations.  Math. Control Signals Systems, 21(3):229–264, 2009.
  * [56] S. Trenn.  Switched differential algebraic equations.  In F. Vasca and L. Ianelli, editors, Dynamics and Control of Switched Electronic Systems. Advances in Industrial Control., pages 189–216. Springer-Verlag, 2012.
  * [57] S. Trenn.  Solution concepts for linear DAEs: a survey.  In A. Ilchmann and T. Reis, editors, Surveys in Differential-Algebraic Equations I, chapter 3, pages 137–172. Springer, Berlin/Heidelberg, 1st edition, 2013.
  * [58] B. Unger and S. Gugercin.  Kolmogorov $n$-widths for linear dynamical systems.  Adv. Comput. Math., 45(5–6):2273–2286, 2019.
  * [59] K.-T. Wong.  The eigenvalue problem $\lambda Tx+Sx$.  J. Differential Equations, 16(2):270–280, 1974.
  * [60] W. M. Wonham.  On a matrix Riccati equation of stochastic control.  SIAM J. Cont., 6(4):681–697, 1968.
  * [61] L. Wu and W. X. Zheng.  Weighted $\mathcal{H}_{\infty}$ model reduction for linear switched systems with time-varying delay.  Automatica J. IFAC, 45:186–193, 2009.