Model order reduction for random nonlinear dynamical systems and low-dimensional representations for their quantities of interest

TL;DR

This paper explores model order reduction techniques, specifically proper orthogonal decomposition, for efficiently approximating quantities of interest in random nonlinear dynamical systems affected by uncertainty.

Contribution

It introduces a projection-based model reduction approach for large stochastic dynamical systems, providing error analysis and numerical validation.

Findings

Reduced-order models significantly decrease computational complexity.

Error bounds for low-dimensional approximations are established.

Numerical examples demonstrate effectiveness of the method.

Abstract

We examine nonlinear dynamical systems of ordinary differential equations or differential algebraic equations. In an uncertainty quantification, physical parameters are replaced by random variables. The inner variables as well as a quantity of interest are expanded into series with orthogonal basis functions like the polynomial chaos expansions, for example. On the one hand, the stochastic Galerkin method yields a large coupled dynamical system. On the other hand, a stochastic collocation method, which uses a quadrature rule or a sampling scheme, can be written in the form of a large weakly coupled dynamical system. We apply projection-based methods of nonlinear model order reduction to the large systems. A reduced-order model implies a low-dimensional representation of the quantity of interest. We focus on model order reduction by proper orthogonal decomposition. The error of a best approximation located in a low-dimensional subspace is analysed. We illustrate results of numerical computations for test examples.

Tables

Table 1. Table 1: Maximum differences to reference solution in time intervals [ 0 , T ] 0 𝑇 [0,T] and [ 0 , 3 ​ T ] 0 3 𝑇 [0,3T] for expected value as well as variance of random output voltage in transistor amplifier example.

		expected value		variance
	ROM-dimension	$[0,T]$	$[0,3T]$	$[0,T]$	$[0,3T]$
Galerkin	FOM	$4.1\cdot {10}^{-3}$	$4.1\cdot {10}^{-3}$	$4.0\cdot {10}^{-4}$	$4.0\cdot {10}^{-4}$
	$25$	$4.5\cdot {10}^{-2}$	$4.5\cdot {10}^{-2}$	$8.3\cdot {10}^{-3}$	$8.3\cdot {10}^{-3}$
	$50$	$4.5\cdot {10}^{-2}$	$4.5\cdot {10}^{-2}$	$8.3\cdot {10}^{-3}$	$8.3\cdot {10}^{-3}$
Collocation	FOM	$5.4\cdot {10}^{-3}$	$5.4\cdot {10}^{-3}$	$3.5\cdot {10}^{-4}$	$3.5\cdot {10}^{-4}$
	$10$	$7.3\cdot {10}^{-3}$	$4.3\cdot {10}^{-2}$	$1.1\cdot {10}^{-2}$	$5.6\cdot {10}^{-2}$
	$25$	$4.2\cdot {10}^{-3}$	$4.2\cdot {10}^{-3}$	$3.4\cdot {10}^{-4}$	$3.4\cdot {10}^{-4}$
	$50$	$4.2\cdot {10}^{-3}$	$5.6\cdot {10}^{-3}$	$3.3\cdot {10}^{-4}$	$1.9\cdot {10}^{-3}$

Equations

$$\begin{array}[]{rcl}E(p)\dot{x}(t,p)&=&A(p)x(t,p)+F(x(t,p),p)+B(p)u(t)\\[4.30554pt] y(t,p)&=&C(p)x(t,p)\\ \end{array}$$

$$x(t_0,p)=x_0(p)\text{mbox}forp\in \Pi$$

$$\mathbb{E}\left[f\right]:={\int }_{\Pi }f(p)\rho (p)\text{mbox}dp$$

$${\mathcal{L}}^2(\Pi ,\rho ):=\left\{f:\Pi \to \text{mathbbm}R:f\text{mbox}measurableand\mathbb{E}\left[f^2\right]<\infty \right\}$$

$$<f,g>:=\mathbb{E}\left[fg\right]={\int }_{\Pi }f(p)g(p)\rho (p)\text{mbox}dp\text{mbox}forf,g\in {\mathcal{L}}^2(\Pi ,\rho ).$$

$${\Vert f\Vert }_{{\mathcal{L}}^2(\Pi ,\rho )}:=\sqrt{<f,f>}.$$

$$x(t,p)=\sum _{i=1}^{\infty }{v}_i(t){\Phi }_i(p)\text{mbox}andy(t,p)=\sum _{i=1}^{\infty }{w}_i(t){\Phi }_i(p),$$

$$v_{i,j}(t)=<x_j(t,\cdot ),{\Phi }_i(\cdot )>\text{mbox}and{w}_i(t)=<y(t,\cdot ),{\Phi }_i(\cdot )>,$$

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

$${\Vert y(t,\cdot )-y^{(m)}(t,\cdot )\Vert }_{{\mathcal{L}}^2(\Pi ,\rho )}=\sqrt{\sum _{i=m+1}^{\infty }{w}_i(t{)}^2}$$

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

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

$$\bar{y}(t,p)=\sum _{i=1}^r{\bar{w}}_i(t){\Psi }_i(p)$$

$$\mathrm{span}\{{\Psi }_1,\dots ,{\Psi }_r\}\subset \mathrm{span}\{{\Phi }_1,\dots ,{\Phi }_m\}$$

$$\begin{array}[]{rcl}\left\|y(t,\cdot)-\bar{y}(t,\cdot)\right\|_{\mathcal{L}^{2}(\Pi,\rho)}&\leq&\left\|y(t,\cdot)-y^{(m)}(t,\cdot)\right\|_{\mathcal{L}^{2}(\Pi,\rho)}\\[4.30554pt] &&\mbox{}+\left\|y^{(m)}(t,\cdot)-\hat{y}(t,\cdot)\right\|_{\mathcal{L}^{2}(\Pi,\rho)}\\[4.30554pt] &&\mbox{}+\left\|\hat{y}(t,\cdot)-\bar{y}(t,\cdot)\right\|_{\mathcal{L}^{2}(\Pi,\rho)}\\ \end{array}$$

$$\begin{array}[]{rcl}\hat{E}\dot{\hat{v}}(t)&=&\hat{A}\hat{v}(t)+\hat{F}(\hat{v}(t))+\hat{B}u(t)\\[4.30554pt] \hat{w}(t)&=&\hat{C}\hat{v}(t)\\ \end{array}$$

$$\begin{array}[]{ll}\hat{A}=\mathbb{E}[S(\cdot)\otimes A(\cdot)]\in\mathbbm{R}^{mn\times mn},&\hat{B}=\mathbb{E}[s(\cdot)\otimes B(\cdot)]\in\mathbbm{R}^{mn},\\[2.15277pt] \hat{C}=\mathbb{E}[S(\cdot)\otimes C(\cdot)]\in\mathbbm{R}^{m\times mn},&\hat{E}=\mathbb{E}[S(\cdot)\otimes E(\cdot)]\in\mathbbm{R}^{mn\times mn},\\ \end{array}$$

$${\hat{F}}_i(\hat{v}):=\mathbb{E}\left[F\left(\sum _{j=1}^m{\hat{v}}_j{\Phi }_j(\cdot ),\cdot \right){\Phi }_i(\cdot )\right]$$

$${\hat{F}}_i(\hat{v})\approx \sum _{\ell =1}^k{\gamma }_{\ell }F\left(\sum _{j=1}^m{\hat{v}}_j{\Phi }_j(p^{(\ell )}),p^{(\ell )}\right){\Phi }_i(p^{(\ell )})$$

$${\hat{w}}_i(t)=\sum _{\ell =1}^k{\gamma }_{\ell }{\Phi }_i(p^{(\ell )})C(p^{(\ell )})x(t,p^{(\ell )})$$

$$\begin{array}[]{rcl}\hat{E}\dot{\hat{x}}(t)&=&\hat{A}\hat{x}(t)+\hat{F}(\hat{x}(t))+\hat{B}u(t)\\[4.30554pt] \hat{w}(t)&=&\hat{C}\hat{x}(t)\\ \end{array}$$

$$\begin{array}[]{ll}\hat{A}=\begin{pmatrix}A(p^{(1)})&&\\ &\ddots&\\ &&A(p^{(k)})\\ \end{pmatrix}\in\mathbbm{R}^{kn\times kn},&\hat{B}=\begin{pmatrix}B(p^{(1)})\\ \vdots\\ B(p^{(k)})\\ \end{pmatrix}\in\mathbbm{R}^{kn},\\[21.52771pt]$$

\hat{E}=\begin{pmatrix}E(p^{(1)})&&\\ &\ddots&\\ &&E(p^{(k)})\\ \end{pmatrix}\in\mathbbm{R}^{kn\times kn}.&\\ \end{array}$$

$$\hat{F}(\hat{x}(t))=\left(\begin{array}{c}F(\hat{x}(t,p^{(1)}),p^{(1)})\\ ⋮\\ F(\hat{x}(t,p^{(k)}),p^{(k)})\end{array}\right)\in \text{mathbbm}{R}^{kn}.$$

$$\begin{array}[]{rcl}\bar{E}\dot{\bar{v}}(t)&=&\bar{A}\bar{v}(t)+\bar{F}(\bar{v}(t))+\bar{B}u(t)\\[4.30554pt] \bar{w}(t)&=&\bar{C}\bar{v}(t)\\ \end{array}$$

$$\bar{A}=T_{\mathrm{l}}^{\top }\hat{A}{T}_{\mathrm{r}},\bar{B}=T_{\mathrm{l}}^{\top }\hat{B},\bar{C}=\hat{C}{T}_{\mathrm{r}},\bar{E}=T_{\mathrm{l}}^{\top }\hat{A}{T}_{\mathrm{r}},\bar{F}(\bar{v})=T_{\mathrm{l}}^{\top }\hat{F}(T_{\mathrm{r}}\bar{v}).$$

$$V=\left(\hat{v}(t_0),\hat{v}(t_1),\dots ,\hat{v}(t_{\ell -1})\right)\in \text{mathbbm}{R}^{mn\times \ell }.$$

$$V=U\left(\begin{array}{c}\Sigma \\ 0\end{array}\right)Q^{\top }$$

$$T_{\mathrm{r}}=(u_1,u_2,\dots ,u_r)\in \text{mathbbm}{R}^{mn\times r}$$

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

$${\Vert \bar{y}(t,\cdot )-\hat{y}(t,\cdot )\Vert }_{{\mathcal{L}}^2(\Pi ,\rho )}=\sqrt{\sum _{i=1}^m{\left({\bar{w}}_i(t)-{\hat{w}}_i(t)\right)}^2}\text{mbox}foreacht,$$

Full text

**Model order reduction for

random nonlinear dynamical systems and

low-dimensional representations for

their quantities of interest**

Roland Pulch

Institut für Mathematik und Informatik, Universität Greifswald,

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

Email: roland.pulch@uni-greifswald.de

Abstract

[TABLE]

1 Introduction

Mathematical modelling often generates dynamical systems of ordinary differential equations (ODEs) or differential algebraic equations (DAEs) in science and engineering. A quantity of interest (QoI) is defined by the state variables or inner variables. Multiple sources of uncertainties may affect the included parameters like modelling errors and measurement errors, for example. Several techniques have been developed to quantify the effects of those uncertainties in the model predictions, see [21, 40, 43]. A common approach is the substitution of uncertain parameters by random variables. Numerical methods often apply orthogonal expansions with unknown coefficient functions and predetermined random-dependent basis functions. The stochastic Galerkin method (intrusive method) produces a larger coupled dynamical system. A stochastic collocation technique (non-intrusive method) using a quadrature rule or a sampling method can be written as a single large weakly coupled dynamical system, see [31].

A low-dimensional approximation of the random QoI often exists, where only a few basis functions are required for a sufficiently accurate representation. In a specific context, the representation can be interpreted as a sparse approximation. Such low-dimensional representations have been computed by different methods like least angle regression [5], sparse grid quadrature [8], compressed sensing [11] and $\ell_{1}$-minimisation [19, 20]. Stochastic reduced bases were investigated for random linear systems of algebraic equations in [25, 37].

Model order reduction (MOR) becomes favourable in the stochastic Galerkin method and the stochastic collocation approach due to the high dimensionality of the systems. In [32], it was shown that an MOR of the Galerkin system implies a low-dimensional approximation of the QoI in the case of linear dynamical systems. In [33], this strategy was carried over to an MOR of the single auxiliary system of the stochastic collocation method in the linear case. In this paper, we consider nonlinear dynamical systems with random parameters, where the QoI still depends linearly on state variables or inner variables. Thus an MOR generates a low-dimensional representation of the random QoI again.

More precisely, MOR is the tool to identify an approximation of the QoI with a low number of basis functions. We do not investigate a reduction of the number of random parameters, cf. [35].

Several efficient MOR methods exist for linear dynamical systems, see [1, 3, 16, 38]. Error bounds are available by Hardy norms of transfer functions in the frequency domain. Yet MOR for nonlinear dynamical systems still represents a challenging task. Typically, projection-based MOR schemes are applied like the proper orthogonal decomposition (POD), see [1, 22], or the trajectory piecewise linear approach, see [24, 36]. We focus on the POD method, which identifies a low-dimensional representation of the QoI. Furthermore, the representation can be improved by the computation of a best approximation in the low-dimensional subspace. We prove an error bound for this best approximation. In the first place, our aim is the identification of a sufficiently accurate approximation with as few basis functions as possible. In the second place, the MOR methods should decrease the computational effort.

The paper is organised as follows. Section 2 incloses the problem setup. Numerical methods yield large dynamical systems formulated in Section 3. The main part is given by Section 4, where we apply MOR and show error bounds. We demonstrate the results of numerical experiments for two illustrative examples in Section 5.

2 Problem definition

We describe the task consisting in the identification of low-dimensional representations for QoIs from random nonlinear dynamical systems.

2.1 Deterministic model

Let a nonlinear dynamical system be given in the form

[TABLE]

with matrices and functions depending on physical parameters $p\in\Pi\subseteq\mathbbm{R}^{q}$. The sizes of the matrices are $A,E\in\mathbbm{R}^{n\times n}$, $B\in\mathbbm{R}^{n\times n_{\rm in}}$, $C\in\mathbbm{R}^{n_{\rm out}\times n}$. The system involves a nonlinear function $F:\mathbbm{R}^{n}\times\Pi\rightarrow\mathbbm{R}^{n}$. For non-singular matrices $E$, the system consists of ODEs with the state variables $x:[t_{0},t_{\rm end}]\times\Pi\rightarrow\mathbbm{R}^{n}$. For singular matrices $E$, a system of DAEs is given with the inner variables $x$. We consider initial value problems

[TABLE]

with a predetermined function $x_{0}:\Pi\rightarrow\mathbbm{R}^{n}$. In the case of DAEs, the initial values have to satisfy consistency conditions and typically depend on the physical parameters of the system.

An input $u:[t_{0},t_{\rm end}]\rightarrow\mathbbm{R}^{n_{\rm in}}$ is supplied to the system (1). An output $y:[t_{0},t_{\rm end}]\times\Pi\rightarrow\mathbbm{R}^{n_{\rm out}}$ is defined by the state variables or inner variables and the matrix $C$. Without loss of generality, we restrict the analysis to the case of single-input-single-output (SISO) with $n_{\rm in}=n_{\rm out}=1$. On the one hand, the results also apply to dynamical systems with a general nonlinear dependence of the right-hand side on the input. Moreover, the theory is applicable to autonomous dynamical systems. On the other hand, the linear dependence of the output on the state variables or inner variables is essential in this paper.

2.2 Stochastic model

We adopt a common approach in uncertainty quantification (UQ), see [40, 43], for example. Assuming that the parameters $p\in\Pi$ of the system (1) are uncertain, they are replaced by independent random variables $p:\Omega\rightarrow\Pi$ on some probability space $(\Omega,\mathcal{A},P)$ with event space $\Omega$, sigma-algebra $\mathcal{A}$ and probability measure $P$. A joint probability density function $\rho:\Pi\rightarrow\mathbbm{R}$ is available in the case of traditional probability distributions. For a measurable function $f:\Pi\rightarrow\mathbbm{R}$, the expected value reads as

[TABLE]

provided that the integral is finite. The Hilbert space

[TABLE]

is equipped with the inner product

[TABLE]

The accompanying norm reads as

[TABLE]

Concerning the dynamical system (1), we assume that $x_{1}(t,\cdot),\ldots,x_{n}(t,\cdot),y(t,\cdot)\in\mathcal{L}^{2}(\Pi,\rho)$ pointwise for $t\in[t_{0},t_{\rm end}]$.

Now let a complete orthonormal system $(\Phi_{i})_{i\in\mathbb{N}}\subset\mathcal{L}^{2}(\Pi,\rho)$ be given. It holds that $<\Phi_{i},\Phi_{j}>=0$ for $i\neq j$ and $<\Phi_{i},\Phi_{j}>=1$ for $i=j$. We assume that the first basis function is always the constant function $\Phi_{1}\equiv 1$. It follows that the expansions

[TABLE]

where the coefficient functions $v_{i}:[t_{0},t_{\rm end}]\rightarrow\mathbbm{R}^{n}$ and $w_{i}:[t_{0},t_{\rm end}]\rightarrow\mathbbm{R}$ are defined by

[TABLE]

converge in $\mathcal{L}^{2}(\Pi,\rho)$ pointwise for $t$ and $j=1,\ldots,n$ with $x=(x_{1},\ldots,x_{n})^{\top}$.

2.3 Low-dimensional orthogonal representations

Our aim is to identify a low-dimensional approximation using just a few basis functions. Several methods to obtain sparse or low-dimensional representations have already been designed and investigated, see [5, 8, 9, 11, 19, 20].

In practice, the series (5) have to be truncated to finite approximations. A truncation yields

[TABLE]

for an integer $m\geq 1$. We investigate the output of the random dynamical system (1) as QoI. The truncation error reads as

[TABLE]

for each $t$ due to the orthonormality of the basis functions.

Often orthonormal polynomials are chosen as basis functions with respect to the theory of the polynomial chaos (PC) expansions,

see [6, 12, 43]. Therein, the multivariate basis polynomials are just the products of the univariate orthonormal polynomials. Each traditional probability distribution implies its family of univariate orthonormal polynomials, see [41]. To obtain an initial set of basis functions, often all polynomials up to a total degree $d$ are included. Hence the number of basis functions is

[TABLE]

see [43, p. 65]. Even if the total degree is moderate (say $3\leq d\leq 5$), the number of basis polynomials becomes large for high numbers $q$ of random parameters.

A numerical method yields approximations $\hat{w}_{1},\ldots,\hat{w}_{m}$ of the exact coefficient functions in the truncated expansion (7). The investigated QoI reads as

[TABLE]

Assuming a large number $m$, our aim is to identify a sufficiently accurate low-dimensional representation of (9)

[TABLE]

with a new orthonormal basis $\{\Psi_{1},\ldots,\Psi_{r}\}$ satisfying

[TABLE]

and new coefficient functions $\bar{w}_{1},\ldots,\bar{w}_{r}$ for some $r\ll m$. A special case is given by the selection of a subset $\{\Psi_{1},\ldots,\Psi_{r}\}\subset\{\Phi_{1},\ldots,\Phi_{m}\}$. In this situation, the approximation (10) can be interpreted as a sparse representation, see [5].

The error of the complete approach is estimated by

[TABLE]

for each $t$. The upper bound consists of three terms: the truncation error, the error of the numerical method and the additional error of the low-dimensional approximation. We assume that the first and second term are sufficiently small due to choosing a large initial set of basis functions and a sufficiently accurate numerical method. The third term is analysed in Section 4.5.

3 Numerical Techniques

We apply two well-known classes of numerical methods for the computation of the unknown coefficient functions (6) in the orthogonal expansions (5) now.

3.1 Stochastic Galerkin method

The stochastic Galerkin technique (intrusive method) represents a general approach, which can be applied to all types of differential equations including random variables, see [40, 43]. Linear or nonlinear dynamical systems were considered in [2, 27, 28, 29], for example. In our problem, the Galerkin method changes the nonlinear dynamical system (1) into the larger coupled system

[TABLE]

with $\hat{v}=(\hat{v}_{1}^{\top},\ldots,\hat{v}_{m}^{\top})^{\top}\in\mathbbm{R}^{mn}$ and $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{m})^{\top}\in\mathbbm{R}^{m}$. To define the involved matrices, we employ auxiliary arrays: the symmetric matrix $S(p):=(\Phi_{i}(p)\Phi_{j}(p))\in\mathbbm{R}^{m\times m}$ and the column vector $s(p):=(\Phi_{i}(p))\in\mathbbm{R}^{m}$. Now the matrices of the linear parts in (11) read as

[TABLE]

using the expected value (3) componentwise and Kronecker products. The nonlinear part is given by $\hat{F}:=(\hat{F}_{1}^{\top},\ldots,\hat{F}_{m}^{\top})^{\top}$ with $\hat{F}_{i}\in\mathbbm{R}^{n}$ and

[TABLE]

for $i=1,\ldots,m$.

The expected values (12) represent probabilistic integrals, which often cannot be calculated analytically. An exception are functions $F$ in (1) consisting of polynomials with low degrees. We require approximations by quadrature formulas for the case of general nonlinear functions $F$ in (1). A quadrature rule or a sampling method is determined by its nodes $\{p^{(1)},\ldots,p^{(k)}\}\subset\Pi$ and its weights $\{\gamma_{1},\ldots,\gamma_{k}\}\subset\mathbbm{R}$. The associated approximation becomes

[TABLE]

for $i=1,\ldots,m$. Thus the computational effort for one evaluation of $\hat{F}$ comprises $k$ evaluations of $F$, whereas the summation is negligible.

Concerning initial value problems of the system (11), the initial condition (2) of the original system (1) has to be expanded as well. Consistent initial values are required in the case of DAEs.

3.2 Stochastic collocation method

Alternatively, we consider a stochastic collocation (non-intrusive method) employing a quadrature rule or a sampling scheme, see [29, 42]. This approach is determined uniquely by the nodes $\{p^{(1)},\ldots,p^{(k)}\}$ and the weights $\{\gamma_{1},\ldots,\gamma_{k}\}$. The approximations of the coefficient functions belonging to the QoI (9) become

[TABLE]

for $i=1,\ldots,m$. Therein, $C(p^{(\ell)})$ denotes the evaluations of the output matrix from (1) for $\ell=1,\ldots,k$. Thus $k$ separate initial value problems (1),(2) have to be solved for the different nodes, which yields the variables $x$.

We write this approach as a single auxiliary system, which exhibits the form (11) again. Due to a different meaning of the inner variables, we reformulate

[TABLE]

with $\hat{x}(t)=(x(t,p^{(1)})^{\top},\ldots,x(t,p^{(k)})^{\top})^{\top}$. Consequently, the matrices read as

[TABLE]

The matrix $\hat{C}\in\mathbbm{R}^{m\times kn}$ follows from the formula (14). The meaning of the outputs is the same as in the stochastic Galerkin system (11). The nonlinear part becomes

[TABLE]

The initial values follow from (2), i.e., $x(t_{0},p^{(\ell)})=x_{0}(p^{(\ell)})$ for $\ell=1,\ldots,k$. This auxiliary system was already derived and applied for linear dynamical systems in [30, 31]. A similar single system was constructed in the case of a non-intrusive approach for dynamical systems with both stochastic noise and random parameters in [26].

In the stochastic collocation method, the system (15) is large and weakly coupled, because it consists of the separate original systems (1), which are connected just by the supply of the same input and the definition of the outputs.

4 Model order reduction

Although the original dynamical system (1) may be small or medium-sized, the dynamical systems (11) and (15) from the numerical methods are large. Thus they represent excellent candidates for an MOR. Linear stochastic Galerkin systems were reduced successfully in [23, 30, 34, 35, 44].

4.1 Projection-based model order reduction

The nonlinear dynamical systems (11) and (15) represent non-parametric formulations. Without loss of generality, we consider the system (11) as full-order model (FOM). MOR yields a reduced-order model (ROM) in form of a smaller nonlinear dynamical system

[TABLE]

with matrices $\bar{A},\bar{E}\in\mathbbm{R}^{r\times r}$, $\bar{B}\in\mathbbm{R}^{r}$, $\bar{C}\in\mathbbm{R}^{m\times r}$ and a function $\bar{F}:\mathbbm{R}^{r}\rightarrow\mathbbm{R}^{r}$. The state variables or inner variables are $\bar{v}:[t_{0},t_{\rm end}]\rightarrow\mathbbm{R}^{r}$, which require initial conditions $\bar{v}(t_{0})=\bar{v}_{0}$. The input of (16) coincides with the input of (11). The aim is to achieve outputs satisfying $\bar{w}(t)\approx\hat{w}(t)$ for all $t$.

In projection-based MOR, see [1, 13], projection matrices $T_{\rm l},T_{\rm r}\in\mathbbm{R}^{mn\times r}$ of full rank are used. It follows that

[TABLE]

A transient simulation of the dynamical systems often requires implicit integration schemes. The smaller dimensionality of the system (16) causes a lower computation work in the linear algebra algorithms. However, the definition of the reduced nonlinear function (17) implies that a function evaluation of $\bar{F}$ still takes an evaluation of the original function $\hat{F}$. Significant reductions of the computational effort can be achieved only if one replaces the function evaluations of $\hat{F}$ by cheaper approximations, which is also called hyper-reduction in this context. Several techniques have been proposed for this purpose like (discrete) empirical interpolation, missing point estimation or piecewise linear approaches. We refer to [7] and the references therein. A hyper-reduction is not within the scope of this paper. Our main aim is the identification of sufficiently accurate low-dimensional representations.

Furthermore, systems of DAEs often require more sophisticated MOR methods in comparison to systems of ODEs, see [4, 17]. Indeed the situation becomes more critical for higher-index DAEs (index larger than one).

4.2 Proper orthogonal decomposition

The projection-based ROM (16),(17) is determined uniquely by the matrices $T_{\rm l}$ and $T_{\rm r}$. We apply the POD to identify a projection matrix. A numerical integration scheme yields a solution of an initial value problem for the larger system (11) or (15) in $[t_{0},t_{\rm end}]$. In grid points $t_{0}<t_{1}<\cdots<t_{\ell-1}=t_{\rm end}$, we obtain approximations, which are called snapshots in this context. We collect them in a matrix

[TABLE]

Often it holds that $\ell\ll mn$. Now a singular value decomposition (SVD), see [15, p. 76] yields

[TABLE]

with orthogonal matrices $U\in\mathbbm{R}^{mn\times mn}$, $Q\in\mathbbm{R}^{\ell\times\ell}$ and a diagonal matrix $\Sigma\in\mathbbm{R}^{\ell\times\ell}$ with the singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{\ell}\geq 0$. Let $u_{1},\ldots,u_{mn}$ be the columns of $U$. The right-hand projection matrix is defined as

[TABLE]

for any $r\leq\ell$. Thus just the singular vectors associated with the $r$ dominant singular values are included in the MOR scheme. The left-hand projection matrix is chosen by a Galerkin-type approach as $T_{\rm l}=T_{\rm r}$.

The POD technique requires a solution of an initial value problem of the FOM to generate the snapshots. Thus saving computational effort is achieved only if the ROM can be reused for other numerical simulations. These additional numerical simulations may involve

• •

other input signals $u(t)$,

• •

other initial values $\bar{v}(t_{0})$,

• •

longer time intervals $[t_{0},t^{\prime}]$ with $t^{\prime}>t_{\rm end}$.

The third case will be examined for a test example in Section 5.

4.3 Low-dimensional subspaces and approximation

The following derivation is applicable in the case of any projection-based MOR as discussed in Section 4.1. The output of the ROM (16) produces an approximation of the QoI in the original system (1) by

[TABLE]

for $t\in[t_{0},t_{\rm end}]$ and $p\in\Pi$. The $\mathcal{L}^{2}(\Pi,\rho)$-error of this reduction with respect to (9) reads as

[TABLE]

due to the orthonormality of the basis and Parseval’s equality. Since a linear dependence of the QoI $y$ on the state variables or inner variables $x$ is assumed in (1), the same derivation applies as in [32]. In (16), the part $\bar{w}=\bar{C}\bar{v}$ with the output matrix $\bar{C}=(\bar{c}_{ij})$ yields

[TABLE]

The approximation (20) exhibits the formulation

[TABLE]

Thus new basis functions

[TABLE]

are defined. We obtain the low-dimensional representation (10) with $\bar{w}_{j}=\bar{v}_{j}$ for $j=1,\ldots,r$. Hence the coefficients are already identified by the MOR scheme.

Let $\Phi=(\Phi_{1},\ldots,\Phi_{m})^{\top}$ and $\Psi=(\Psi_{1},\ldots,\Psi_{r})^{\top}$. It holds that $\Psi=\bar{C}^{\top}\Phi$. The set of functions (22) is linearly independent in most of the cases. However, they are not orthogonal in general. An orthonormal basis can be constructed from the original basis by an SVD of the output matrix

[TABLE]

including orthogonal matrices $U\in\mathbbm{R}^{m\times m}$, $Q\in\mathbbm{R}^{r\times r}$ and a diagonal matrix $\Sigma\in\mathbbm{R}^{r\times r}$ with the singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{r}\geq 0$. Thus we obtain the alternative basis functions

[TABLE]

with coefficients from the matrix $U$. If follows that $\{\Psi_{1}^{*},\ldots,\Psi_{r}^{*}\}$ is an orthonormal basis spanning the same subspace as $\{\Psi_{1},\ldots,\Psi_{r}\}$. The alternative basis can be deflated to $\{\Psi_{1}^{*},\ldots,\Psi_{k}^{*}\}$ for a $k\leq r$ to remove a (numerical) rank deficiency. The details are explained in [32].

The original orthonormal basis $\{\Phi_{1},\ldots,\Phi_{m}\}$ owns the advantage that the first coefficient yields the expected value due to $\Phi_{1}\equiv 1$. The first and second moment can be obtained as well in the novel orthonormal basis. Expected value and variance read as

[TABLE]

due to (24). Hence just the SVD (23) has to be computed at the beginning. The moments do not require significant additional work. Yet the only advantage of the basis $\{\Psi_{1}^{*},\ldots,\Psi_{r}^{*}\}$ in comparison to $\{\Phi_{1},\ldots,\Phi_{m}\}$ is the reduced dimensionality.

4.4 Best approximation

A projection-based MOR (16),(17) identifies a low-dimensional representation (10) including its time-dependent coefficient functions. In addition, we obtain a best approximation within the subspace spanned by the new basis (22), where just information from the right-hand projection matrix $T_{\rm r}$ is used.

The best approximation

[TABLE]

with $w^{*}=(w_{1}^{*},\ldots,w_{r}^{*})^{\top}\in\mathbbm{R}^{r}$ is defined by the optimisation problem

[TABLE]

with $w=(w_{1},\ldots,w_{r})^{\top}\in\mathbbm{R}^{r}$ pointwise for $t\geq t_{0}$. This best approximation can be computed from the solution of the FOM by a linear least squares problem.

Theorem 1

Let $\hat{w}=(\hat{w}_{1},\ldots,\hat{w}_{m})^{\top}$ be the output of the dynamical system (11) or (15) for an initial value problem. Let the output matrix $\bar{C}=\hat{C}T_{\rm r}$ have full rank. The coefficients $w^{*}=(w_{1}^{*},\ldots,w_{r}^{*})$ of the best approximation (25) minimising (26) are the solution of the linear least squares problem

[TABLE]

pointwise for $t\geq t_{0}$.

Proof:

We omit the dependence on time for notational convenience. The Hilbert space norm can be expressed by inner products

[TABLE]

We obtain

[TABLE]

Using (22), basic calculations yield

[TABLE]

with $\bar{C}=(\bar{c}_{ij})$. It follows that

[TABLE]

The degrees of freedom are $\bar{w}\in\mathbbm{R}^{r}$, whereas $\hat{w}\in\mathbbm{R}^{m}$ is constant. A necessary condition for a minimum is a vanishing gradient. Thus we achieve

[TABLE]

which is equivalent to the normal equation, see [39, p. 232], associated with the linear least squares problem (27). $\Box$

The least squares problem (27) implies the formula

[TABLE]

provided that the output matrix exhibits full rank. Due to (28), the best approximation (25) is continuous in time, because the solution of the dynamical system (11) or (15) is assumed to be continuous. The computation of the best approximation requires to solve the FOM. However, the application of the formula (28) afterwards is cheap, because the transformation matrix is time-invariant. A QR-decomposition of the output matrix can be reused at all time points.

4.5 Error analysis

The difference $\hat{y}-\bar{y}$ between the QoI (9) from a numerical method and the approximation (10) is exactly the error of the MOR approach. Hence this error depends on the individual choice of the projection-based MOR method.

A more detailed examination is feasible for the best approximation using the POD method of Section 4.2. Again the Galerkin system (11) is considered without loss of generality, since the analysis also applies to the collocation system (15). We obtain the following property of the approximation quality with respect to the state variables or the inner variables.

Lemma 1

Let $\hat{v}\in C^{1}$ be the exact solution of an initial value problem of the dynamical system (11). Let $T_{\rm r}\in\mathbbm{R}^{mn\times r}$ be the projection matrix from the POD approach with $\ell>r$ snapshots. For each $t\in[t_{0},t_{\rm end}]$, the solution $v^{*}(t)\in\mathbbm{R}^{r}$ of the linear least squares problem

[TABLE]

satisfies the estimate

[TABLE]

with the time step size $\Delta t=\max\{t_{j}-t_{j-1}\;:\;j=1,\ldots,\ell-1\}$ and the singular value $\sigma_{r+1}$ from the POD.

Proof:

The intermediate value theorem of differential calculus yields componentwise $\hat{v}_{i}(t)=\hat{v}_{i}(t_{j})+(t-t_{j})\dot{\hat{v}}_{i}(\xi)$ for each $t\in[t_{j},t_{j+1}]$ and $i=1,\ldots,mn$ with intermediate values $\xi$. It follows that

[TABLE]

Since $v^{*}$ represents an optimum, we obtain

[TABLE]

for any $\breve{v}\in\mathbbm{R}^{r}$. The first term can be bounded by

[TABLE]

and the estimate (30) for each $t\in[t_{j},t_{j+1}]$. For the second term, we apply the equality $\hat{v}(t_{j})=Ve_{j}$ with the matrix (18) and a canonical unit vector $e_{j}\in\mathbbm{R}^{\ell}$. Let $V_{r}\in\mathbbm{R}^{mn\times\ell}$ be the closest rank-$r$ approximation identified by the SVD, see [15, p. 79]. It holds that

[TABLE]

where the diagonal matrix $\Sigma_{r}\in\mathbbm{R}^{\ell\times\ell}$ contains only the $r$ dominant singular values. Furthermore, we obtain $V_{r}=T_{\rm r}\widetilde{\Sigma}_{r}Q^{\top}$ with a modified diagonal matrix $\widetilde{\Sigma}_{r}$, because $T_{\rm r}$ consists of the first $r$ columns of $U$. Thus we choose the vector $\breve{v}=\widetilde{\Sigma}_{r}Q^{\top}e_{j}$. It follows that

[TABLE]

due to the error estimate for the spectral norm in [15, p. 79]. $\Box$

The solutions of the least squares problems (29) are smooth again, because $C^{1}$-solutions of the FOM are assumed. However, we do not require this smoothness of the optimum in the following. Now we show an error bound on the QoI.

Theorem 2

Let $\hat{v}\in C^{1}$ be the solution and $\hat{y}$ be the QoI of an initial value problem of the dynamical system (11). The best approximation $y^{*}$ with respect to the subspace, which is identified by the POD with $\ell>r$ snapshots, satisfies the error bound

[TABLE]

for each $t\in[t_{0},t_{\rm end}]$ with the time step size $\Delta t=\max\{t_{j}-t_{j-1}\;:\;j=1,\ldots,\ell-1\}$ and the singular value $\sigma_{r+1}$ from the POD.

Proof:

It holds that $\Psi=\bar{C}^{\top}\Phi$ with $\Psi:=(\Psi_{1},\ldots,\Psi_{r})^{\top}$ and $\Phi:=(\Phi_{1},\ldots,\Phi_{m})^{\top}$ due to (22). The best approximation exhibits a representation

[TABLE]

with $w^{*}=(w_{1}^{*},\ldots,w_{r}^{*})^{\top}$. The POD yields the projection matrix $T_{\rm r}$. The orthonormality of the basis functions $\Phi$ allows for the application of Parseval’s equality. For each $t\in[t_{0},t_{\rm end}]$, it follows that

[TABLE]

for any $\breve{w}\in\mathbbm{R}^{r}$. Thus Lemma 1 yields the bound (31) by the choice $\breve{w}=v^{*}$ from the least squares problem (29). $\Box$

The spectral norm of the output matrix is often uncritical in the estimate (31). For example, if the single output of (1) is just a state variable or inner variable, then it follows that $\|\hat{C}\|_{2}=1$. The error bound (31) becomes small in the case of a fast decay of the singular values and a sufficiently small time step size. On the one hand, the singular values are computed a priori, which yields the first term of the estimate. On the other hand, the second term is not computable, because the involved derivatives are unknown.

5 Illustrative examples

Now we apply an MOR approach to the stochastic Galerkin system as well as a stochastic collocation formulation for two test examples. All numerical calculations are performed within the software package MATLAB.

5.1 Scrapie model

We use a model of the scrapie disease from [10, p. 37]. Reaction kinetics yields an autonomous system of ODEs

[TABLE]

for three molecular concentrations $x_{1},x_{2},x_{3}$. The system (32) owns the form (1) with $B=0$ and an identity matrix $E$. As output, we examine the three concentrations separately. The five parameters are reaction constants, where we arrange the nominal values $p^{*}=(10^{-5},0.1,1,10^{-4},0.1)^{\top}$ from [10]. Initial values (2) are chosen as $x(0,p)=(1,0,0.1)^{\top}$ for all $p$. The time interval $[0,500]$ is considered.

In the stochastic modelling, we replace the parameters by independent uniformly distributed random variables with 10% variation around the nominal values $p^{*}$. The PC expansion (5) includes multivariate polynomials, which are the products of univariate Legendre polynomials. Our truncated expansion (7) involves all polynomials up to total degree three, which implies $m=56$ basis functions. The stochastic Galerkin method yields a dynamical system (11) with $mn=168$ state variables. Since just quadratic nonlinearities are given in the right-hand side of the system (32), the right-hand sides of (11) can be evaluated exactly except for roundoff errors. In a stochastic collocation approach, we consider a tensor-product formula of the Gauss-Legendre quadrature, see [39, p. 171], with $k=3^{5}=243$ nodes. The accompanying dynamical system (15) exhibits $kn=729$ state variables.

The trapezoidal rule yields the numerical solutions of all initial value problems in this test example. Variable time step sizes are determined by a local error control with relative tolerance $\varepsilon_{\rm rel}$ and absolute tolerance $\varepsilon_{\rm abs}$. In the stochastic Galerkin method, we apply the choice $\varepsilon_{\rm rel}=10^{-4}$, $\varepsilon_{\rm abs}=10^{-6}$ for the computation of the snapshots, which causes 149 steps and thus 150 snapshots including the initial values. These snapshots are also used to obtain an approximation of the expected values (first coefficient) as well as the standard deviations (other coefficients) of the three random concentrations shown in Figure 1. Now the POD method requires the SVD (19). Figure 2 (left) illustrates the computed singular values, which decay rapidly. The projection matrices $T_{\rm l},T_{\rm r}$ and the ROM (16) follow from the POD technique for user-defined reduced dimensions. Since the dimensionality of the FOM (11) is relatively small, we cannot expect saving computational effort by an MOR. Nevertheless, we discuss the MOR to show the feasibility of the approach for the identification of a low-dimensional representation.

Now both the FOM (11) and the ROMs (16) are integrated with accuracy requirements $\varepsilon_{\rm rel}=10^{-3}$, $\varepsilon_{\rm abs}=10^{-6}$. Initial values are transformed via $\bar{v}(0)=T_{\rm r}^{\top}\hat{v}(0)$. We obtain the approximation (20) from solving the ROM as well as the best approximation (25) for each concentration as QoI. The $\mathcal{L}^{2}(\Pi,\rho)$-error (21) between FOM and ROM can be determined pointwise in time. We evaluate the approximations at 200 equidistant time points in the total time interval using interpolation in time. Figure 3 depicts the $\mathcal{L}^{2}(\Pi,\rho)$-error (21) for different reduced dimensions, where the maximum error of all time points is determined. We recognise that the MOR error decreases rapidly for increasing dimensions at the beginning and then stagnates. The error becomes dominated by the quality of the snapshots, which depends on the density of the associated time points, and the error of the time integration. In contrast, the error of the best approximation decreases further for increasing dimensions and achieves a much lower magnitude.

We repeat the strategy for the stochastic collocation method now. Each system (32) is integrated separately for the nodes of the quadrature with different step sizes for tolerances $\varepsilon_{\rm rel}=10^{-4}$, $\varepsilon_{\rm abs}=10^{-6}$. Interpolation yields 200 snapshots of the weakly coupled system (15) at equidistant time points. The result of the POD is shown in Figure 2 (right), which demonstrates a fast decay of the singular values again. We impose the accuracy requirements $\varepsilon_{\rm rel}=10^{-3}$, $\varepsilon_{\rm abs}=10^{-6}$ on the integration of both the FOM (15) and the ROMs. Initial values are obtained by $\bar{x}(0)=T_{\rm r}^{\top}\hat{x}(0)$. Figure 4 illustrates the maximum $\mathcal{L}^{2}(\Pi,\rho)$-errors (21) for the approximations (20) and (25) of the QoIs in the case of different reduced dimensions. A comparison to Figure 3 shows that the quality of the approximations is slightly worse in this stochastic collocation approach. We remark that the accuracy of the stochastic Galerkin method or the stochastic collocation method with respect to the exact random QoI in (1) is not compared, because this item is out of scope. We determine the error of the MOR separately in the Galerkin technique and the collocation scheme.

5.2 Transistor amplifier circuit

We consider the electric circuit of a transistor amplifier depicted in Figure 5. A mathematical modelling yields a dynamical system (1) of DAEs ($n=5$) with differential index one, which is given in [18, p. 377]. The system is nonlinear, where the mapping $F$ models the bipolar transistor of the circuit and thus involves exponential functions. Three capacitances and six resistances are included in the matrices $E$ and $A$, respectively. The matrix $B\in\mathbbm{R}^{5\times 2}$ inserts the inputs $u=(U_{\rm in},U_{\rm op})^{\top}$. The operating voltage $U_{\rm op}$ is constant. We put the value $U_{\rm op}$ into the matrix $B$ and the input becomes $u=(U_{\rm in},\widetilde{U}_{\rm op})^{\top}$ with $\widetilde{U}_{\rm op}\equiv 1$, because $U_{\rm op}$ is a parameter in the matrix $B$ now. A single time-varying input is supplied by the voltage source $U_{\rm in}$, which we select as the sinusoidal input

[TABLE]

The unknowns of the system consist of the five node voltages. The QoI is defined as the output voltage $U_{\rm out}$ at the fifth node. Thus the matrix $C$ in (1) is just a unit vector. Figure 6 shows the numerical solution of an initial value problem (1),(2) for constant physical parameters from [18].

In this test example, all numerical solutions of initial value problems are computed by the backward differentiation formulas (BDF), see [39, p. 531]. A local error control with tolerances $\varepsilon_{\rm rel}$ and $\varepsilon_{\rm abs}$ yields adaptive time step sizes as well as adaptive orders (1 to 5).

For an uncertainty quantification, we choose all capacitances, all resistances and the operating voltage as random parameters ($q=10$) with independent uniform distributions varying 1% around the constant parameters from above. Concerning the orthogonal expansion (7), we apply all basis polynomials up to total degree three and obtain $m=286$ terms.

The stochastic Galerkin method generates a DAE system (11) with $mn=1430$ inner variables. We require a quadrature formula (13) to evaluate the nonlinear right-hand side approximately. We use a sparse grid quadrature, see [14], which is adapted from the Gauss-Legendre rule with level 3 and $k=2441$ nodes. Initial values for the system (11) are determined as in the previous example. Now the numerical solution with tolerances $\varepsilon_{\rm rel}=10^{-5}$, $\varepsilon_{\rm abs}=10^{-6}$ yields 222 snapshots including the initial values in the time interval $[0,T]$. On the one hand, the snapshots imply an approximation of the coefficient functions in the truncated expansion (7) of the output voltage. Figure 7 illustrates the maximum coefficients occurring in the discrete time points. We observe different orders of magnitudes, which indicates a potential for a sufficiently accurate low-dimensional approximation of the QoI. On the other hand, a POD of the snapshots reveals the singular values within (19) shown by Figure 8 (left). We obtain the associated projection matrices $T_{\rm l},T_{\rm r}$ for user-defined reduced dimensions.

Now we repeat the procedure for a stochastic collocation technique, where we apply the sparse grid quadrature from above. The auxiliary system (15) becomes a DAE with differential index one and $kn=12205$ inner variables. Concerning the time integrations, the same tolerances are used as in the stochastic Galerkin approach. The numerical solution of the initial value problem of (15) is interpolated onto 200 equidistant time points in $[0,T]$, which yields 201 snapshots. Alternatively, if the dynamical systems (1) are solved separately for the nodes of the quadrature (as done in Section 5.1), then the interpolated snapshots cause failures within the transient simulation of some reduced systems in this example. Figure 8 (right) shows the dominating singular values from (19) associated with the snapshots from the collocation system (15).

Both the stochastic Galerkin system (11) and the collocation system (15) are reduced by the POD method now. Initial value problems of FOMs and ROMs are solved with tolerances $\varepsilon_{\rm rel}=10^{-3}$, $\varepsilon_{\rm abs}=10^{-6}$ in the following. We consider the total time interval $[0,3T]$, whereas the snapshots are located in $[0,T]$ only. Sometimes the transient simulation of an ROM fails, because the true dynamics is not captured. The reasons are that the dimension of an ROM is not large enough or the snapshots do not reveal some required information.

A Monte-Carlo simulation generates a reference solution, where initial value problems of the original DAEs (1) are solved for $10^{4}$ samples of the random parameters. The time integrations are done with high accuracy requirements $\varepsilon_{\rm rel}=10^{-6}$, $\varepsilon_{\rm abs}=10^{-8}$. Figure 9 depicts the computed expected value as well as standard deviation. We approximate the expected value as well as the variance of the QoI using the FOMs and their ROMs. These statistics are compared to the reference solution in Table 1. The maximum difference is determined on the time interval $[0,T]$ of the snapshots and the longer time interval $[0,3T]$. In the interval of the snapshots, the differences decay monotone for increasing dimensions of the ROMs. The ROMs do not improve from $r=25$ to $r=50$, because the total error is already dominated by the quality of the snapshots and the error of the time integration. In the longer interval, the monotonicity is not given for increasing dimensions, because the snapshots do not reproduce all required information.

We discuss QoIs from the reduced systems of dimension $r=25$ for both Galerkin approach and collocation technique in the time interval $[0,3T]$ further. We investigate the difference of an approximation (20) from an ROM and a best approximation (25) to a solution (9) of an FOM. Figure 10 depicts the computed $\mathcal{L}^{2}(\Pi,\rho)$-errors (21) depending on time. The behaviour of the error coincides in all four variants. The errors are relatively small in the first cycle, where the snapshots are located. Yet the error increases at later cycles. Of course the error of the best approximation is always smaller than the MOR error. In addition, the magnitudes of the $\mathcal{L}^{2}(\Pi,\rho)$-errors agree for Galerkin method and collocation method. Hence the collocation technique becomes favourable in this example, because the computational effort of a time integration is lower.

We note that computation work is not decreased in the used MOR for this test example, because hyper-reduction is not included and thus the nonlinear functions of the FOMs still have to be evaluated in the ROMs. There is a potential for saving computing time by the usage of hyper-reduction as mentioned in Section 4.1.

6 Conclusions

Orthogonal expansions were applied to the solution of random nonlinear dynamical systems. An MOR for the stochastic Galerkin system or an auxiliary system of a stochastic collocation method implied low-dimensional approximations of the expansions. On the one hand, a transient MOR method yields a low-dimensional representation directly. On the other hand, the projection of a solution of the full-order model onto the subspace, which is identified by the MOR, generates a best approximation. Numerical simulations demonstrated that both the Galerkin method and collocation techniques are feasible to determine adequate low-dimensional approximations. However, a transient simulation of a reduced-order model may become critical or fail in the case of complex nonlinear dynamical systems. In contrast, the best approximation is robust and identifies an accurate low-dimensional approximation, while still the full-order model has to be solved. The strategy of hyper-reduction is required if the proposed methods shall save computing time in the transient simulations.