A new type of singular perturbation approximation for stochastic bilinear systems

TL;DR

This paper introduces a novel singular perturbation approximation method for stochastic bilinear systems, providing the first proven $L^2$-error bound for this class, enhancing model order reduction accuracy.

Contribution

It extends existing MOR techniques to stochastic bilinear systems and establishes a new $L^2$-error bound for SPA, which was previously unproven even for deterministic cases.

Findings

Proposed a modified reduced order model with a different reachability Gramian.

Proved an $L^2$-error bound for SPA in stochastic bilinear systems.

The error bound is novel even for deterministic bilinear systems.

Abstract

Model order reduction (MOR) techniques are often used to reduce the order of spatially-discretized (stochastic) partial differential equations and hence reduce computational complexity. A particular class of MOR techniques is balancing related methods which rely on simultaneously diagonalizing the system Gramians. This has been extensively studied for deterministic linear systems. The balancing procedure has already been extended to bilinear equations [1], an important subclass of nonlinear systems. The choice of Gramians in [1] is referred to be the standard approach. In [18], a balancing related MOR scheme for bilinear systems called singular perturbation approximation (SPA) has been described that relies on the standard choice of Gramians. However, no error bound for this method could be proved. In this paper, we extend the setting used in [18] by considering a stochastic system with bilinear drift and linear diffusion term. Moreover, we propose a modified reduced order model and choose a different reachability Gramian. Based on this new approach, an $L^2$-error bound is proved for SPA which is the main result of this paper. This bound is new even for deterministic bilinear systems.

Equations

$$dx(t)$$

$$y(t)$$

$${\Vert u\Vert }_{L_T^2}^2:={\int }_0^T{\Vert u(t)\Vert }_2^{2}dt<\infty$$

$$\lambda \left(A\otimes I+I\otimes A+\sum _{k=1}^m{N}_k\otimes N_k+\sum _{i,j=1}^v{H}_i\otimes H_j{k}_{ij}\right)\subset {\mathbb{C}}_{-}.$$

$$A^T{P}^{-1}+P^{-1}A+\sum _{k=1}^m{N}_k^T{P}^{-1}{N}_k+\sum _{i,j=1}^v{H}_i^T{P}^{-1}{H}_j{k}_{ij}$$

$$A^{T}Q+QA+\sum _{k=1}^m{N}_k^{T}Q{N}_k+\sum _{i,j=1}^v{H}_i^{T}Q{H}_j{k}_{ij}$$

$$(A,B,C,H_i,N_k)\mapsto (\tilde{A},\tilde{B},\tilde{C},{\tilde{H}}_i,{\tilde{N}}_k):=(SA{S}^{-1},SB,C{S}^{-1},S{H}_i{S}^{-1},S{N}_k{S}^{-1}),$$

$$\tilde{A}=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right],\tilde{B}=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right],{\tilde{N}}_k=\left[\begin{array}{cc}{N}_{k,11}& N_{k,12}\\ N_{k,21}& N_{k,22}\end{array}\right],{\tilde{H}}_i=\left[\begin{array}{cc}{H}_{i,11}& H_{i,12}\\ H_{i,21}& H_{i,22}\end{array}\right],\tilde{C}=\left[\begin{array}{cc}{C}_1& C_2\end{array}\right],$$

$$\displaystyle\tilde{x}=\left[\begin{array}[]{c}x_{1}\\ x_{2}\end{array}\right]\text{ and }\Sigma=\left[\begin{array}[]{cc}\Sigma_{1}&\\ &\Sigma_{2}\end{array}\right],$$

$$d{x}_r$$

$$y_r(t)$$

$$\bar{A}:=A_{11}-A_{12}{A}_{22}^{-1}{A}_{21},\bar{B}:=B_1-A_{12}{A}_{22}^{-1}{B}_2,\bar{C}:=C_1-C_2{A}_{22}^{-1}{A}_{21},$$

$$\bar{D}:=-C_2{A}_{22}^{-1}{B}_2,{\bar{E}}_k:=-N_{k,12}{A}_{22}^{-1}{B}_2,{\bar{F}}_i:=-H_{i,12}{A}_{22}^{-1}{B}_2,$$

$${\bar{H}}_i:=H_{i,11}-H_{i,12}{A}_{22}^{-1}{A}_{21},{\bar{N}}_k:=N_{k,11}-N_{k,12}{A}_{22}^{-1}{A}_{21},$$

$${\left(\mathbb{E}{\Vert y-y_r\Vert }_{L_T^2}^2\right)}^{\frac{1}{2}}\le 2({\tilde{\sigma }}_1+{\tilde{\sigma }}_2+\dots +{\tilde{\sigma }}_{\nu }){\Vert u\Vert }_{L_T^2}\exp \left(0.5{\Vert u^0\Vert }_{L_T^2}^2\right),$$

$$A^T{\Sigma }^{-1}+{\Sigma }^{-1}A+\sum _{k=1}^m{N}_k^T{\Sigma }^{-1}{N}_k+\sum _{i,j=1}^v{H}_i^T{\Sigma }^{-1}{H}_j{k}_{ij}$$

$$A^T\Sigma +\Sigma A+\sum _{k=1}^m{N}_k^T\Sigma N_k+\sum _{i,j=1}^v{H}_i^T\Sigma H_j{k}_{ij}$$

$$y_{\mp }(t)$$

$$d{x}_{-}$$

$$d{x}_{+}=[A{x}_{\pm }+2Bu-\left[\begin{array}{c}0\\ c_0\end{array}\right]+\sum _{k=1}^m{N}_k{x}_{\pm }{u}_k]dt+\sum _{i=1}^v[H_i{x}_{\pm }-\left[\begin{array}{c}0\\ c_i\end{array}\right]]d{M}_i.$$

$${\left(\mathbb{E}{\Vert y-y_r\Vert }_{L_T^2}^2\right)}^{\frac{1}{2}}\le 2\sigma {\Vert u\Vert }_{L_T^2}\exp \left(0.5{\Vert u^0\Vert }_{L_T^2}^2\right).$$

$$\mathbb{E}\left[x_{-}^T(t)\Sigma x_{-}(t)\right]=$$

$$+{\int }_0^t\sum _{i,j=1}^v\mathbb{E}\left[{\left(H_i{x}_{\mp }+\left[\begin{array}{c}0\\ c_i\end{array}\right]\right)}^T\Sigma \left(H_j{x}_{\mp }+\left[\begin{array}{c}0\\ c_j\end{array}\right]\right)\right]k_{ij}ds.$$

$$\sum _{k=1}^{m}2x_{-}^T(s)\Sigma N_k{x}_{\mp }(s)u_k(s)$$

$$\le \sum _{k=1}^m{\Vert {\Sigma }^{\frac{1}{2}}{x}_{-}(s)u_k^0(s)\Vert }_2^2+{\Vert {\Sigma }^{\frac{1}{2}}{N}_k{x}_{\mp }(s)\Vert }_2^2$$

$$=x_{-}^T(s)\Sigma x_{-}(s){\Vert u^0(s)\Vert }_2^2+\sum _{k=1}^m{x}_{\mp }^T(s)N_k^T\Sigma N_k{x}_{\mp }(s),$$

$$2x_{-}^T(s)\Sigma A{x}_{\mp }(s)$$

$$=x_{\mp }^T(s)(A^T\Sigma +\Sigma A)x_{\mp }(s)-2{\left[\begin{array}{c}0\\ h(s)\end{array}\right]}^T\Sigma A{x}_{\mp }(s),$$

$$\mathbb{E}\left[x_{-}^T(t)\Sigma x_{-}(t)\right]\le$$

$$+\mathbb{E}{\int }_0^{t}2x_{-}^T\Sigma \left[\begin{array}{c}0\\ c_0\end{array}\right]+\sum _{i,j=1}^v{\left(2H_i{x}_{\mp }+\left[\begin{array}{c}0\\ c_i\end{array}\right]\right)}^T\Sigma \left[\begin{array}{c}0\\ c_j\end{array}\right]k_{ij}ds$$

$$+{\int }_0^t\mathbb{E}\left[x_{-}^T\Sigma x_{-}\right]{\Vert u^0\Vert }_2^{2}ds-\mathbb{E}{\int }_0^{t}2{\left[\begin{array}{c}0\\ h\end{array}\right]}^T\Sigma A{x}_{\mp }ds.$$

$${\left(2H_i{x}_{\mp }+\left[\begin{array}{c}0\\ c_i\end{array}\right]\right)}^T\Sigma \left[\begin{array}{c}0\\ c_j\end{array}\right]={\left(2H_i{x}_{\mp }+\left[\begin{array}{c}0\\ c_i\end{array}\right]\right)}^T\left[\begin{array}{c}0\\ {\Sigma }_2{c}_j\end{array}\right]$$

$$={\left(2\left[\begin{array}{cc}{H}_{i,21}& H_{i,22}\end{array}\right](x-\left[\begin{array}{c}{x}_r\\ -h\end{array}\right])+c_i\right)}^T{\Sigma }_2{c}_j={\left(2\left[\begin{array}{cc}{H}_{i,21}& H_{i,22}\end{array}\right]x-c_i\right)}^T{\Sigma }_2{c}_j,$$

$$-2{\left[\begin{array}{c}0\\ h\end{array}\right]}^T\Sigma A{x}_{\mp }$$

$$=-2h^T{\Sigma }_2(\left[\begin{array}{cc}{A}_{21}& A_{22}\end{array}\right]x+B_{2}u),$$

$$\mathbb{E}\left[x_{-}^T(t)\Sigma x_{-}(t)\right]\le$$

$$+\mathbb{E}{\int }_0^{t}2x_2^T{\Sigma }_2{c}_0+\sum _{i,j=1}^v{\left(2\left[\begin{array}{cc}{H}_{i,21}& H_{i,22}\end{array}\right]x-c_i\right)}^T{\Sigma }_2{c}_j{k}_{ij}ds$$

$$-\mathbb{E}{\int }_0^{t}2h^T{\Sigma }_2(\left[\begin{array}{cc}{A}_{21}& A_{22}\end{array}\right]x+B_{2}u)ds.$$

Full text

1 Introduction

Many phenomena in real life can be described by partial differential equations (PDEs). For an accurate mathematical modeling of these real world applications, it is often required to take random effects into account. Uncertainties in a PDE model can, for example, be represented by an additional noise term leading to stochastic PDEs (SPDEs) [11, 15, 27, 28].

It is often necessary to numerically approximate time-dependent SPDEs since analytic solutions do not exist in general. Discretizing in space can be considered as a first step. This can, for example, be done by spectral Galerkin [17, 19, 20] or finite element methods [2, 21, 22]. This usually leads to large-scale SDEs. Solving such complex SDE systems causes large computational cost. In this context, model order reduction (MOR) is used to save computational time by replacing high dimensional systems by systems of low order in which the main information of the original system should be captured.

1.1 Literature review

Balancing related MOR schemes were developed for deterministic linear systems first. Famous representatives of this class of methods are balanced truncation (BT) [3, 25, 26] and singular perturbation approximation (SPA) [14, 23].

BT was extended in [5, 8] and SPA was generalized in [32] to stochastic linear systems. With this first extension, however, no $L^{2}$-error bound can be achieved [6, 12]. Therefore, an alternative approach based on a different reachability Gramian was studied for stochastic linear systems leading to an $L^{2}$-error bound for BT [12] and for SPA [31].

BT [1, 5] and SPA [18] were also generalized to bilinear systems, which we refer to as the standard approach for these systems. Although bilinear terms are very weak nonlinearities, they can be seen as a bridge between linear and nonlinear systems. This is because many nonlinear systems can be represented by bilinear systems using a so-called Carleman linearization. Applications of these equations can be found in various fields [10, 24, 33]. The standard approach for bilinear system has the drawback that no $L^{2}$-error bound could be shown so far. A first error bound for the standard ansatz was recently proved in [4], where an output error bound in $L^{\infty}$ was formulated for infinite dimensional bilinear systems. Based on the alternative choice of Gramians in [12], a new type of BT for bilinear systems was considered [30] providing an $L^{2}$-error bound under the assumption of a possibly small bound on the controls.

A more general setting extending both the stochastic linear and the deterministic bilinear case was investigated in [29]. There, BT was studied and an $L^{2}$-error bound was proved overcoming the restriction of bounded controls in [30]. In this paper, we consider SPA for the same setting as in [29] in order to generalize the work in [18]. Moreover, we modify the reduced order model (ROM) in comparison to [18] and show an $L^{2}$-error bound which closes the gap in the theory in this context.

For further extensions of balancing related MOR techniques to other nonlinear systems, we refer to [7, 34].

1.2 Setting and ROM

Let every stochastic process appearing in this paper be defined on a filtered probability space $\left(\Omega,\mathcal{F},\left(\mathcal{F}_{t}\right)_{t\geq 0},\mathbb{P}\right)$111We assume that $\left(\mathcal{F}_{t}\right)_{t\geq 0}$ is right-continuous and $\mathcal{F}_{0}$ contains all sets $A$ with $\mathbb{P}(A)=0$.. Suppose that $M=\left(M_{1},\ldots,M_{v}\right)^{T}$ is an $\left(\mathcal{F}_{t}\right)_{t\geq 0}$-adapted and $\mathbb{R}^{v}$-valued mean zero Lévy process with $\mathbb{E}\left\|M(t)\right\|^{2}_{2}=\mathbb{E}\left[M^{T}(t)M(t)\right]<\infty$ for all $t\geq 0$. Moreover, we assume that for all $t,h\geq 0$ the random variable $M\left(t+h\right)-M\left(t\right)$ is independent of $\mathcal{F}_{t}$.

We consider a large-scale stochastic control system with bilinear drift that can be interpreted as a spatially-discretized SPDE. We investigate the system

[TABLE]

We assume that $A,N_{k},H_{i}\in\mathbb{R}^{n\times n}$ ($k\in\left\{1,\ldots,m\right\}$ and $i\in\left\{1,\ldots,v\right\}$), $B\in\mathbb{R}^{n\times m}$ and $C\in\mathbb{R}^{p\times n}$. Moreover, we define $x(t-):=\lim_{s\uparrow t}x(s)$. The control $u=\left(u_{1},\ldots,u_{m}\right)^{T}$ is assumed to be deterministic and square integrable, i.e.,

[TABLE]

for every $T>0$. By [27, Theorem 4.44] there is a matrix $K=\left(k_{ij}\right)_{i,j=1,\ldots,v}$ such that $\mathbb{E}[M(t)M^{T}(t)]=Kt$. $K$ is called covariance matrix of $M$.

In this paper, we study SPA to obtain a ROM. SPA is a balancing related method and relies on defining a reachability Gramian $P$ and an observability Gramian $Q$. These matrices are selected, such that $P$ characterizes the states in (1a) and $Q$ the states in (1b) which barely contribute to the system dynamics, see [29] for estimates on the reachability and observability energy. The estimates in [29] are global, whereas the standard choice of Gramians leads to results being valid in a small neighborhood of zero only [5, 16].

In order to ensure the existence of these Gramians, throughout the paper it is assumed that

[TABLE]

Here, $\lambda\left(\cdot\right)$ denotes the spectrum of a matrix. The reachability Gramian $P$ and the observability Gramian $Q$ are, according to [29], defined as the solutions to

[TABLE]

where the existence of a positive definite solution to (3) goes back to [12, 31].

We approximate the large scale system (1) by a system which has a much smaller state dimension $r\ll n$. This reduced order model (ROM) is supposed be chosen, such that the corresponding output $y_{r}$ is close to the original one, i.e., $y_{r}\approx y$ in some metric. In order to be able to remove both the unimportant states in (1a) and (1b) simultaneously, the first step of SPA is a state space transformation

[TABLE]

where $S=\Sigma^{-\tfrac{1}{2}}X^{T}L_{Q}^{T}$ and $S^{-1}=L_{P}Y\Sigma^{-\tfrac{1}{2}}$. The ingredients of the balancing transformation are computed by the Cholesky factorizations $P=L_{P}L_{P}^{T}$, $Q=L_{Q}L_{Q}^{T}$, and the singular value decomposition $X\Sigma Y^{T}=L_{Q}^{T}L_{P}$. This transformation does not change the output $y$ of the system, but it guarantees that the new Gramians are diagonal and equal, i.e., $SPS^{T}=S^{-T}QS^{-1}=\Sigma=\mathop{\operator@font diag}\nolimits(\sigma_{1},\ldots,\sigma_{n})$ with $\sigma_{1}\geq\ldots\geq\sigma_{n}$ being the Hankel singular values (HSVs) of the system.

We partition the balanced coefficients of (1) as follows:

[TABLE]

where $A_{11},N_{k,11},H_{i,11}\in\mathbb{R}^{r\times r}$ ($k\in\left\{1,\ldots,m\right\}$ and $i\in\left\{1,\ldots,v\right\}$), $B_{1}\in\mathbb{R}^{r\times m}$ and $C_{1}\in\mathbb{R}^{p\times r}$ etc. Furthermore, we partition the state variable $\tilde{x}$ of the balanced system and the diagonal matrix of HSVs

[TABLE]

where $x_{1}$ takes values in $\mathbb{R}^{r}$ ($x_{2}$ accordingly), $\Sigma_{1}$ is the diagonal matrix of large HSVs and $\Sigma_{2}$ contains the small ones.

Based on the balanced full model (1) with matrices as in (5), the ROM is obtained by neglecting the state variables $x_{2}$ corresponding to the small HSVs. The ROM using SPA is obtained by setting $dx_{2}(t)=0$ and furthermore neglecting the diffusion and bilinear term in the equation related to $x_{2}$. The resulting algebraic constraint can be solved and leads to $x_{2}(t)=-A_{22}^{-1}(A_{21}x_{1}(t)+B_{2}u(t))$. Inserting this expression into the equation for $x_{1}$ and into the output equation, the reduced system is

[TABLE]

with matrices defined by

[TABLE]

where $x_{r}(0)=0$ and the time dependence in (11a) is omitted to shorten the notation. This straight forward ansatz is based on observations from the deterministic case ($N_{k}=H_{i}=0$), where $x_{2}$ represents the fast variables, i.e., $\dot{x}_{2}(t)\approx 0$ after a short time, see [23].

This ansatz for stochastic systems might, however, be false, no matter how small the HSVs corresponding to $x_{2}$ are. Despite the fact that for the motivation, a maybe less convincing argument is used, this leads to a viable MOR method for which an error bound can be proved. An averaging principle would be a mathematically well-founded alternative to this naive approach. Averaging principles for stochastic systems have for example been investigated in [35, 36]. A further strategy to derive a ROM in this context can be found in [9].

Moreover, notice that system (11) is not a bilinear system anymore due to the quadratic term in the control $u$. This is an essential difference to the ROM proposed in [18].

1.3 Main result

The work in this paper on SPA for system (1) can be interpreted as a generalization of the deterministic bilinear case [18]. This extension builds a bridge between stochastic linear systems and stochastic nonlinear systems such that SPA can possibly be applied to many more stochastic equations and applications.

In this paper, we provide an alternative to [29], where BT was studied. We extend the work of [18] combined with a modification of the ROM and the choice of a new Gramian defined through (3). Based on this, we obtain an error bound that was not even available for the deterministic bilinear case. This is the main result of this paper and is formulated in the following theorem. Its proof requires new techniques that cannot be found in the literature so far.

Theorem 1**.**

Let $y$ be the output of the full model (1) with $x(0)=0$ and $y_{r}$ be the output of the ROM (11) with zero initial state. Then, for all $T>0$, it holds that

[TABLE]

where $\tilde{\sigma}_{1},\tilde{\sigma}_{2},\ldots,\tilde{\sigma}_{\nu}$ are the distinct diagonal entries of $\Sigma_{2}=\mathop{\operator@font diag}\nolimits(\sigma_{r+1},\ldots,\sigma_{n})=\mathop{\operator@font diag}\nolimits(\tilde{\sigma}_{1}I,\tilde{\sigma}_{2}I,\ldots,\tilde{\sigma}_{\nu}I)$ and $u^{0}=(u^{0}_{1},\dots,u_{m}^{0})^{T}$ is the control vector with components defined by $u_{k}^{0}\equiv\begin{cases}0&\text{if }N_{k}=0,\\ u_{k}&\text{else}.\end{cases}$

Theorem 1 is proved in Section 2.3. We observe that an exponential term enters the bound in Theorem 1 which is due to the bilinearity in the drift. Setting $N_{k}=0$ for all $k=1,\ldots,m$ the exponential becomes a one which is the bound of the stochastic linear case [31]. The result in Theorem 1 tells us that the ROM (11) yields a very good approximation if the truncated HSVs (diagonal entries of $\Sigma_{2}$) are small and the vector $u^{0}$ of control components with a non-zero $N_{k}$ is not too large. The exponential in the error bound can be an indicator that SPA performs badly if $u^{0}$ is very large.

The remainder of the paper deals with the proof of Theorem 1.

2 $L^{2}$-error bound for SPA

The proof of the error bound in Theorem 1 is divided into two parts. We first investigate the error that we encounter by removing the smallest HSV from the system in Section 2.1. In this reduction step, the structure from the full model (1) to the ROM (11) changes. Therefore, when removing the other HSVs from the system, another case needs to be studied in Section 2.2. There, an error bound between two ROM is achieved which are neighboring, i.e., the larger ROM has exactly one HSV more than the smaller one. The results of Sections 2.1 and 2.2 are then combined in Section 2.3 in order to prove the general error bound.

For simplicity, let us from now on assume that system (1) is already balanced and has a zero initial condition ($x_{0}=0$). Thus, (3) and (4) become

[TABLE]

i.e., $P=Q=\Sigma=\mathop{\operator@font diag}\nolimits(\sigma_{1},\ldots,\sigma_{n})>0$.

2.1 Error bound of removing the smallest HSV

We introduce the variable $x_{\mp}=\left[\begin{smallmatrix}{x}_{1}-x_{r}\\ x_{2}+A_{22}^{-1}(A_{21}x_{r}+B_{2}u)\end{smallmatrix}\right]$ since the corresponding output

[TABLE]

is the output error between the full and the reduced system (11). We aim to find an equation for $x_{\mp}$. This is done through the state variable $x_{-}=\left[\begin{smallmatrix}{x}_{1}-x_{r}\\ x_{2}\end{smallmatrix}\right]$. The differential $d(x_{1}-x_{r})$ is obtained by subtracting the state equation (11a) of the reduced system from the first $r$ rows of (1a). The corresponding right side is then rewritten using $x_{\mp}$. Moreover, the right side of the differential of $x_{2}$, compare with the last $n-r$ rows of (1a), is also formulated with the help of $x_{\mp}$. This results in

[TABLE]

where $c_{0}(t):=\sum_{k=1}^{m}[N_{k,21}x_{r}(t)-N_{k,22}A_{22}^{-1}(A_{21}x_{r}(t)+B_{2}u(t))]u_{k}(t)$ and $c_{i}(t):=H_{i,21}x_{r}(t)-H_{i,22}A_{22}^{-1}(A_{21}x_{r}(t)+B_{2}u(t))$ for $i=1,\ldots,v$.

We furthermore introduce the reverse state to $x_{\mp}$ in terms of the signs. This is $x_{\pm}=\left[\begin{smallmatrix}{x}_{1}+x_{r}\\ x_{2}-A_{22}^{-1}(A_{21}x_{r}+B_{2}u)\end{smallmatrix}\right]$. Using the state $x_{+}=\left[\begin{smallmatrix}{x}_{1}+x_{r}\\ x_{2}\end{smallmatrix}\right]$, with a differential obtained by combining (1a) and (11a) again, and expressing its right side with $x_{\pm}$, we have

[TABLE]

We will see that the proof of the error bound can be reduced to the task of finding suitable estimates for $\mathbb{E}[x_{-}^{T}(t)\Sigma x_{-}(t)]$ and $\mathbb{E}[x_{+}^{T}(t)\Sigma^{-1}x_{+}(t)]$. This idea was also used to determine an error bound for BT [29]. However, the proof for SPA requires different techniques to find the estimates.

Theorem 2**.**

Let $y$ be the output of the full model (1) with $x(0)=0$, $y_{r}$ be the output of the ROM (11) with $x_{r}(0)=0$ and $\Sigma_{2}=\sigma I$, $\sigma>0$, in (10). Then, it holds that

[TABLE]

Proof.

We derive a suitable upper bound for $\mathbb{E}[x_{-}^{T}(t)\Sigma x_{-}(t)]$ first applying Ito’s formula. Hence, Lemma 4 and Equation (15) yield

[TABLE]

We find an estimate for the terms related to $N_{k}$, that is

[TABLE]

where $u^{0}$ is defined as in Theorem 1. Moreover, adding a zero, we rewrite

[TABLE]

where $h(s)=A_{22}^{-1}(A_{21}x_{r}(s)+B_{2}u(s))$ With (18) and (19), (17) becomes

[TABLE]

Taking the partitions of $x_{-}$ and $\Sigma$ into account, we see that $x_{-}^{T}\Sigma\left[\begin{smallmatrix}{0}\\ c_{0}\end{smallmatrix}\right]=x_{2}^{T}\Sigma_{2}c_{0}$. Furthermore, the partitions of $x_{\mp}$ and $H_{i}$ yield

[TABLE]

since $\left[\begin{smallmatrix}{H}_{i,21}&H_{i,22}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{r}\\ -h\end{smallmatrix}\right]=c_{i}$. Using the partition of $A$, it holds that

[TABLE]

because $\left[\begin{smallmatrix}{A}_{21}&A_{22}\end{smallmatrix}\right]\left[\begin{smallmatrix}{-}x_{r}\\ h\end{smallmatrix}\right]=B_{2}u$. We insert (13) and (14) into inequality (20) and exploit the relations in (21) and (22). Hence,

[TABLE]

We define the function $\alpha_{-}(t):=\mathbb{E}\int_{0}^{t}2x_{2}^{T}\Sigma_{2}c_{0}+\sum_{i,j=1}^{v}\left(2\left[\begin{smallmatrix}{H}_{i,21}&H_{i,22}\end{smallmatrix}\right]x-c_{i}\right)^{T}\Sigma_{2}c_{j}k_{ij}ds-\mathbb{E}\int_{0}^{t}2h^{T}\Sigma_{2}(\left[\begin{smallmatrix}{A}_{21}&A_{22}\end{smallmatrix}\right]x+B_{2}u)ds$ and apply Lemma 6 implying

[TABLE]

Since $\Sigma$ is positive definite, we obtain and upper bound for the output error by

[TABLE]

Defining the term $\alpha_{+}(t):=\mathbb{E}\int_{0}^{t}2x_{2}^{T}\Sigma_{2}^{-1}c_{0}+\sum_{i,j=1}^{v}\left(2\left[\begin{smallmatrix}{H}_{i,21}&H_{i,22}\end{smallmatrix}\right]x-c_{i}\right)^{T}\Sigma_{2}^{-1}c_{j}k_{ij}ds-\mathbb{E}\int_{0}^{t}2h^{T}\Sigma_{2}^{-1}(\left[\begin{smallmatrix}{A}_{21}&A_{22}\end{smallmatrix}\right]x+B_{2}u)ds$ and exploiting the assumption that $\Sigma_{2}=\sigma I$, leads to

[TABLE]

The remaining step is to find a bound for the right side of (23) that does not depend on $\alpha_{+}$ anymore. For that reason, a bound for the expression $\mathbb{E}[x_{+}^{T}(t)\Sigma^{-1}x_{+}(t)]$ is derived next using Ito’s lemma again. From (16) and Lemma 4, we obtain

[TABLE]

Analogously to (18), it holds that

[TABLE]

Additionally, we rearrange the term related to $A$ as follows

[TABLE]

Moreover, we have

[TABLE]

We plug in the above results into (24) which gives us

[TABLE]

From inequality (12) and the Schur complement condition on definiteness, it follows that

[TABLE]

We multiply (28) with $\left[\begin{smallmatrix}{x}_{\pm}\\ 2u\end{smallmatrix}\right]^{T}$ from the left and with $\left[\begin{smallmatrix}{x}_{\pm}\\ 2u\end{smallmatrix}\right]$ from the right. Hence,

[TABLE]

Applying this result to (25) yields

[TABLE]

We first of all see that $x_{+}^{T}\Sigma^{-1}\left[\begin{smallmatrix}{0}\\ c_{0}\end{smallmatrix}\right]=x_{2}^{T}\Sigma_{2}^{-1}c_{0}$ using the partitions of $x_{+}$ and $\Sigma$. With the partition of $H_{i}$, we moreover have

[TABLE]

In addition, it holds that

[TABLE]

Plugging the above relations into (30) leads to

[TABLE]

We add $2\mathbb{E}\int_{0}^{t}\sum_{i,j=1}^{v}c_{i}^{T}\Sigma_{2}^{-1}c_{j}k_{ij}ds$ to the right side of (31) and preserve the inequality since this term is a nonnegative due to Lemma 5. This results in

[TABLE]

Gronwall’s inequality in Lemma 6 yields

[TABLE]

We find an estimate for the following expression:

[TABLE]

Combining (32) with (33), we obtain

[TABLE]

Comparing this result with (23) implies

[TABLE]

∎

We proceed with the study of an error bound between two ROM that are neighboring.

2.2 Error bound for neighboring ROMs

In this section, we investigate the output error between two ROMs, in which the larger ROM has exactly one HSV than the smaller one. This concept of neighboring ROMs was first introduced in [31] but in the much simpler stochastic linear setting.

The reader might wonder why a second case is considered besides the one in Section 2.1 since one might just start with a full model that has the same structure as the ROM (11). The reason is that is not clear how the Gramians need to be chosen for (11). In order to investigate the error between two ROMs by SPA, a finer partition than the one in (5) is required. We partition the matrices of the balanced full system (1) as follows:

[TABLE]

The partitioned balanced solution to (1a) and the Gramians are then of the form

[TABLE]

We introduce the ROM of truncating $\Sigma_{3}$ first. According to the procedure described in Section 1.2, the reduced system is obtained by setting $dx_{3}$ equal to zero, neglecting the bilinear and the diffusion term in this equation. The solution $\tilde{x}_{3}$ of the resulting algebraic constraint is an approximation for $x_{3}$. One can solve for this approximating variable and obtains $\tilde{x}_{3}=-A_{33}^{-1}(A_{31}x_{1}+A_{32}x_{2}+B_{3}u)$. Inserting this result for $x_{3}$ in the equations for $x_{1}$, $x_{2}$ and into the output equation (1b) leads to

[TABLE]

where $\left[\begin{smallmatrix}{x}_{1}(0)\\ x_{2}(0)\end{smallmatrix}\right]=\left[\begin{smallmatrix}{0}\\ 0\end{smallmatrix}\right]$ and

[TABLE]

We aim to determine the error between this ROM and the reduced system of neglecting $\Sigma_{2}$ and $\Sigma_{3}$. This is

[TABLE]

where $x_{r}(0)=0$,

[TABLE]

and we define

[TABLE]

In order to find a bound for the error between (38b) and (39b), state variables analogously to $x_{\mp}$ and $x_{\pm}$ in Section 2.1 are constructed in the following and corresponding equations are derived. For simplicity, we use a similar notation again and define

[TABLE]

One can see that these states are obtained by combining the states appearing on the right sides of (38a) and (39a). Furthermore, the output of $\hat{x}_{\mp}$ leads to the output error

[TABLE]

which is a direct consequence of (38b) and (39b).

Now, we find the differential equations for $\hat{x}_{\mp}$ and $\hat{x}_{\mp}$. Using (40), we find that

[TABLE]

Applying the first line of (42), we obtain the following equation

[TABLE]

where $\hat{c}_{0}=\sum_{k=1}^{m}\left[\begin{smallmatrix}{N}_{k,21}&{N}_{k,22}&{N}_{k,23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{r}\\ -h_{1}\\ -h_{2}\end{smallmatrix}\right]u_{k}$ and $\hat{c}_{i}=\left[\begin{smallmatrix}{H}_{i,21}&{H}_{i,22}&{H}_{i,23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{r}\\ -h_{1}\\ -h_{2}\end{smallmatrix}\right]$ for $i=1,\ldots,v$. We supplement (39a) with (43) and combine this with (38a). Hence, we obtain

[TABLE]

where $\hat{x}_{-}=\left[\begin{smallmatrix}{x}_{1}-x_{r}\\ x_{2}\end{smallmatrix}\right]$ and furthermore

[TABLE]

where $\hat{x}_{+}=\left[\begin{smallmatrix}{x}_{1}+x_{r}\\ x_{2}\end{smallmatrix}\right]$. We now state the output error between the systems (38) and (39) for the case that the ROM are neighboring, i.e., the larger model has exactly one HSV more than the smaller one.

Theorem 3**.**

Let $\bar{y}$ be the output of the ROM (38), $\bar{y}_{r}$ be the output of the ROM (39) and $\Sigma_{2}=\sigma I$, $\sigma>0$, in (37). Then, it holds that

[TABLE]

Proof.

We make use of equations (44) and (45) in order to prove this bound. We set $\hat{\Sigma}=\left[\begin{smallmatrix}{\Sigma}_{1}&\\ &\Sigma_{2}\end{smallmatrix}\right]$ as a submatrix of $\Sigma$ in (37). Lemma 4 now yields

[TABLE]

We see that the right side of (46) contains the submatrices $\hat{A},\hat{B},\hat{H},\hat{N}$ and $\hat{\Sigma}$. In order to be able to refer to the full matrix inequality (13), we find upper bounds for certain terms in the following involving the full matrices $A,B,H,N$ and $\Sigma$. With the same estimate as in (18) and the control vector $u^{0}$ defined in Theorem 1, we have

[TABLE]

Adding the term $\sum_{k=1}^{m}\left(\left[\begin{smallmatrix}{N}_{k,31}&{N}_{k,32}&{N}_{k,33}\end{smallmatrix}\right]\hat{x}_{\mp}(s)\right)^{T}\Sigma_{3}\left[\begin{smallmatrix}{N}_{k,31}&{N}_{k,32}&{N}_{k,33}\end{smallmatrix}\right]\hat{x}_{\mp}(s)$ to the right side of this inequality results in

[TABLE]

Moreover, it holds that

[TABLE]

We derive $\left[\begin{smallmatrix}{A}_{31}&{A}_{32}&A_{33}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]=-B_{3}u$ by the definition of $\tilde{x}_{3}$. Moreover, it can be seen from the second line of (42) that $\left[\begin{smallmatrix}{A}_{31}&{A}_{32}&A_{33}\end{smallmatrix}\right]\hat{x}_{\mp}=0$. Hence,

[TABLE]

It remains to find a suitable upper bound related to the expression depending on $\hat{H}_{i}$. We first of all see that

[TABLE]

The term $\sum_{i,j=1}^{v}\left(\left[\begin{smallmatrix}{H}_{i,31}&{H}_{i,32}&{H}_{i,33}\end{smallmatrix}\right]\hat{x}_{\mp}(s)\right)^{T}\Sigma_{3}\left[\begin{smallmatrix}{H}_{j,31}&{H}_{j,32}&{H}_{j,33}\end{smallmatrix}\right]\hat{x}_{\mp}(s)k_{ij}$ is nonnegative through Lemma 5. Adding this term to the right side of the above equation yields

[TABLE]

Applying (47), (48) and (49) to (46), results in

[TABLE]

Using that $\hat{c}_{i}=\left[\begin{smallmatrix}{H}_{i,21}&{H}_{i,22}&{H}_{i,23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{r}\\ -h_{1}\\ -h_{2}\end{smallmatrix}\right]$, we have

[TABLE]

It can be seen further that

[TABLE]

taking the first line of (42) into account. Inserting (51) and (52) into (50) and using the fact that $2\hat{x}_{-}^{T}\hat{\Sigma}\left[\begin{smallmatrix}{0}\\ \hat{c}_{0}\end{smallmatrix}\right]=2x_{2}\Sigma_{2}\hat{c}_{0}$ leads to

[TABLE]

where we set $\hat{\alpha}_{-}(t):=\mathbb{E}\int_{0}^{t}2x_{2}^{T}\Sigma_{2}\hat{c}_{0}+\left(2\left[\begin{smallmatrix}{H}_{i,21}&H_{i,22}&H_{i,23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]-\hat{c}_{i}\right)^{T}\Sigma_{2}\hat{c}_{j}ds-\mathbb{E}\int_{0}^{t}2h_{1}^{T}\Sigma_{2}(\left[\begin{smallmatrix}{A}_{21}&A_{22}&A_{23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]+B_{2}u)ds$. With (13) and (41), we obtain

[TABLE]

Applying Lemma 6 to this inequality yields

[TABLE]

Since the above left side of the inequality is positive, we obtain

[TABLE]

We exploit that $\Sigma_{2}=\sigma I$. Hence, we have

[TABLE]

where we set $\hat{\alpha}_{+}(t):=\mathbb{E}\int_{0}^{t}2x_{2}^{T}\Sigma_{2}^{-1}\hat{c}_{0}+\left(2\left[\begin{smallmatrix}{H}_{i,21}&H_{i,22}&H_{i,23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]-\hat{c}_{i}\right)^{T}\Sigma_{2}^{-1}\hat{c}_{j}ds-\mathbb{E}\int_{0}^{t}2h_{1}^{T}\Sigma_{2}^{-1}(\left[\begin{smallmatrix}{A}_{21}&A_{22}&A_{23}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]+B_{2}u)ds$. In order to find a suitable bound for the right side of (54), Ito’s lemma is applied to $\mathbb{E}[\hat{x}_{+}^{T}(t)\hat{\Sigma}^{-1}\hat{x}_{+}(t)]$. Due to (45) and Lemma 4, we obtain

[TABLE]

Analogously to (47), it holds that

[TABLE]

Furthermore, we see that

[TABLE]

Since $\left[\begin{smallmatrix}{A}_{31}&{A}_{32}&A_{33}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{1}\\ x_{2}\\ \tilde{x}_{3}\end{smallmatrix}\right]=\left[\begin{smallmatrix}{A}_{31}&{A}_{32}&A_{33}\end{smallmatrix}\right]\left[\begin{smallmatrix}{x}_{r}\\ -h_{1}\\ -h_{2}\end{smallmatrix}\right]=-B_{3}u$ by the definition of $\tilde{x}_{3}$ and the second line of (42), we obtain $\left[\begin{smallmatrix}{A}_{31}&{A}_{32}&A_{33}\end{smallmatrix}\right]\hat{x}_{\pm}=-2B_{3}u$. Thus,

[TABLE]

Finally, we see that

[TABLE]

applying Lemma 5. With (56), (57) and (58) inequality (55) becomes

[TABLE]

Similar to (29), we obtain

[TABLE]

This leads to

[TABLE]

In the following (60) is expressed by terms depending on $\Sigma_{2}$. We obtain $\hat{x}_{+}^{T}\hat{\Sigma}^{-1}\left[\begin{smallmatrix}{0}\\ \hat{c}_{0}\end{smallmatrix}\right]=x_{2}^{T}\Sigma_{2}^{-1}\hat{c}_{0}$ exploiting the partitions of $\hat{x}_{+}$ and $\hat{\Sigma}$. The terms depending on $\hat{H}_{i}$ become

[TABLE]

adding $\sum_{i,j=1}^{v}\hat{c}_{i}^{T}\Sigma_{2}^{-1}\hat{c}_{j}k_{ij}$ which is positive due to Lemma 5. Furthermore, using the first line of (42), it holds that

[TABLE]

We insert (61) and (62) into (60) and obtain

[TABLE]

With Lemma 6, analogously to (34), we find

[TABLE]

The relations (54) and (63) yield the claim. ∎

2.3 Proof of Theorem 1

We apply the results in Theorems 2 and 3. We remove the HSVs step by step and exploit the triangular inequality in order to bound the error between the outputs $y$ and $y_{r}$. We have

[TABLE]

where $y_{r_{i}}$ are the outputs of the ROMs with dimensions $r_{i}$ defined by $r_{i+1}=r_{i}+m(\tilde{\sigma}_{i})$ for $i=1,2\ldots,\nu-1$. Here, $m(\tilde{\sigma}_{i})$ denotes the multiplicity of $\tilde{\sigma}_{i}$ and $r_{1}=r$. In the reduction step from $y$ to $y_{r_{\nu}}$ only the smallest HSV $\tilde{\sigma}_{\nu}$ is removed from the system. Hence, by Theorem 2, we have

[TABLE]

The ROMs of the outputs $y_{r_{j}}$ and $y_{r_{j-1}}$ are neighboring according to Section 2.2, i.e., only the HSV $\tilde{\sigma}_{r_{j-1}}$ is removed in the reduction step. By Theorem 3, we obtain

[TABLE]

for $j=2,\ldots,\nu$. This provides the claimed result.

3 Conclusions

In this paper, we investigated a large-scale stochastic bilinear system. In order to reduce the state space dimension, a model order reduction technique called singular perturbation approximation was extended to this setting. This method is based on Gramians proposed in [29] that characterize how much a state contributes to the system dynamics. This choice of Gramians as well as the structure of the reduced system is different than in [18]. With this modification, we provided a new $L^{2}$-error bound that can be used to point out the cases in which the reduced order model by singular perturbation approximation delivers a good approximation to the original model. This error bound is new even for deterministic bilinear systems.

Appendix A Supporting Lemmas

In this appendix, we state three important results and the corresponding references that we frequently use throughout this paper.

Lemma 4**.**

Let $a,b_{1},\ldots,b_{v}$ be $\mathbb{R}^{d}$-valued processes, where $a$ is $\left(\mathcal{F}_{t}\right)_{t\geq 0}$-adapted and almost surely Lebesgue integrable and the functions $b_{i}$ are integrable with respect to the mean zero square integrable Lévy process $M=(M_{1},\ldots,M_{v})^{T}$ with covariance matrix $K=\left(k_{ij}\right)_{i,j=1,\ldots,v}$. If the process $x$ is given by

[TABLE]

then, we have

[TABLE]

Proof.

We refer to [31, Lemma 5.2] for a proof of this lemma. ∎

Lemma 5**.**

Let $A_{1},\ldots,A_{v}$ be $d_{1}\times d_{2}$ matrices and $K=(k_{ij})_{i,j=1,\ldots,v}$ be a positive semidefinite matrix, then

[TABLE]

is also positive semidefinite.

Proof.

The proof can be found in [31, Proposition 5.3]. ∎

Lemma 6** (Gronwall lemma).**

Let $T>0$, $z,\alpha:[0,T]\rightarrow\mathbb{R}$ be measurable bounded functions and $\beta:[0,T]\rightarrow\mathbb{R}$ be a nonnegative integrable function. If

[TABLE]

then it holds that

[TABLE]

for all $t\in[0,T]$.

Proof.

The result is shown as in [13, Proposition 2.1]. ∎