Stability-preserving model order reduction for linear stochastic Galerkin systems

TL;DR

This paper investigates methods to preserve stability in model order reduction of linear stochastic Galerkin systems, ensuring reliable reduced models for high-dimensional uncertainty quantification in engineering applications.

Contribution

It compares two stability-preserving transformations for Galerkin systems, analyzing their properties and demonstrating their effectiveness through numerical examples.

Findings

Both transformation techniques effectively preserve stability.

Numerical results show accurate reduced models for mechanical and electrical systems.

Different properties of the two methods influence their numerical performance.

Abstract

Mathematical modeling often yields linear dynamical systems in science and engineering. We change physical parameters of the system into random variables to perform an uncertainty quantification. The stochastic Galerkin method yields a larger linear dynamical system, whose solution represents an approximation of random processes. A model order reduction (MOR) of the Galerkin system is advantageous due to the high dimensionality. However, asymptotic stability may be lost in some MOR techniques. In Galerkin-type MOR methods, the stability can be guaranteed by a transformation to a dissipative form. Either the original dynamical system or the stochastic Galerkin system can be transformed. We investigate the two variants of this stability-preserving approach. Both techniques are feasible, while featuring different properties in numerical methods. Results of numerical computations are demonstrated for two test examples modeling a mechanical application and an electric circuit, respectively.

Tables

Table 1. Table 1: Properties of stochastic Galerkin system in mass-spring-damper example.

dimension	11400
number of outputs	1140
# non-zeros in $\hat{E}$	13110
# non-zeros in $\hat{A}$	50958
spectral abscissa of $(\hat{E},\hat{A})$	$-0.0036$
spectral abscissa of $\hat{A}+{\hat{A}}^{\top }$	42.69

Table 2. Table 2: Number of asymptotically stable systems out of 100 ROMs with dimensions 1 , … , 100 1 … 100 1,\ldots,100 and computation times for projection matrix ( 10 ) in stabilization using frequency domain integrals for mass-spring-damper example.

# nodes in quadrature		10		20		30		40		untransf.
# stable ROMs		79		91		96		100		38
computation time (s)		19.0		38.0		56.8		75.6		–

Table 3. Table 3: Number of non-zero entries (nnz) in matrices of stochastic Galerkin systems for mass-spring-damper example.

	nnz in $\hat{E}$		nnz in $\hat{A}$
untransformed Galerkin system	$\mathrm{13\hspace{0.17em}110}$	(0.01%)	$\mathrm{50\hspace{0.17em}958}$	(0.04%)
Galerkin system in (ii)	$\mathrm{30\hspace{0.17em}000\hspace{0.17em}000}$	(23.1%)	$\mathrm{29\hspace{0.17em}999\hspace{0.17em}562}$	(23.1%)

Table 4. Table 4: Computation times for MOR method producing reduced dimension r = 100 𝑟 100 r=100 and additional stabilization techniques in mass-spring-damper example.

	computation time (s)
Arnoldi method for original system	2.4
Arnoldi method for system in (ii)	65.4
technique (i)	75.6
technique (ii)	2591.9
technique (iii)	0.1

Table 5. Table 5: Properties of regularized stochastic Galerkin system in band-pass filter example.

dimension	6900
number of outputs	300
# non-zeros in ${\hat{E}}_{\mathrm{reg}}$	21912
# non-zeros in ${\hat{A}}_{\mathrm{reg}}$	21912
spectral abscissa of $({\hat{E}}_{\mathrm{reg}},{\hat{A}}_{\mathrm{reg}})$	$-4.94\cdot {10}^{-4}$
spectral abscissa of $\left({\hat{E}}_{\mathrm{reg}}^{\top }{\hat{A}}_{\mathrm{reg}}\right)+{\left({\hat{E}}_{\mathrm{reg}}^{\top }{\hat{A}}_{\mathrm{reg}}\right)}^{\top }$	$1.28\cdot {10}^{-5}$

Table 6. Table 6: Number of asymptotically stable systems out of 100 ROMs with dimensions 1 , … , 100 1 … 100 1,\ldots,100 in MOR of stochastic Galerkin system for band-pass filter example.

	FOM system	# stable ROMs
i)	DAE	22
ii)	regularized system (ODE)	59
iii)	transformed as in Section 4.1 (ODE)	100
iv)	transformed as in Section 4.3 (ODE)	100

Equations

$$\begin{array}{rl}E\dot{x}(t)& =Ax(t)+Bu(t)\\ y(t)& =Cx(t)\end{array}$$

$$H(s)=C(sE-A{)}^{-1}B,$$

$$\Vert H{\Vert }_{{\mathcal{H}}_2}=\sqrt{\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\Vert H(\mathrm{i}\omega ){\Vert }_{\mathrm{F}}^2\mathrm{d}\omega }$$

$$\begin{array}{rl}\bar{E}\dot{\bar{x}}(t)& =\bar{A}\bar{x}(t)+\bar{B}u(t)\\ \bar{y}(t)& =\bar{C}\bar{x}(t)\end{array}$$

$$\bar{A}=W^{\top }AV,\bar{B}=W^{\top }B,\bar{C}=CV,\bar{E}=W^{\top }EV,$$

$$\underset{t\ge 0}{\sup }\Vert y(t)-\bar{y}(t){\Vert }_{\infty }\le {\Vert H-\bar{H}\Vert }_{{\mathcal{H}}_2}\Vert u{\Vert }_{{\mathcal{L}}^2[0,\infty )}$$

$$\dot{z}(t)=L^{-1}A{L}^{-\top }z(t)+L^{-1}Bu(t)$$

$$S:=\tilde{A}+{\tilde{A}}^{\top }=L^{-1}A{L}^{-\top }+{\left(L^{-1}A{L}^{-\top }\right)}^{\top }=L^{-1}\left(A+A^{\top }\right)L^{-\top }.$$

$$2\mathrm{Re}(\lambda )=\lambda +\overline{\lambda }=\left(\lambda +\overline{\lambda }\right)v^{\ast }v=v^{\ast }\left(\tilde{A}+{\tilde{A}}^{\top }\right)v=v^{\ast }Sv<0.$$

$$A^{\top }ME+E^{\top }MA<0,$$

$$\begin{array}{rl}{E}^{\top }ME\dot{x}(t)& =E^{\top }MAx(t)+E^{\top }MBu(t)\\ y(t)& =Cx(t).\end{array}$$

$$W=MEV.$$

$$A^{\top }ME+E^{\top }MA+F=0$$

$$\Vert \tilde{M}-M\Vert <\frac{1}{\Vert A^{\top }\Vert \cdot \Vert E\Vert +\Vert A\Vert \cdot \Vert E^{\top }\Vert }$$

$$\begin{array}{rl}E(\mu )\dot{x}(t,\mu )& =A(\mu )x(t,\mu )+B(\mu )u(t)\\ y(t,\mu )& =C(\mu )x(t,\mu ).\end{array}$$

$$x(0,\mu )=x_0(\mu )$$

$$\mathbb{E}[f]={\int }_{\Omega }f(\mu (\omega ))\mathrm{d}P(\omega )={\int }_{\mathcal{M}}f(\mu )\rho (\mu )\mathrm{d}\mu$$

$$⟨f,g⟩={\int }_{\mathcal{M}}f(\mu )g(\mu )\rho (\mu )\mathrm{d}\mu$$

$${\mathcal{L}}^2(\mathcal{M},\rho )=\left\{f:\mathcal{M}\to \text{mathbbm}R:f\mathrm{measurable},\mathbb{E}[f^2]<\infty \right\}.$$

$${\Phi }_i(\mu )={\Psi }_{i_1}^{(1)}({\mu }_1)\cdot {\Psi }_{i_2}^{(2)}({\mu }_2)\cdot \dots \cdot {\Psi }_{i_{q-1}}^{(q-1)}({\mu }_{q-1})\cdot {\Psi }_{i_q}^{(q)}({\mu }_q),$$

$$y(t,\mu )=\sum _{i=1}^{\infty }{w}_i(t){\Phi }_i(\mu )$$

$$x(t,\mu )=\sum _{i=1}^{\infty }{v}_i(t){\Phi }_i(\mu )$$

$$\mathbb{E}[y(t,\cdot )]=w_1(t)\text{mbox}and\mathrm{Var}[y(t,\cdot )]=\sum _{i=2}^{\infty }{w}_i(t{)}^2$$

$$m=\frac{(d+q)!}{d!q!}.$$

$$y^{(m)}(t,\mu )=\sum _{i=1}^m{w}_i(t){\Phi }_i(\mu )\text{mbox}and{x}^{(m)}(t,\mu )=\sum _{i=1}^m{v}_i(t){\Phi }_i(\mu ).$$

$$\begin{array}{rl}\hat{E}\dot{\hat{v}}(t)& =\hat{A}\hat{v}(t)+\hat{B}u(t)\\ \hat{w}(t)& =\hat{C}\hat{v}(t)\end{array}$$

$$S(\mu )=({\Phi }_i(\mu ){\Phi }_j(\mu ){)}_{i,j=1,\dots ,m}\text{mbox}ands(\mu )=({\Phi }_i(\mu ){)}_{i=1,\dots ,m}.$$

$$\hat{A}=\mathbb{E}[S\otimes A],\hat{B}=\mathbb{E}[s\otimes B],\hat{C}=\mathbb{E}[S\otimes C],\hat{E}=\mathbb{E}[S\otimes E]$$

$${\hat{y}}^{(m)}(t,\mu )=\sum _{i=1}^m{\hat{w}}_i(t){\Phi }_i(\mu ),$$

$$\begin{array}{rl}{z}^{\top }\hat{E}z& =\sum _{i,j=1}^m{z}_i^{\top }{\hat{E}}_{ij}{z}_j=\sum _{i,j=1}^m{z}_i^{\top }\mathbb{E}[{\Phi }_i{\Phi }_{j}E]z_j=\mathbb{E}\left[\sum _{i,j=1}^m{z}_i^{\top }E{z}_j{\Phi }_i{\Phi }_j\right]\\ & =\mathbb{E}\left[{\left(\sum _{i=1}^m{z}_i{\Phi }_i\right)}^{\top }E\left(\sum _{j=1}^m{z}_j{\Phi }_j\right)\right]\ge 0,\end{array}$$

Full text

**Stability-preserving model order reduction

for linear stochastic Galerkin systems**

Roland Pulch

Institute of Mathematics and Computer Science, University of Greifswald,

Walther-Rathenau-Str. 47, D-17489 Greifswald, Germany.

Email: roland.pulch@uni-greifswald.de

Abstract

[TABLE]

1 Introduction

Numerical simulation of mathematical models represents the main issue in scientific computing. We consider linear dynamical systems, which play an important role in mechanics and electrical engineering, for example. Furthermore, uncertainty quantification becomes more and more relevant in many fields of applications, see [17], for instance. A common approach is to replace uncertain parameters by random variables, see [37, 39]. Statistics of the stochastic model can be computed by sampling methods or quadrature rules. Alternatively, the stochastic Galerkin method changes the random-dependent linear dynamical system into a larger deterministic linear dynamical system.

The dimension of the stochastic Galerkin system becomes huge in the case of large numbers of random variables. Methods of model order reduction (MOR) are able to decrease the complexity. Transient solutions of a reduced system allow for an efficient numerical simulation. Several MOR methods are available for general linear dynamical systems, see [1, 2, 4, 33]. MOR of linear stochastic Galerkin systems was also examined in several previous works [9, 19, 32, 26, 27, 40]. MOR of nonlinear stochastic Galerkin systems was considered in [30].

However, even though the linear stochastic Galerkin system is asymptotically stable, the reduced Galerkin system often looses this stability in some MOR techniques. We investigate stability-preserving strategies in the case of Galerkin-type projection-based MOR like the Arnoldi method or proper orthogonal decomposition, for example. Galerkin-type MOR can be applied to any linear dynamical system (not only stochastic Galerkin systems). A dissipativity property guarantees the preservation of stability. If a general linear dynamical system does not satisfy the dissipativity property, then it can be transformed into a dissipative structure, see [7, 24, 29]. The crucial part to identify a transformation consists in the solution of a Lyapunov equation. Direct methods, see [11], or approximate methods, see [22, 23, 35, 38], yield the numerical solutions of Lyapunov equations.

We examine the stability-preserving approach in the case of linear stochastic Galerkin systems consisting of ordinary differential equations. Several variants are feasible. The high-dimensional Galerkin-projected system is transformed or, vice versa, the original systems are transformed followed by a Galerkin projection. We analyze the two strategies and another variant.

In addition, network approaches produce models consisting of differential-algebraic equations in industrial applications. Thus we extend the stability-preserving techniques to this class of problems. The Lyapunov equations have no solution now. Therefore, we use a regularization technique, which was also employed in [20].

We apply the analyzed techniques to mathematical models of two test examples: a mass-spring-damper system and an electric circuit of a band-pass filter.

2 Stability preservation in reduction

We review a concept for stability preservation in Galerkin-type projection-based MOR for general linear dynamical systems.

2.1 Linear dynamical systems

We consider linear dynamical systems in the form

[TABLE]

with constant matrices $A,E\in\mathbbm{R}^{n\times n}$, $B\in\mathbbm{R}^{n\times n_{\rm in}}$, and $C\in\mathbbm{R}^{n_{\rm out}\times n}$. The state variables or inner variables are $x:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{n}$. Inputs $u:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{n_{\rm in}}$ are supplied to the system. The outputs $y:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{n_{\rm out}}$ are defined as quantities of interest (QoI). Initial value problems are given by $x(0)=x_{0}$.

If the mass matrix $E$ is non-singular, then the system (1) consists of ordinary differential equations (ODEs). If the mass matrix $E$ is singular, then the system (1) represents differential-algebraic equations (DAEs). Furthermore, we assume that the system satisfies the following stability condition.

We specify some common notions.

Definition 1**.**

Given a matrix pencil $(E,A)$ with $E,A\in\mathbbm{R}^{n\times n}$, the set of eigenvalues is $\Lambda=\left\{\lambda\in\mathbbm{C}\,:\,\det(\lambda E-A)=0\right\}$ and the spectral abscissa reads as $\alpha(E,A)=\max\left\{{\rm Re}(\lambda)\,:\,\lambda\in\Lambda\right\}$. The spectral abscissa of a single matrix $A$ is the spectral abscissa of the matrix pencil $(I,A)$ with the identity matrix $I\in\mathbbm{R}^{n\times n}$.

Definition 2**.**

A matrix pencil $(E,A)$ with $E,A\in\mathbbm{R}^{n\times n}$ is called regular, if there is (at least one) $\lambda\in\mathbbm{C}$ such that $\det(\lambda E-A)\neq 0$.

Definition 3**.**

A linear dynamical system (1) is called asymptotically stable, if the involved matrix pencil $(E,A)$ has a spectral abscissa satisfying $\alpha(E,A)<0$.

The asymptotic stability implies that the generalized eigenvalue problem has only a finite number of eigenvalues, which all exhibit a negative real part. We assume that the system (1) satisfies this stability condition. In the case of ODEs, the matrix pencil is always regular, even if the system is unstable. In the case of DAEs, the asymptotic stability implies that the matrix pencil $(E,A)$ is regular.

2.2 Transfer functions and Hardy norms

The input-output behavior of a linear dynamical system (1) can be described in the frequency domain, see [1] for ODEs and [6] for DAEs.

Definition 4**.**

The transfer function of a linear dynamical system (1) is the mapping $H:\mathbbm{C}\backslash\Lambda\rightarrow\mathbbm{C}^{n_{\rm out}\times n_{\rm in}}$ with

[TABLE]

where $\Lambda$ represents the set of eigenvalues from Definition 1.

The transfer function is always a rational function, whose poles are the eigenvalues. The regularity of the matrix pencil guarantees a finite set of poles. The magnitude of a transfer function can be measured by Hardy norms, see [1].

Definition 5**.**

The $\mathcal{H}_{2}$-norm of a transfer function reads as

[TABLE]

including the Frobenius matrix norm $\|\cdot\|_{\rm F}$, the angular frequency $\omega$, and ${\rm i}=\sqrt{-1}$.

Since we assume asymptotically stable linear dynamical systems (1), the transfer function is defined on the complete imaginary axis. In the case of ODEs, the $\mathcal{H}_{2}$-norm is always finite. In the case of DAEs, the existence of the $\mathcal{H}_{2}$-norm is not guaranteed. Nevertheless, the $\mathcal{H}_{2}$-norm is quite often finite for DAEs of index 1 or 2. The norm (3) may also exist for unstable systems, if there is no pole on the imaginary axis.

2.3 Projection-based model order reduction

Projection matrices $V,W\in\mathbbm{R}^{n\times r}$ of full rank are specified with $r\ll n$. Concerning the full-order model (FOM) in (1), the reduced-order model (ROM) reads as

[TABLE]

with state variables or inner variables $\bar{x}:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{r}$. Initial values $\bar{x}(0)=\bar{x}_{0}$ are supposed. We obtain the matrices via

[TABLE]

which is also called a Petrov-Galerkin-type MOR. A Galerkin-type projection-based MOR is characterized by $W=V$, where just one projection matrix has to be determined. Important examples are the one-sided Arnoldi method and the proper orthogonal decomposition (POD), see [1].

Each linear dynamical system is described by a transfer function (2) in the frequency domain. The difference between the transfer functions $H$ of FOM and $\bar{H}$ of ROM quantifies the error of the MOR. Hardy norms like the $\mathcal{H}_{2}$-norm from Definition 5, for example, can be applied to their transfer functions. However, a small error in the frequency domain implies a small error in the time domain only if both systems are asymptotically stable. It holds that

[TABLE]

with the maximum vector norm $\|\cdot\|_{\infty}$ and the $\mathcal{L}^{2}[0,\infty)$-norm in time provided that all initial values are zero, see [3]. The $\mathcal{H}_{2}$-norm of Definition 5 is a strong measure. Hence there are neither a priori error bounds nor cheap a posteriori error bounds available in MOR techniques. Just an approximation of the $\mathcal{H}_{2}$-norm in (6) can be computed posterior, where the computational effort is dominated by evaluations of the transfer function in the FOM. The balanced truncation method yields an a priori error bound in the $\mathcal{H}_{\infty}$-norm, see [1, p. 212].

2.4 Dissipative systems

In balanced truncation, see [1], the ROM (4) is always asymptotically stable provided that the system (1) is asymptotically stable. Yet the asymptotic stability may be lost in the ROM (4) within other MOR methods like Krylov subspace techniques, see [10], and POD, for example.

Using Galerkin-type MOR, the stability is guaranteed for some classes of linear dynamical systems. In the case of ODEs, we define the following type of system.

Definition 6**.**

A linear dynamical system (1) is called dissipative, if

$E$* is symmetric as well as positive definite, and* 2. 2.

$A+A^{\top}$* is negative definite.*

The above condition represents a dissipativity of the matrix $A$ as shown in [21]. Other definitions of dissipative systems are used in the literature. We prove a property, which was shown for an explicit system of ODEs in [24].

Theorem 1**.**

If the linear dynamical system (1) is dissipative with respect to Definition 6, then it is also asymptotically stable as in Definition 3.

Proof:

Let $E=LL^{\top}$ be the Cholesky decomposition of the mass matrix. The system (1) is equivalent to the explicit system of ODEs

[TABLE]

with $z=L^{\top}x$ and $y=CL^{-\top}z$. The symmetric part of the involved system matrix $\tilde{A}$ reads as

[TABLE]

It holds that $v^{*}Sv=(L^{-\top}v)^{*}(A+A^{\top})(L^{-\top}v)$ for any $v\in\mathbbm{C}^{n}\backslash\{0\}$. Thus $S$ is negative definite. Let $\lambda\in\mathbbm{C}$ be an eigenvalue of $\tilde{A}$ and $v$ be an associated eigenvector satisfying $v^{*}v=1$. It follows that

[TABLE]

Thus the systems (7) and (1) are asymptotically stable. $\Box$

In contrast, the asymptotic stability does not imply the dissipativity of Definition 6, even if the mass matrix $E$ is symmetric and positive definite. An example will be presented in Section 6.1. Now we consider Galerkin-type MOR.

Theorem 2**.**

If the linear dynamical system (1) is dissipative, then a Galerkin-type MOR yields a dissipative reduced system (4). Hence the reduced system is asymptotically stable.

The proof can be found in [29], for example.

2.5 Transformations

The asymptotic stability of a linear dynamical system is invariant with respect to basis transformations. In contrast, if a system of ODEs is not dissipative, then it can be converted to an equivalent dissipative form by a basis transformation in the state space, see [24]. Alternatively, a basis transformation is feasible in the image space only, see [7]. We require a symmetric positive definite solution $M\in\mathbbm{R}^{n\times n}$ of the Lyapunov inequality

[TABLE]

which means that the matrix on the left-hand side of (8) is negative definite. If the mass matrix $E$ is non-singular and the matrix pencil $(E,A)$ satisfies the Definition 3 of asymptotic stability, then an infinite set of solutions $M$ exists.

We change the system of ODEs (1) into the equivalent system

[TABLE]

The transformed system exhibits the desired dissipativity property.

Theorem 3**.**

If the asymptotically stable linear dynamical system (1) has a non-singular mass matrix $E$, then the transformed system (9) with $M$ satisfying (8) is dissipative in view of Definition 6.

We use a Galerkin-type MOR with a projection matrix $V$ to the transformed system (9). This approach can be written as a Petrov-Galerkin-type MOR applied to the original system (1) with matrices (5) and the projection matrix

[TABLE]

Thus we do not require to calculate the transformed system (9) explicitly. Alternatively, we compute the projection matrix (10). However, the original Galerkin-type MOR of the system (1) is not equivalent to the MOR of the system (9).

2.6 Numerical solution of Lyapunov inequality

We solve the Lyapunov inequality (8) using a Lyapunov equation

[TABLE]

including a predetermined symmetric positive definite matrix $F\in\mathbbm{R}^{n\times n}$. This matrix represents a degree of freedom, because any choice yields a symmetric positive definite solution of the Lyapunov inequality. A simple admissible choice is the identity matrix $I_{n}\in\mathbbm{R}^{n\times n}$. Moreover, we do not need to solve the Lyapunov equation (11) with a high accuracy, because a rough approximation $\widetilde{M}$ often still satisfies the Lyapunov inequality (8). In [28], it was shown that any approximation $\widetilde{M}$ of the exact solution $M$ of the Lyapunov equation (11) for $F=I_{n}$ with the property

[TABLE]

in some subordinate matrix norm also solves the Lyapunov inequality (8). The condition (12) is just sufficient and not necessary. However, approximations $\widetilde{M}$ satisfying (12) may have a relative error of up to 100%, see [28], which motivates that rough estimates can solve the problem.

There are direct methods to compute a solution $M$ of (11) or a symmetric decomposition $M=LL^{\top}$, see [11, 22]. Their computational effort is typically $\mathcal{O}(n^{3})$. In the high-dimensional case, we have to use approximate methods to decrease the computation work. The following techniques are available:

• i)

projection methods (Krylov subspace techniques, POD, etc.), see [13, 38],

• ii)

alternating direction implicit (ADI) iteration, see [15, 23],

• iii)

frequency domain integrals, see [5, 28],

and others. In the cases (i) and (ii), the methods yield an approximation $\widetilde{M}=ZZ^{\top}$ with a low-rank factor $Z\in\mathbbm{R}^{n\times k}$ ($k\ll n$). Thus the transformation is given by a singular matrix $\widetilde{M}$. It follows that the mass matrix $\bar{E}$ of the reduced system (4) may become singular or ill-conditioned, as shown in [29]. In contrast, the method (iii) from [28] computes the projection matrix (10), where the underlying approximation $\widetilde{M}$ is always non-singular. However, the matrix $\widetilde{M}$ is never computed but a matrix-matrix product with this approximation. In the frequency domain integral approach, the projection matrix $V$ has to be determined by the original linear dynamical system (1).

3 Stochastic Galerkin systems

We illustrate the concept of the stochastic Galerkin method and define our problem under investigation.

3.1 Random linear dynamical systems

We include parameters in the linear dynamical systems and obtain

[TABLE]

The matrices $A,E\in\mathbbm{R}^{n\times n}$, $B\in\mathbbm{R}^{n\times n_{\rm in}}$, and $C\in\mathbbm{R}^{n_{\rm out}\times n}$ depend on parameters $\mu\in\mathcal{M}\subseteq\mathbbm{R}^{q}$. We assume that the dimension $n$ is independent of the number of parameters $q$. The values $\mu$ may represent physical parameters or artificial parameters.

Thus the state variables or inner variables $x:[0,t_{\rm end}]\times\mathcal{M}\rightarrow\mathbbm{R}^{n}$ depend on time as well as the parameters. The inputs $u:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{n_{\rm in}}$ are independent of the parameters, whereas the outputs $y:[0,t_{\rm end}]\times\mathcal{M}\rightarrow\mathbbm{R}^{n_{\rm out}}$ vary with respect to the parameters. We consider a single output ($n_{\rm out}=1$) without loss of generality.

We assume that the system (13) is either an ODE for all $\mu\in\mathcal{M}$ or a DAE for all $\mu\in\mathcal{M}$. Let the system be asymptotically stable for each parameter with respect to Definition 3. Furthermore, initial value problems

[TABLE]

are considered including a function $x_{0}:\mathcal{M}\rightarrow\mathbbm{R}^{n}$. The initial values may be independent of the parameters. The initial values have to be consistent in the case of DAEs.

We suppose that the parameters in the linear dynamical system (13) are affected by uncertainties. A common approach is to substitute the parameters by independent random variables $\mu:\Omega\rightarrow\mathcal{M}$ on a probability space $(\Omega,\mathcal{F},P)$ with event space $\Omega$, sigma-algebra $\mathcal{F}$ and probability measure $P$, see [37, 39]. We apply traditional probability distributions like uniform, Gaussian, beta, etc. Hence a joint probability density function $\rho:\mathcal{M}\rightarrow\mathbbm{R}$ is available. This approach yields a stochastic model.

A measurable function $f:\mathcal{M}\rightarrow\mathbbm{R}$ depending on the random variables exhibits the expected value

[TABLE]

provided that the integral is finite. The expected value (15) implies the inner product

[TABLE]

for two functions in the Hilbert space

[TABLE]

The associated norm is $\|f\|_{\mathcal{L}^{2}(\mathcal{M},\rho)}=\sqrt{\langle f,f\rangle}$ as usual.

3.2 Polynomial chaos expansions

In most cases, a complete orthogonal basis $(\Phi_{i})_{i\in\mathbbm{N}}$ of polynomials $\Phi_{i}:\mathcal{M}\rightarrow\mathbbm{R}$ exists. These multivariate polynomials are the products

[TABLE]

where $(\Psi_{\ell}^{(j)})_{\ell\in\mathbbm{N}_{0}}$ is the family of univariate orthogonal polynomials with respect to the $j$th random variable. The degree of $\Psi_{\ell}^{(j)}$ is exactly $\ell\geq 0$. Let this basis also be normalized. Each traditional probability distribution implies its own family of orthogonal basis polynomials, see [39]. The basis is complete in the case of uniform, beta, and Gaussian distribution, for example. Yet the polynomials do not span the complete Hilbert space (17) in the case of a log-normal distribution, see [8].

We assume that the QoI of the system (13) is in the space (17) for each time point. If the parameter domain $\mathcal{M}$ is compact, then the continuity of the QoI is sufficient for belonging to (17) pointwise in time. If the parameter domain $\mathcal{M}$ is unbounded, then integrability conditions have to be satisfied with respect to the probability distribution. Consequently, the QoI can be expanded into the series

[TABLE]

with coefficient functions $w_{i}:[0,t_{\rm end}]\rightarrow\mathbbm{R}$, which is called a (generalized) polynomial chaos expansion (PCE). The series (19) converges in the norm of (17) pointwise in time.

Likewise, the state variables exhibit the PCE

[TABLE]

with coefficient functions $v_{i}:[0,t_{\rm end}]\rightarrow\mathbbm{R}^{n}$, provided that each state variable is in the Hilbert space (17).

We assume that the first basis polynomial is the unique constant polynomial $\Phi_{1}\equiv 1$. The orthonormality $\langle\Phi_{i},\Phi_{j}\rangle=\delta_{ij}$ of the basis functions imply the formulas

[TABLE]

for the expected value and the variance of the QoI in each time point.

3.3 Stochastic Galerkin method

We truncate the series (19),(20) to a finite sum. Typically, all basis polynomials up to some total degree $d$ are included. There is a bijective mapping between the integers $i$ and the multi-indices $i_{1},\ldots,i_{q}$ with the number $q$ of random parameters. In view of (18), we obtain the finite index set $\{i\in\mathbbm{N}\,:\,i_{1}+i_{2}+\cdots+i_{q-1}+i_{q}\leq d\}$. The cardinality of this index set is

[TABLE]

If the number of random variables is large, then the number of basis polynomials becomes huge even for moderate degrees, say $d=3$. Figure 1 illustrates the growth of the number of basis polynomials.

Without loss of generality, the truncated series read as

[TABLE]

Inserting the approximations (21) into the linear dynamical system (13) yields a residual. The Galerkin approach requires that the residual is orthogonal to the space spanned by the basis polynomials $\Phi_{1},\ldots,\Phi_{m}$ with respect to the inner product (16). Basic calculations produce a larger coupled linear dynamical system

[TABLE]

for $\hat{v}=(\hat{v}_{1}^{\top},\ldots,\hat{v}_{m}^{\top})^{\top}$ and $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{m})^{\top}$. Hence we obtain a linear dynamical system of dimension $\hat{n}=mn$ with $m$ outputs. The matrices exhibit the sizes $\hat{A},\hat{E}\in\mathbbm{R}^{\hat{n}\times\hat{n}}$, $\hat{B}\in\mathbbm{R}^{\hat{n}\times n_{\rm in}}$, and $\hat{C}\in\mathbbm{R}^{m\times\hat{n}}$. To define the matrices, we introduce the auxiliary arrays

[TABLE]

It follows that

[TABLE]

using Kronecker products. Therein, the probabilistic integration (15) operates separately in each component of the matrices. Often it holds that ${\rm rank}(\hat{C})=m$. More details on the stochastic Galerkin method for linear dynamical systems are given in [25, 32].

The mass matrix $\hat{E}$ may be singular even though all matrices $E$ are non-singular due to properties of the spectrum, see [36]. However, this loss of invertibility hardly occurs. Yet the mass matrix $\hat{E}$ may become ill-conditioned. We assume that $\hat{E}$ is always non-singular in the case of ODEs (13). Likewise, the stochastic Galerkin system (22) may be unstable even though the systems (13) are asymptotically stable for all parameters. Academic examples are given in [31]. Again this loss of stability is hardly observed in practice. Furthermore, the passivity of stochastic Galerkin systems was investigated for models of electric circuits in [16].

Initial conditions $\hat{v}(0)=\hat{v}_{0}$ are derived from the initial values (14) of the original dynamical system (13) by an own truncated PCE. If the initial values (14) are identical to zero, then the choice $\hat{v}(0)=0$ is obvious. The approximation of the QoI reads as

[TABLE]

where the outputs $\hat{w}_{1},\ldots,\hat{w}_{m}$ of the stochastic Galerkin system (22) yield the coefficients.

If the entries of the matrices are polynomials depending on $\mu$ in the system (13), then the matrices (24) of the stochastic Galerkin system can be calculated analytically. Consequently, we obtain the matrices exactly (except for round-off errors) independent of the number $q$ of random variables. In contrast, stochastic collocation techniques, which are non-intrusive methods, cf. [37, 39], induce a quadrature error or sampling error. This error typically grows for fixed numbers of collocation points and increasing dimensions $q$. Although the stochastic Galerkin approach represents an intrusive method, the effort of coding the algorithms is not extensive, because just constant matrices have to be specified for linear time-invariant systems.

We prove a property, which will be used later.

Lemma 1**.**

If the matrices $E(\mu)$ are symmetric and positive definite for almost all $\mu\in\mathcal{M}$, then the stochastic Galerkin projection $\hat{E}$ is also symmetric and positive definite.

Proof:

The Galerkin-projected matrix consists of the blocks $\hat{E}_{ij}=\mathbb{E}[\Phi_{i}\Phi_{j}E]$ for $i,j=1,\ldots,m$, see (24). Hence the symmetry is obvious. Let $z=(z_{1}^{\top},\ldots,z_{m}^{\top})^{\top}\in\mathbbm{R}^{mn}$. We obtain

[TABLE]

because the integrand is almost everywhere non-negative in the probabilistic integration (15). The basis functions $(\Phi_{i})_{i\in\mathbbm{N}}$ are linearly independent. Thus $z\neq 0$ implies that the above finite sum is non-zero on a subset $\mathcal{U}\subset\mathcal{M}$ satisfying $P(\{\omega:\mu(\omega)\in\mathcal{U}\})>0$ for the probability measure $P$. It follows that $z^{\top}\hat{E}z>0$ for $z\neq 0$. $\Box$

Likewise, this relation applies to the case of negative definite matrices.

4 Stability preservation

We investigate three strategies to preserve the asymptotic stability in an MOR of a stochastic Galerkin system. Figure 2 illustrates different possibilities.

4.1 Transformation of stochastic Galerkin system

This approach follows the steps (c),(b),(f) in Figure 2. If the linear dynamical systems (13) are asymptotically stable for all parameters, then the stochastic Galerkin system (22) is usually asymptotically stable. However, if the system (22) is not dissipative, then stability may be lost in some MOR methods by step (e). In this case, we can transform the system to a dissipative form as demonstrated in Section 2.5 and Section 2.6. Consequently, the reduced system is stable due to Theorem 2. Remark that we do not have to calculate the transformed matrices of the high-dimensional system explicitly. Just an appropriate projection matrix has to be determined via (10).

In this approach, the critical part is the solution of the high-dimensional Lyapunov equation (11). A direct method of linear algebra would require $\mathcal{O}(m^{3}n^{3})$ operations. Thus we are restricted to approximate methods or iteration schemes. The possibilities are listed in Section 2.6.

4.2 Transformation of parameter-dependent system

Now the succession (a),(d),(f) is performed with respect to Figure 2. Consequently, we transform the linear dynamical systems (13) first.

4.2.1 Transformation

The stochastic Galerkin system is not always asymptotically stable, even if all parameter-dependent systems (13) are asymptotically stable. However, we obtain a positive result in the case of dissipativity.

Theorem 4**.**

If the linear dynamical systems (13) are dissipative for almost all $\mu\in\mathcal{M}$, then the stochastic Galerkin system (22) is also dissipative. Consequently, a Galerkin-type reduction of (22) yields an asymptotically stable system.

Proof:

Lemma 1 implies that the mass matrix $\hat{E}$ is symmetric and positive definite. The matrix $\hat{A}+\hat{A}^{\top}$ is the Galerkin projection of $A(\cdot)+A(\cdot)^{\top}$ due to

[TABLE]

see (24). Since $A(\mu)+A(\mu)^{\top}$ is negative definite for almost all $\mu\in\mathcal{M}$, Lemma 1 shows that $\hat{A}+\hat{A}^{\top}$ is negative definite. Hence both conditions of Definition 6 are satisfied. Theorem 2 guarantees the stability of a reduced system. $\Box$

If the original systems (13) are not dissipative for almost all $\mu\in\mathcal{M}$, then we transform the systems appropriately. The following transformation has to be applied for almost all $\mu$ simultaneously (even if some system is already dissipative) to preserve continuity and smoothness of the functions.

We obtain the parameter-dependent Lyapunov equation

[TABLE]

for $\mu\in\mathcal{M}$. Although the matrix $F$ may depend on the parameters, we choose a constant matrix, because no parameter-aware strategy with benefits is known yet. The Lyapunov equations (26) yield a unique family $M(\mu)$ of symmetric positive definite matrices.

Furthermore, the stochastic Galerkin system (22) obtained by step (c) is not equivalent to a stochastic Galerkin system obtained by steps (a),(d). On the one hand, the stochastic Galerkin method is invariant with respect to transformations by constant non-singular matrices. On the other hand, we use the parameter-dependent matrix $E(\mu)^{\top}M(\mu)$ in the transformation (a) to the dissipative form (9).

4.2.2 Polynomial system matrices

Now we assume that the matrices $A(\mu),B(\mu),E(\mu)$ of the system (13) involve only polynomials in the variable $\mu$, which is often given in practice. Thus the matrices $\hat{A},\hat{B},\hat{E}$ of the stochastic Galerkin system (22) can be calculated analytically in the case of traditional probability distributions. We do not consider the output matrix $C(\mu)$, because it is not transformed. However, the matrix-valued function $M(\mu)$ satisfying (26) consists of rational functions in the variable $\mu$. In the case of low dimensions $n$, the function $M(\mu)$ can be calculated explicitly by a computer algebra software. Alternatively, the solution $M(\mu)$ can be evaluated for a finite set of parameters $\mu$.

We require the Galerkin projection of the transformed matrices

[TABLE]

in the dissipative system (9). The entries of these transformed matrices are rational functions of $\mu$. A quadrature rule yields numerical approximations $\tilde{A},\tilde{B},\tilde{E}$ of the matrices $\hat{A}^{\prime},\hat{B}^{\prime},\hat{C}^{\prime}$. Let $\{\mu_{1},\ldots,\mu_{k}\}\subset\mathcal{M}$ be the nodes and $\{\gamma_{1},\ldots,\gamma_{k}\}\subset\mathbbm{R}$ be the weights. The approximation of the blocks in $\hat{A}^{\prime}$ reads as

[TABLE]

for $i,j=1,\ldots,m$. Using the auxiliary array (23), the approximation becomes

[TABLE]

Likewise, the quadrature yields $\tilde{B}$ and $\tilde{E}$. The computational effort is dominated by determining $k$ numerical solutions of the Lyapunov equations (26).

Lemma 2**.**

If all weights are positive, then a quadrature of kind (27) yields approximations $\tilde{A},\tilde{E}$, where $\tilde{E}$ is symmetric positive semi-definite and $\tilde{A}+\tilde{A}^{\top}$ is negative semi-definite.

Proof:

The symmetry of $\tilde{E}$ is obvious. Let $z\in\mathbbm{R}^{mn}\backslash\{0\}$ with $z=(z_{1}^{\top},\ldots,z_{m}^{\top})^{\top}$ and

[TABLE]

We obtain using the notation (28)

[TABLE]

Since $E^{\prime}(\mu_{\ell})$ is positive definite and the weights $\gamma_{\ell}$ are positive for each $\ell$, all terms of the sum are non-negative. Hence $\tilde{E}$ is positive semi-definite. The negative semi-definiteness of $\tilde{A}+\tilde{A}^{\top}$ is concluded by the same treatment. $\Box$

In addition, it is very likely that one or more terms are positive in the above sum. Thus we assume that these matrices are definite. It follows that the approximate Galerkin system has the dissipativity condition of Definition 6.

Concerning the positivity of the weights, we address three classes of multivariate quadrature schemes or sampling methods, which can all be found in [37]:

• i)

tensor-product formulas: Univariate quadrature rules often exhibit exclusively positive weights like Gaussian quadrature, for example. The tensor-product formulas inherit the positivity, because their weights are the products of the weights in the univariate case.

• ii)

sparse grid quadrature: Often both positive and negative weights occur. Sparse grids with purely positive weights typically require more nodes for the same accuracy.

• iii)

(quasi) Monte-Carlo methods: In these sampling techniques, the weights are $\gamma_{\ell}=\frac{1}{k}$ for all $\ell=1,\ldots,k$ and thus positive.

4.2.3 Properties

The strategy of this section looks appealing in the case of low dimensions $n$, because the Lyapunov equations can be solved cheap. However, the number of random parameters is typically large in this case, because otherwise the stochastic Galerkin system is not high-dimensional and an MOR is obsolete. A large number of random variables implies a quadrature in a high-dimensional space. Using a computationally feasible number of nodes, the quadrature error may still be too large such that the exact solution of the approximate Galerkin system yields poor approximations (25) of the random-dependent QoI. Therefore, the above approach is critical.

Furthermore, there is a major loss of sparsity in this technique. Let the entries in the matrices be polynomials of low degrees depending on $\mu$ in the system (13). Even if the matrices are dense, then the stochastic Galerkin system (22) exhibits sparse matrices due to the orthogonality of the basis polynomials. However, the dense matrices of the transformed systems (9) are rational functions depending on $\mu$. Inner products (16) of their entries and orthogonal polynomials are non-zero and thus the matrices of (22) become dense in the Galerkin projection. In contrast to the approach of Section 4.1, we have to calculate the matrices of the alternative Galerkin system explicitly to determine a projection matrix $V$ of the MOR.

4.3 Transformation using reference parameter

Again the steps (c),(b),(f) are considered in Figure 2, where the transformation is done different from Section 4.1. We derive an additional technique to decrease the computational effort and to omit quadrature errors. A single reference parameter $\mu^{*}\in\mathcal{M}$ is selected like the mean value $\mu^{*}=\mathbb{E}[\mu]$, for example. We directly solve the Lyapunov equation (26) only for $\mu^{*}$ using some matrix $F$. The solution $M^{*}=M(\mu^{*})\in\mathbbm{R}^{n\times n}$ is symmetric and positive definite. We define the larger transformation matrix

[TABLE]

using the identity matrix $I_{m}\in\mathbbm{R}^{m\times m}$ and the Kronecker product. Obviously, the matrix (29) is symmetric and positive definite again. We employ this matrix to transform the stochastic Galerkin system (22) into the form (9). Since the matrix (29) is block-diagonal, matrix-matrix multiplications with $\hat{M}$ are cheap. We do not need to compute the transformed system matrices. Alternatively, we directly compute the projection matrix (10), where a projection matrix $V$ is determined by the original Galerkin system.

Any random variables can be decomposed into the form $\mu(\omega)=\mu^{*}+\Delta\mu(\omega)$. We consider the random variables

[TABLE]

with the real parameter $\theta\in[0,1]$.

Theorem 5**.**

Let the random variables in a linear system of ODEs (13) be of the form (30). There is a positive constant $\theta_{0}\in(0,1]$ such that the transformation of a stochastic Galerkin system (22) using the matrix (29) yields a dissipative system of type (9) for all $\theta\in[0,\theta_{0}]$.

Proof:

The transformed mass matrix $\hat{E}^{\top}\hat{M}\hat{E}$ is symmetric and positive definite for all $\theta$, since just the non-singularity of the original mass matrix $\hat{E}$ is required. In the limit, we obtain

[TABLE]

due to the orthogonality of the basis polynomials. Hence the matrix becomes block-diagonal with identical blocks. The factor $M^{*}$ satisfies the Lyapunov equation (26) for the chosen positive definite matrix $F$. The symmetric part of the matrix (31) is negative definite, because it holds that

[TABLE]

The two conditions of Definition 6 are satisfied and thus the system is dissipative in the limit. Since the limit (31) is reached continuously with respect to the parameter $\theta$, it follows that a sufficiently small perturbation of $\theta=0$ still yields matrices with the required properties. $\Box$

Theorem 5 shows that the stochastic Galerkin system becomes dissipative for all sufficiently small $\theta>0$. However, this implication does not guarantee that the system is dissipative for $\theta=1$, which reproduces our desired choice of random variables in (30). Nevertheless, the computational effort is low such that it is worth to try this approach. Even if this transformed stochastic Galerkin system is not dissipative, a loss of stability may happen less often in an MOR.

4.4 Nonlinear dynamical systems

We outline a stabilization concept for nonlinear dynamical systems. Let an autonomous system of ODEs be given in the form

[TABLE]

with a parameter-dependent mass matrix $E\in\mathbbm{R}^{n\times n}$ and a nonlinear smooth right-hand side $f:\mathbbm{R}^{n}\times\mathcal{M}\rightarrow\mathbbm{R}^{n}$. The QoI $y$ is defined by the linear or nonlinear function $g:\mathbbm{R}^{n}\times\mathcal{M}\rightarrow\mathbbm{R}^{n_{\rm out}}$. Let $n_{\rm out}=1$. We assume that a family of asymptotically stable stationary solutions $x^{*}:\mathcal{M}\rightarrow\mathbbm{R}^{n}$ exists, i.e., $f(x^{*}(\mu),\mu)=0$ for all $\mu\in\mathcal{M}$. The definition of stable stationary solutions can be found in [34, p. 22]. As in [31], the system (32) is replaced by the equivalent system

[TABLE]

which exhibits the constant asymptotically stable stationary solution $\tilde{x}^{*}=0$. Let $\tilde{J}(\mu)\in\mathbbm{R}^{n\times n}$ be the Jacobian matrix of $\tilde{f}$ evaluated at $\tilde{x}^{*}=0$ and $\mu\in\mathcal{M}$. The stochastic Galerkin method yields a larger nonlinear system of ODEs

[TABLE]

including the mass matrix $\hat{E}$ from (24). The definition of the nonlinear function $\hat{f}:\mathbbm{R}^{mn}\rightarrow\mathbbm{R}^{mn}$ is given in [30, 31]. This function inherits the smoothness of $f$. Now $\hat{v}^{*}=0$ is a stationary solution of (34). Although this stationary solution may be unstable, this loss of stability hardly occurs in practice. Thus we assume that the stationary solution is asymptotically stable again. Let $\hat{J}\in\mathbbm{R}^{mn\times mn}$ be the Jacobian matrix of $\hat{f}$ evaluated at $\hat{v}^{*}=0$. It follows that the spectral abscissa of the matrix pencil $(\hat{E},\hat{J})$ is negative, see Definition 1.

A Galekin-type projection-based MOR yields a small dynamical system

[TABLE]

with a projection matrix $V$ and the mass matrix $\bar{E}=V^{\top}\hat{E}V$. The reduced system owns the stationary solution $\bar{v}^{*}=0$, which may be unstable.

Stability-preserving techniques can be derived. We consider the method of Section 4.1, where the high-dimensional Lyapunov equation (11) is solved including the mass matrix $\hat{E}$ from (34) and $A=\hat{J}$. The Galerkin system (34) is transformed and reduced to (35), while the asymptotic stability of the stationary solution is guaranteed as shown in [29, p. 35] for general nonlinear implicit systems of ODEs. Concerning the stability-preserving approach of Section 4.2, parameter-dependent Lyapunov equations (26) are solved including $E(\mu)$ from (32) and $A(\mu)=\tilde{J}(\mu)$ with the Jacobian matrix associated to the stationary solution of (33). This technique was used for explicit ODEs in [31]. The application of the stability-preserving method from Section 4.3 is straightforward now.

5 Differential-algebraic equations

If the linear dynamical system (1) features a singular mass matrix $E$, then a system of DAEs is given. MOR methods are also available for DAEs, see [6]. Yet the Lyapunov equations (11) do not have a solution, which is crucial in our stability-preserving technique. Let $x\in\mathbbm{R}^{n}$. It follows that

[TABLE]

Choosing an $x\neq 0$ in the kernel of the matrix $E$ implies $Ex=0$ and thus a contradiction to the definiteness of the matrix $F$ appears. We require an alternative strategy now.

5.1 Regularization

We apply a regularization of an asymptotically stable DAE system (1), which was also used in [20]. The asymptotic stability guarantees a regular matrix pencil as specified by Definition 2. The regularized system matrices read as

[TABLE]

introducing parameters $\alpha,\beta>0$. The matrix $E_{\rm reg}$ is non-singular for all $\alpha>0$. Choosing $\alpha=\beta^{2}$, it follows that the linear dynamical system (1) with $A_{\rm reg},E_{\rm reg}$ is asymptotically stable for all sufficiently small parameters $\beta$, see [20]. An advantage of this regularization technique is that the sparsity pattern of $sE-A$ coincides with the sparsity pattern of $sE_{\rm reg}-A_{\rm reg}$ for each $s\in\mathbbm{C}$. Hence no loss of sparsity happens in the relevant operations.

We also choose a small parameter $\beta$ (and $\alpha=\beta^{2}$) to ensure that the difference between the original DAE system and the regularized ODE system is small. Assuming that the DAE system together with the defined outputs has a transfer function with finite $\mathcal{H}_{2}$-norm (3), error bounds were derived with respect to this norm in [28].

5.2 Stochastic Galerkin projection

The steps (a)-(d) refer to the flowchart in Figure 2. The following two approaches are equivalent provided that the same parameters $\alpha,\beta$ are chosen in a regularization (36):

• i)

Regularize the parameter-dependent system (13) and then project the systems of ODEs in the stochastic Galerkin method.

• ii)

Project the parameter-dependent DAE system (13) in the stochastic Galerkin method and then regularize the Galerkin system (22).

In both cases, the stochastic Galerkin projection (c) yields the same matrices $\hat{A}$ and $\hat{E}$ for the system (22). The mass matrix $\hat{E}$ is non-singular due to the regularization.

The transformations (a) and (b) to a dissipative representation are done only for a regularized system, since a system of ODEs is required for each Lyapunov equation. The transformation (b) involves the matrices from (i) or (ii). However, in the transformation (a), we have to regularize the systems of DAEs (13) immediately. Furthermore, the drawbacks of the succession (a),(d), see Section 4.2.3, also apply in this case.

In the transformation (b), the approximate solution of the high-dimensional Lyapunov equations (11) may become critical. Small regularization parameters $\alpha$ nearly coincide with the singular case. This problem is less pronounced in direct methods to solve the Lyapunov equations.

We can also employ the technique from Section 4.3 in this context. The linear dynamical system (13) is regularized for a reference parameter $\mu^{*}$ and some $\alpha,\beta>0$. Solving the Lyapunov equation (26) including $\mu^{*}$ yields the high-dimensional transformation matrix (29). Now the regularized Galerkin system from strategy (ii), where the same values $\alpha,\beta$ are used, is transformed based on (29). Again the dissipativity of the transformed representation cannot be guaranteed a priori.

6 Illustrative examples

We investigate a system of ODEs as well as a system of DAEs in this section. We computed on a FUJITSU Esprimo P920 Intel(R) Core(TM) i5-4570 CPU with 3.20 GHz (4 cores) and operation system Microsoft Windows 7. The software package MATLAB [18] (version R2018a) was used for all computations.

6.1 Mass-spring-damper system

Figure 3 depicts a mechanical configuration, which consists of 5 masses, 7 springs and 5 dampers. The single input is the excitation at the bottom spring, whereas the single output is the position of the top mass. A mathematical modeling yields a linear dynamical system (13) of $n=10$ ODEs of first order including $q=17$ parameters. The system is asymptotically stable for all positive parameters. The system matrices are affine-linear functions (polynomials of degree one) of the parameters. The Bode plot of the system is shown for a specific choice of the parameters in Figure 4. This system represents an extension of a test example used in [14], where 4 masses were included.

In the stochastic modeling, we replace the parameters by independent uniformly distributed random variables, which vary 10% around their mean values. The mean values are the constant parameters from above. In the truncated orthogonal expansion (21), we include all multivariate Legendre polynomials up to degree three. We obtain $m=1140$ basis functions. The stochastic Galerkin system (22) is calculated exactly except for round-off errors. Table 1 illustrates the properties of this example. The system is asymptotically stable. The mass matrix is symmetric and positive definite. Yet the system is not dissipative.

We employ the one-sided Arnoldi method, see [1], to perform an MOR of the stochastic Galerkin system (22). The single real expansion point $s=0.7$ is used in this Krylov subspace technique. The reduced systems are arranged for dimensions $r=1,\ldots,100$. It follows that just 38 ROMs are stable.

Now we investigate the stabilization techniques from Section 4. The following cases are discussed:

• i)

transformation of stochastic Galerkin system (Section 4.1),

• ii)

transformation of original systems (Section 4.2), and

• iii)

transformation based on reference parameter (Section 4.3).

In (i) and (iii), we reuse the projection matrix $V$ determined for the original stochastic Galerkin system, because the Arnoldi method is invariant with respect to a basis transformation in the image space. In (ii), we repeat the Arnoldi algorithm for the alternative stochastic Galerkin system, because a smaller reduction error is expected.

We choose the identity matrix ($F=I$) as degree of freedom in each Lyapunov equation. In all three approaches, the stability is achieved for each ROM. Figure 5 illustrates the relative errors with respect to the $\mathcal{H}_{2}$-norm (3), i.e.,

[TABLE]

with the transfer functions $H,\bar{H}$ of FOM and ROM, respectively.

In the technique (i), we use the method from [28], where an integral is discretized by a quadrature rule in the frequency domain. Therein, the solution of the high-dimensional Lyapunov equation is not computed but the associated projection matrix (10). The computation work is characterized by an $LU$-decomposition of the system matrix ${\rm i}\omega_{j}\hat{E}-\hat{A}$ in each node $\omega_{j}$ for the angular frequency $\omega$. We employ the (univariate) Gauss-Legendre quadrature. Table 2 shows the number of stable ROMs for different numbers of nodes. We observe that 40 nodes are sufficient to stabilize all reduced systems. Figure 5 (i) depicts the relative $\mathcal{H}_{2}$-errors (37) in this case. The error of the MOR exhibits the same magnitude as in the original stochastic Galerkin system.

In the technique (ii), we use a sparse grid quadrature of Smolyak-type with level 3 based on the Clenshaw-Curtis rule. The scheme exhibits 7209 nodes in the 17-dimensional space and negative weights arise. Table 3 illustrates the loss of sparsity in this approach. Again all ROMs become asymptotically stable. The relative errors (37) of the MOR are demonstrated in Figure 5 (ii). The errors between the ROMs and the stochastic Galerkin projection of the transformed systems, which we call the inherent errors, decay for increasing dimensions. However, the errors between these ROMs and the stochastic Galerkin system of the untransformed systems (13), which we call the comparison, are large and stagnate. This property indicates that the quadrature error is too large, even though a high number of nodes is used, because the original Galerkin system can be considered sufficiently accurate.

In the technique (iii), we define the mean value of the random variables as reference parameter. We calculate the matrices of the transformed Galerkin system to analyze the definiteness. Figure 6 depicts the maximum eigenvalue of the symmetric part of $\hat{E}^{\top}\hat{M}\hat{A}$ with $\hat{M}$ from (29) in dependence on the parameter $\theta$ from (30). We observe that the negative definiteness is lost for $\theta\geq\theta_{0}\approx 0.04$. Although our relevant case $\theta=1$ is not dissipative, the stability is preserved in all ROMs. Moreover, the error of the MOR is not compromised.

We also measured computation times illustrated by Table 4. The ROMs of dimension $r=100$ are determined using the Arnoldi method and the three stability-preserving techniques. The computation time of a stabilization technique consists of all steps to calculate the reduced matrices (5) except for the determination of the projection matrix $V$, which is done by the Arnoldi method. In the technique (i), $k=40$ nodes are used in the Gaussian quadrature. On the one hand, the technique (ii) causes a relatively large computational effort due to the high number of nodes in the sparse grid quadrature. Also the Arnoldi algorithm is more expensive due to the loss of sparsity in the matrices. On the other hand, the technique (iii) requires a negligible computation work, which shows that it is worth trying this approach.

6.2 Band-pass filter

We examine the electric circuit of a band-pass filter shown in Figure 7. A single input voltage is supplied and a single output voltage drops at a load conductance. Modified nodal analysis [12] yields a linear system of DAEs of dimension $n=23$. This DAE system exhibits the index one. The physical parameters are 7 capacitances, 7 inductances and 9 conductances ($q=23$). Figure 8 depicts the Bode plot of the system for a constant selection of the parameters. Furthermore, the system matrices are affine-linear functions of the parameters.

In the stochastic modeling, we introduce independent uniformly distributed random variables varying 20% around their mean values. The truncated series (21) include all basis polynomials up to degree two, i.e., $m=300$. The stochastic Galerkin system (22) consists of linear DAEs, which have a finite $\mathcal{H}_{2}$-norm (3) for the defined outputs.

In the regularization (36), we choose the parameters $\alpha=10^{-10}$ and $\beta=10^{-5}$. Table 5 shows the properties of the regularized system. The system is asymptotically stable. Both $\hat{E}$ and $\hat{E}_{\rm reg}$ are not symmetric. Thus the spectral abscissa, see Definition 1, of the symmetric part of $\hat{A}_{\rm reg}$ is irrelevant. Multiplication by $\hat{E}_{\rm reg}^{\top}$ from the left yields a system with symmetric positive definite mass matrix, which is still not dissipative.

The one-sided Arnoldi method yields the projection matrices of the MOR for dimensions $r=1,\ldots,100$. Therein, we choose the real number $s=10^{6}$ as expansion point. A loss of stability occurs in the reduction of both the DAE system and the regularized system, as shown in Table 6. We apply the stabilization techniques of Section 4.1 and Section 4.3 to the regularized system. The Lyapunov equations are solved using the identity matrix as input matrix. We solve the Lyapunov equations by a direct method of linear algebra. Thus a critical behavior for small regularization parameters, as mentioned in Section 5.2, is avoided. The matrices of a transformed Galerkin system are never computed explicitly. All ROMs become asymptotically stable in both approaches.

Finally, we compute approximations of the $\mathcal{H}_{2}$-norms for the difference between the DAE system and the reduced systems. Figure 9 illustrates the relative errors (37) of the reductions. We recognize that the errors are nearly identical in the two stabilization techniques. The errors stagnate for reduced dimensions $r>80$ in (ii)-(iv), because the total error is dominated by the error of the regularization in this part. Most important, the stabilization approaches do not compromise the error of the MOR.

7 Conclusions

We examined stability-preserving model order reduction of linear stochastic Galerkin systems using transformations to dissipative forms. Three approaches were analyzed. The transformation of the conventional Galerkin system represents an adequate method. The Galerkin projection of the transformed original systems features severe drawbacks in the case of large numbers of random parameters. The approach using a cheap transformation matrix based on a reference parameter is promising. Although the dissipativity property cannot be guaranteed in the transformed system, the numerical results of test examples demonstrate that preservation of stability is achieved. Moreover, the error of the model order reduction does not increase in this stabilization.