Computing the stochastic $H^\infty$-norm

TL;DR

This paper introduces a method to compute the stochastic $H^inity$-norm of linear systems using a Riccati-type matrix equation, providing a new approach without frequency domain or Hamiltonian conditions.

Contribution

It develops a novel computational technique for the stochastic $H^inity$-norm based on a parametrized algebraic Riccati equation, extending deterministic methods.

Findings

Provides a Riccati-based algorithm for stochastic $H^inity$-norm computation.

Establishes theoretical characterization without frequency domain or Hamiltonian conditions.

Enables practical computation for stochastic control systems.

Abstract

The stochastic $H^\infty$-norm is defined as the $L^2$-induced norm of the input-output operator of a stochastic linear system. Like the deterministic $H^\infty$-norm it is characterised by a version of the bounded real lemma, but without a frequency domain description or a Hamiltonian condition. Therefore, we base its computation on a parametrised algebraic Riccati-type matrix equation.

Tables

Table 1. Table 1: Averaged computing times (in sec) for random systems.

$n$	10	20	40	80	160
LMI	0.11s	0.99s	35.72s	2030s	-
Alg1	4.43s	7.98s	24.43s	156.6s	1156s

Table 2. Table 2: Computing times (in sec) and results for heat equation

$n$	25	36	49	64	81	100
LMI	3.83s	36.92s	306.3s	1810s	8631s	-
Alg1	9.56s	13.25s	26.29s	73.38s	129.1s	177.3s
$\Vert {𝕃}_n\Vert$	0.4724	0.4694	0.4669	0.4647	0.4628	0.4611

Equations

$$dx$$

$$\Vert v(\cdot ){\Vert }_{L_w^2}^2:=\mathcal{E}\left({\int }_0^{\infty }\Vert v(t){\Vert }^{2}dt\right)<\infty ,$$

$$\mathcal{E}(\Vert x(t,x_0,0){\Vert }^2)\stackrel{t\to \infty }{⟶}0,$$

$$\Vert \mathbb{L}\Vert$$

$${\mathcal{R}}_{\gamma }(X)=A^{T}X+XA+N^{T}XN-C^{T}C$$

$$-(B^{T}X-D^{T}C{)}^T({\gamma }^{2}I-D^{T}D{)}^{-1}(B^{T}X-D^{T}C).$$

$$({\mathcal{R}}_{\gamma }{)}_X^{\prime }(\Delta )$$

$$\displaystyle\left[\begin{array}[]{cc}(\mathcal{L}_{A}+\Pi_{N})(X)-C^{T}C&XB-C^{T}D\\ B^{T}X-D^{T}C&\gamma^{2}I-D^{T}D\end{array}\right]>0.$$

$$\Vert \mathbb{L}\Vert$$

$$AP+P{A}^T+NP{N}^T=-B{B}^T.$$

$$\displaystyle\left[\begin{array}[]{cc}(\mathcal{L}_{A}+\Pi_{N})(X)-C^{T}C&XB-C^{T}D\\ B^{T}X-D^{T}C&\gamma^{2}I-D^{T}D\end{array}\right]\geq 0.$$

$$0\le \mathrm{trace}(-B^{T}XB)\le m^2{\gamma }^2\Vert B^T{P}^{†}B{\Vert }_2.$$

$${\mathcal{R}}_{\gamma }^{\det }(X)$$

$${\mathcal{R}}_{{\gamma }_1}^{\det }(X)$$

$$X_{k+1}$$

$$\sigma (({\mathcal{R}}_{\gamma }{)}_{X_k}^{\prime })\subset {\mathbb{C}}_{-},{\mathcal{R}}_{\gamma }(X_k)\le 0,\text{ and }{X}_k\ge X_{k+1}.$$

$$P_0=I,P_{k+1}=-{\mathcal{L}}_{A_{X_k}}^{-1}{\Pi }_N(P_k),{\rho }_k=\frac{\mathrm{trace}(P_k{P}_{k+1})}{\mathrm{trace}(P_k{P}_k)}$$

$$A_{X_k}^T\Delta +\Delta A_{X_k}+N^T\Delta N$$

$$\Delta$$

$${\Delta }_{j+1}$$

$$dx$$

$$y$$

$${\mathcal{R}}_{\gamma }(X)$$

$$P(X)$$

$$S(X)$$

$$Q_{\gamma }(X)$$

$${\mathcal{L}}_A:X\mapsto A^{T}X+XA\text{ and }{\Pi }_N:X\mapsto N^{T}XN.$$

$${\mathcal{L}}_A^{\ast }:X\mapsto AX+X{A}^T\text{ and }{\Pi }_N^{\ast }:X\mapsto NX{N}^T.$$

$$\displaystyle\left[\begin{array}[]{cc}A^{T}X+XA+N^{T}XN-C_{\varepsilon}^{T}C_{\varepsilon}&XB-C_{\varepsilon}^{T}D_{\varepsilon}\\ B^{T}X-D_{\varepsilon}^{T}C_{\varepsilon}&\gamma^{2}I-D_{\varepsilon}^{T}D_{\varepsilon}\end{array}\right]>0\;.$$

$$0$$

$$\displaystyle\stackrel{{\scriptstyle j\to\infty}}{{\to}}\left[\begin{array}[]{cc}A^{T}\tilde{X}+\tilde{X}A+N^{T}\tilde{X}N&\tilde{X}B\\ B^{T}\tilde{X}&0\end{array}\right]\geq 0\;,$$

$$0$$

Full text

1 Introduction

The $H^{\infty}$-norm is a fundamental concept for asymptotically stable deterministic linear time invariant systems. It is equal to the input/output norm of a system both in the frequency and the time domain. It is used in robustness analysis and serves as a performance index in $H^{\infty}$ control. In model order reduction, it is an important measure for the quality of the approximation. There are very efficient algorithms for the computation of the $H^{\infty}$-norm, which are based on a Hamiltonian characterization. The most widely used among these was described in [1, 2], but recent progress has been made e.g. in [3, 4, 5, 6].

A stochastic version of the $H^{\infty}$-norm was introduced by Hinrichsen and Pritchard in [7]. It has a similar range of applications as its deterministic $H^{\infty}$ counterpart, but its numerical computation has hardly been considered in the literature. A major obstacle in transferring ideas and algorithms from the deterministic case is the lack of a suitable frequency domain interpretation or a Hamiltonian characterization in the stochastic setup.

In this note we present an algorithm to compute the stochastic $H^{\infty}$-norm, based on a Riccati characterization. According to the stochastic bounded real lemma, [7], the norm is given as the infimum of all $\gamma>0$ for which a given parametrized Riccati equation has a stabilizing solution. We check the solvability of the Riccati equation by a Newton iteration.

The paper is structured as follows. In Section 2 we introduce stochastic systems, define the stochastic $H^{\infty}$-norm and state the stochastic bounded real lemma. We also provide a new version of the non-strict bounded real lemma and give some new bounds for the stabilizing solution, which are proven in appendix A.2. In Section 3 we describe our basic algorithm and discuss ways to make all the steps fast. In Section 4 we report on numerical experiments. In particular, we compare our algorithm with an LMI solver. To keep the notational burden low, we confine ourselves to the case, where only one multiplicative noise term affects the state vector. Our results can easily be extended to more general situations which we hint at in appendix A.1.

2 The stochastic $H^{\infty}$-norm

We consider stochastic linear systems of the form

[TABLE]

where $A,N\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{p\times n}$, $D\in\mathbb{R}^{p\times m}$, and $w=(w(t))_{t\in\mathbb{R}_{+}}$ is a zero mean real Wiener process on a probability space $(\Omega,{\cal F},\mu)$ with respect to an increasing family $({\cal F}_{t})_{t\in\mathbb{R}_{+}}$ of $\sigma$-algebras ${\cal F}_{t}\subset{\cal F}$ (e.g. [8, 9]).

Let $L^{2}_{w}(\mathbb{R}_{+},\mathbb{R}^{q})$ denote the corresponding space of non-anticipating stochastic processes $v$ with values in $\mathbb{R}^{q}$ and norm

[TABLE]

where ${\cal E}$ denotes expectation. For initial data $x(0)=x_{0}$ and input $u\in L^{2}_{w}(\mathbb{R}_{+},\mathbb{R}^{m})$ we denote the solution and the output of (1) by $x(t,x_{0},u)$ and $y(t,x_{0},u)$, respectively.

Definition 2.1

System (1) is called asymptotically mean-square-stable, if

[TABLE]

for all initial conditions $x_{0}$. In this case, for simplicity, we also call the pair $(A,N)$ asymptotically mean-square stable.

If $(A,N)$ is asymptotically mean-square stable, then (1) defines an input-output operator $\mathbb{L}:u\mapsto y$ from $L^{2}_{w}(\mathbb{R}_{+},\mathbb{R}^{m})$ to $L^{2}_{w}(\mathbb{R}_{+},\mathbb{R}^{p})$ via $u\mapsto y(\cdot,0,u)$, see [7]. By $\|\mathbb{L}\|$ we denote the induced operator norm,

[TABLE]

which is an analogue of the deterministic $H^{\infty}$-norm. We therefore call it the stochastic $H^{\infty}$-norm of system (1).

2.1 The stochastic bounded real lemma

The norm (2) can be characterized by the stochastic bounded real lemma. To this end, we define the quadratic (Riccati-type) mapping $\mathcal{R}_{\gamma}:\mathbb{R}^{n\times n}\to\mathbb{R}^{n\times n}$, which depends on the parameter $\gamma>\|D\|_{2}$, by

[TABLE]

Its Fréchet derivative at some $X\in\mathbb{R}^{n\times n}$ is the linear mapping $(\mathcal{R}_{\gamma})^{\prime}_{X}:\mathbb{R}^{n\times n}\to\mathbb{R}^{n\times n}$ given by

[TABLE]

where $A_{X}=A-B(\gamma^{2}I-D^{T}D)^{-1}(B^{T}X-D^{T}C)$.

Writing $\mathcal{L}_{A}:X\mapsto A^{T}X+XA$ and $\Pi_{N}:X\mapsto N^{T}XN$, we have

[TABLE]

The pair $(A,N)$ is asymptotically mean-square stable if and only if \sigma(\mathcal{L}_{A}+\Pi_{N})\subset\mathbb{C}_{-}=\{\lambda\in\mathbb{C}\;\big{|}\;\operatorname{Re}\lambda<0\}, e.g. [10].

Theorem 2.2

[7]** Assume that $(A,N)$ is asymptotically mean-square stable. For $\gamma>\|D\|_{2}$, the following are equivalent.

• (i)

$\|\mathbb{L}\|<\gamma$.

• (ii)

There exists a negative definite solution $X<0$ to the linear matrix inequality

[TABLE]

• (iii)

There exists a negative definite solution $X<0$ to the strict Riccati inequality $\mathcal{R}_{\gamma}(X)>0$.

• (iv)

There exists a solution $X\leq 0$ to the Riccati equation $\mathcal{R}_{\gamma}(X)=0$, such that $\sigma((\mathcal{R}_{\gamma})^{\prime}_{X})\subset\mathbb{C}_{-}$.

Remark 2.3

A solution of the Riccati equation $\mathcal{R}_{\gamma}(X)=0$, with $\sigma((\mathcal{R}_{\gamma})^{\prime}_{X})\subset\mathbb{C}_{-}$ is called a stabilizing solution. If it exists, then it is uniquely defined, and is the largest solution of the inequality $\mathcal{R}_{\gamma}(X)\geq 0$, see [10]. We will write $X_{+}(\gamma)$ for this solution. By Theorem 2.2, the norm $\|L\|$ is the infimum of all $\gamma$ such that $\mathcal{R}_{\gamma}(X)=0$ possesses a stabilizing solution, i.e.

[TABLE]

Under a controllability assumption we can also give a nonstrict version of Theorem 2.2 for asymptotically mean-square stable systems. We define the controllability Gramian $P$ of system (1) as the solution of

[TABLE]

If the system is stable, then $P$ is nonnegative definite, $P\geq 0$.

Corollary 2.4

Assume that $(A,N)$ is asymptotically mean-square stable and $P>0$ in (8). For $\gamma>\|D\|_{2}$, the following are equivalent.

• (i)

$\|\mathbb{L}\|\leq\gamma$.

• (ii)

There exists a solution $X\leq 0$ to the linear matrix inequality

[TABLE]

• (iii)

There exists a solution $X\leq 0$ to the Riccati equation $\mathcal{R}_{\gamma}(X)=0$.

Moreover, if $\|\mathbb{L}\|=\gamma$, then $\mathcal{R}_{\gamma}(X)=0$ has a largest solution $X=X_{+}(\gamma)$, for which $0\in\sigma\left((\mathcal{R}_{\gamma})^{\prime}_{X}\right)\subset\mathbb{C}_{-}\cup i\mathbb{R}$.

This result is slightly stronger than [11, Proposition 9.6] or [10, Corollary 5.3.14], where it was shown that (i) implies (iii) if $(A,B)$ is controllable. In the appendix we give a new simplified proof, which can also be modified to obtain lower bounds for solutions of (11) as follows.

2.2 Inequalities for solutions of the Riccati equation

Lemma 2.5

Assume that $(A,N)$ is asymptotically mean-square stable, and $\gamma>\|D\|_{2}$. Let $P^{\dagger}\geq 0$ be the Moore-Penrose inverse of $P$ given by (8). If $X\leq 0$ satisfies (11), then

[TABLE]

Note that $\operatorname{trace}(-B^{T}XB)$ is monotonically decreasing. Hence, if (12) is violated for some $X\leq 0$ and $\tilde{X}\leq X$, then $\tilde{X}$ cannot be a solution of (11). This bound is particularly easy to check.

Alternatively, we may compare with solutions of Riccati equations from deterministic control. Let $\mathcal{R}^{\det}_{\gamma}$ denote the counterpart of $\mathcal{R}_{\gamma}$ with $N=0$, i.e.

[TABLE]

Lemma 2.6

*Assume that $(A,N)$ is asymptotically mean-square stable, and $\gamma_{1}\geq\gamma>\|\mathbb{L}\|$.

Then the Riccati equation from the deterministic case*

[TABLE]

possesses a smallest solution $X_{-}\leq 0$, and $X_{-}\leq X$ for all solutions $X$ of (11).

3 Computation of the stochastic $H^{\infty}$-norm

To exploit the characterization (7), we need a method to check, whether the Riccati equation $\mathcal{R}_{\gamma}(X)=0$ possesses a stabilizing solution. Given the Fréchet derivative of $R_{\gamma}(X)$ displayed in (3), it is natural to apply Newton’s method to solve the stochastic algebraic Riccati equation from part (iv) of Theorem 2.2. The following result was proven in [11].

Theorem 3.1

Let $(A,N)$ be mean-square stable and assume that $\gamma>\|\mathbb{L}\|$. Consider the Newton iteration

[TABLE]

where we assume $\sigma((\mathcal{R}_{\gamma})_{X_{0}}^{\prime})\subset\mathbb{C}_{-}$. Then the sequence $X_{k}$ converges to $X_{+}$, and for all $k\geq 1$ it holds that

[TABLE]

Moreover, under the given assumptions $X_{0}=0$ is a suitable initial guess (see appendix for a proof).

Lemma 3.2

Let $(A,N)$ be mean-square stable and assume that $\gamma>\|\mathbb{L}\|$. Then $\sigma((\mathcal{R}_{\gamma})_{0}^{\prime})\subset\mathbb{C}_{-}$.

For a given $\gamma>\|D\|_{2}$, we can check whether $\gamma>\|\mathbb{L}\|$ by running the Newton iteration (14) starting from $X_{0}=0$. If all iterates are stabilizing, and the sequence converges with a given level of tolerance, then we conclude that $\gamma\geq\|\mathbb{L}\|$.

Conversely, if $\gamma<\|\mathbb{L}\|$, then either $\sigma((\mathcal{R}_{\gamma})_{X_{k}}^{\prime})\not\subset\mathbb{C}_{-}$ for some $k$, or the sequence $X_{k}$ is monotonically decreasing and unbounded.

If for some $k$ the condition $\sigma((\mathcal{R}_{\gamma})_{X_{k}}^{\prime})\subset\mathbb{C}_{-}$ is violated or the iteration takes more than a fixed number of steps, then we conclude that $\gamma\leq\|\mathbb{L}\|$. Additionally we might test the conditions of Lemma 2.5 or Lemma 2.6 in each step and conclude that $\gamma\leq\|\mathbb{L}\|$ if one of them is not fulfilled. However, in all our examples only the stability condition was relevant.

Using bisection, we can thus compute $\|\mathbb{L}\|$ up to a given precision.

3.1 The basic algorithm

We summarize this approach as our basic algorithm.

The stability test in line 5 and the solution of the linear system in line 6 are central issues. Both concern the generalized Lyapunov mapping $\mathcal{R}_{X_{k}}^{\prime}$. A naive implementation with general purpose eigenvalue and linear system solvers, respectively, would result in an overall complexity of about $\mathcal{O}(n^{6})$. About the same complexity is required for LMI-solvers. It is, however, well known that standard Lyapunov equations of the form $\mathcal{L}_{A_{X_{k}}}(X)=Y$ can be solved in $\mathcal{O}(n^{3})$ operations, using e.g. the Bartels-Stewart algorithm, [12]. Exploiting this in iterative approaches, we can bring down the complexity of Algorithm 1 at least to $\mathcal{O}(n^{3})$. This will be explained briefly in the following two subsections. Moreover, we suggest a way to choose $\gamma_{0}$ and $\gamma_{1}$ in line 1.

In the numerical experiments, we will show that our algorithm outperforms general purpose LMI methods.

3.2 The stability test

The condition $\sigma((\mathcal{R}_{\gamma})_{X_{k}}^{\prime})\subset\mathbb{C}_{-}$ in line 5 holds if and only if $\sigma(A_{X_{k}})\subset\mathbb{C}_{-}$ and $\rho(\mathcal{L}_{A_{X_{k}}}^{-1}\Pi_{N})<1$, where $\rho$ denotes the spectral radius, [10, Theorem 3.6.1]. Hence, we can first check, whether $\sigma(A_{X_{k}})\subset\mathbb{C}_{-}$ and then apply the power method to compute the spectral radius $\rho$ of $\mathcal{L}_{A_{X_{k}}}^{-1}\Pi_{N}$. Note that the mapping $-\mathcal{L}_{A_{X_{k}}}^{-1}\Pi_{N}$ is nonnegative, in the sense that it maps the cone of nonnegative definite matrices to itself, see [10]. Hence, the iterative scheme

[TABLE]

produces a sequence of nonnegative definite matrices $P_{k}$ which generically converge to the dominant eigenvector. In the limit we have $P_{k+1}\approx\rho P_{k}$, i.e. $\rho_{k}\stackrel{{\scriptstyle k\to\infty}}{{\to}}\rho$.

3.3 The generalized Lyapunov equation

In the Newton step in line 6, the generalized Lyapunov equation

[TABLE]

has to be solved for $\Delta$ to obtain $X_{k+1}=X_{k}+\Delta$. Equations of this type have been studied e.g. in [13].

Note that $\Delta=\Delta^{T}\in\mathbb{R}^{n\times n}$ satisfies the fixed point equation

[TABLE]

The condition $\sigma((\mathcal{R}_{\gamma})^{\prime}_{X_{k}})\subset\mathbb{C}_{-}$ implies $\rho(\mathcal{L}_{A_{X_{k}}}^{-1}\Pi_{N})<1$, where $\rho$ denotes the spectral radius. Hence the fixed point iteration

[TABLE]

is convergent. In each step this iteration only requires the solution of a standard Lyapunov equation. at a cost at most in $\mathcal{O}(n^{3})$. The speed of convergence can be improved by using a Krylov subspace approach like gmres or bicgstab. For details see [13]. More recently, also low-rank techniques have been considered in [14, 15, 16].

3.4 Choosing $\gamma_{0}$ and $\gamma_{1}$

For the bisection it is useful to find suitable upper and lower bounds for $\|\mathbb{L}\|$. Let $G(s)=C(sI-A)^{-1}B+D$ be the transfer function of the deterministic system obtained from (1) by replacing $N$ with zero. The $H^{\infty}$-norm $\|G\|_{H^{\infty}}$ equals the input-output norm of this deterministic system. Then from Theorem 2.2 we conclude $\|G\|_{H^{\infty}}\leq\|\mathbb{L}\|$, because the inequality (6) for a given matrix $X<0$ implies that the corresponding linear matrix inequality with $N=0$ holds for the same $X$. Hence, if $\gamma>\|\mathbb{L}\|$, then $\gamma>\|G\|_{H^{\infty}}$. Therefore, we choose $\gamma_{0}=\|G\|_{H^{\infty}}$ and try $\gamma_{1}=2\gamma_{0}$. If the Newton iteration does not converge for $\gamma_{1}$, then we replace $\gamma_{0}$ by $2\gamma_{0}$ and repeat the previous step, until we have $\gamma_{1}>\|\mathbb{L}\|$.

4 Numerical experiments

The following experiments were carried out on a 2011 MacBook Air with a 1.4 GHz Intel Core 2 Duo processor and 4 GB Memory running OS X 10.11.6 using MATLAB® version R2016b.

4.1 Random systems

We first consider random data $(A,N,B,C)$ produced by randn. The matrix $A$ is made stable by mirroring the unstable eigenvalues at $i\mathbb{R}$. Then the spectral radius $\rho$ of $\mathcal{L}_{A}^{-1}\Pi_{N}$ is estimated as described in subsection 3.2 and an update of $N$ is obtained by multiplication with $(2\rho+1)^{-1/2}$. Thus $(A,N)$ is guaranteed to be mean-square stable. We compute the stochastic $H^{\infty}$-norm by our algorithm and compare it with the result obtained by the MATLAB®-function mincx, see appendix A.3. In all our tests, the relative difference of the computed norms lies within the chosen tolerance level. The computing times, however, differ significantly, see Table 1. While for small dimensions $n$ the implementation of the LMI-solver seems to be superior to our implementation, for larger $n$ the algorithmic complexity becomes relevant. For $n>100$ the LMI-solver is impractical.

4.2 A heat transfer problem

This stochastic modification of a heat transfer problem described in [17] was also discussed in [18]. On the unit square $\Omega=[0,1]^{2}$, the heat equation $T_{t}=\Delta T$ for $T=T(t,x)$ is given with Dirichlet condition $T=u_{j}$, $j=1,2,3$, on three of the boundary edges and a stochastic Robin condition $n\cdot\nabla T=(1/2+\dot{w})T$ on the fourth edge (where $\dot{w}$ stands for white noise). We measure the average value $y(t)=\int_{\Omega}T(t,x)\,dx$

A standard 5-point finite difference discretization on a $k\times k$ grid leads to a modified Poisson matrix $A\in\mathbb{R}^{n\times n}$ with $n=k^{2}$ and corresponding matrices $N\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times 3}$ $C=\tfrac{1}{n}[1,\ldots,1]\in\mathbb{R}^{n\times n}$. The $H^{\infty}$-norm of this discretization $\|\mathbb{L}_{n}\|$ approximates the induced input/output norm of the partial differential equation. Table 2 shows the computing times (in seconds) and the computed norms (which coincide) for the two methods. For $k>9$ the LMI-solver took to long to be considered.

Again, we observe that Algorithm 1 allows to treat larger dimensions than the LMI-solver. However, the computing times for our algorithm also grow fairly fast. As an alternative to bisection one might consider extrapolating the spectral radii $\rho(\gamma)=\rho\left((\mathcal{R}_{\gamma})^{\prime}_{X_{+}(\gamma)}\right)$ which are computed in the course of the process for $\gamma>\|\mathbb{L}\|$, or perhaps the spectral abscissae $\alpha(\gamma)=\max\operatorname{Re}\sigma\left((\mathcal{R}_{\gamma})^{\prime}_{X_{+}(\gamma)}\right)$. Then the norm $\|\mathbb{L}\|$ is given as the value of $\gamma$, where $\rho(\gamma)=1$, or $\alpha(\gamma)=0$. Unfortunately, the slopes of $\rho$ and $\alpha$ are very steep as $\gamma$ approaches $\|\mathbb{L}\|$. Thus an extrapolation does not seem promising. The behaviour is visualized for the heat equation system with $n=25$ in Figure 1.

5 Conclusions

We have suggested an algorithm to compute the stochastic $H^{\infty}$-norm. It builds upon several ideas developed in the literature, and is the first algorithm, whose complexity is considerably smaller than that of a general purpose LMI-solver. We chose to present the algorithm for the simplest case of just one multiplicative noise term, which however can easily be generalized to the class described in appendix A.1. Already in the simple case, the stochastic $H^{\infty}$-norm is much harder to compute than the $H^{\infty}$-norm of a deterministic system, and the computing times are still very high. We see it as a challenge to come up with a faster method.

The note also contains some extensions of known results with new proofs, like the nonstrict stochastic bounded real lemma and lower bounds for Riccati solutions.

Appendix A Appendix

A.1 Generalization

System (1) can be generalized in a straight-forward manner to the case of multiple noise terms at the state and the input (see e.g. [10]). Then our system takes the form

[TABLE]

where $N_{x,j}\in\mathbb{R}^{n\times n}$, $N_{u,j}\in\mathbb{R}^{n\times m}$ and the $w_{j}$ are independent Wiener processes. The Riccati operator $\mathcal{R}_{\gamma}$ then takes the form

[TABLE]

[TABLE]

Our basic algorithm and all our considerations carry over to this case literally. Only the expressions for $\mathcal{R}_{\gamma}$ and $(\mathcal{R}_{\gamma})^{\prime}_{X}$ become more technical.

A.2 Proofs

As above, we write

[TABLE]

On the space of symmetric matrices we consider the scalar product $\langle X,Y\rangle=\operatorname{trace}XY$, and note that the corresponding adjoint operators are

[TABLE]

Further facts on Riccati- and Lyapunov-type operators are cited from [10].

Proof of Corollary 2.4: (iii)$\Rightarrow$(ii) follows from the definiteness criterion via the Schur-complement.

(ii)$\Rightarrow$(iii): If (ii) holds, then $\mathcal{R}_{\gamma}(X)\geq 0$, and by [11] there exists a solution $X_{+}\leq 0$ to the equation $\mathcal{R}_{\gamma}(X)=0$.

(ii)$\Rightarrow$(i): If (6) holds and we replace $C$ and $D$ by $C_{\varepsilon}=\left[\begin{smallmatrix}C\\ \varepsilon I\end{smallmatrix}\right]$ and $D_{\varepsilon}=\left[\begin{smallmatrix}D\\ 0\end{smallmatrix}\right]$, then we get

[TABLE]

This implies $\|\mathbb{L}_{\varepsilon}\|<\gamma$ for the corresponding modified input-output operator. By $\|\mathbb{L}_{\varepsilon}\|\to\|\mathbb{L}\|$ as $\varepsilon\to 0$, we obtain $\|\mathbb{L}\|\leq\gamma$.

(i)$\Rightarrow$(ii): If (i) holds, then $\|\mathbb{L}\|<\gamma+\frac{1}{k}$ for all $k\in\mathbb{N}$, $k>0$. Hence there exist stabilizing solutions $X_{k}\leq 0$ of $\mathcal{R}_{\gamma+\frac{1}{k}}(X)=0$. Moreover, $X_{k}$ is the largest solution of (11) with $\gamma$ replaced by $\gamma+\frac{1}{k}$. Hence it follows that $X_{k+1}\leq X_{k}$ for all $k$. If the $X_{k}$ are bounded below, then the sequence $(X_{k})$ converges and the limit satisfies the nonstrict linear matrix inequality in (iii). Thus it suffices to show boundedness. We assume that the sequence is not bounded, i.e. $\|X_{k}\|\to\infty$ for $k\to\infty$. Consider the normalized sequence $\tilde{X}_{k}=\frac{X_{k}}{\|X_{k}\|}$, which – by Bolzano-Weierstrass – has a convergent subsequence $\tilde{X}_{k_{j}}$ with limit $\tilde{X}\neq 0$. Then

[TABLE]

implying $B^{T}\tilde{X}=0$ and $0\neq A^{T}\tilde{X}+\tilde{X}A+N^{T}\tilde{X}N\geq 0$. Since, by assumption $P>0$, we obtain

[TABLE]

which is a contradiction.

Thus, $\mathcal{R}_{\gamma}(X)=0$ has a solution $X_{\infty}$, which is the limit of the largest and stabilizing solutions $X_{k}$ of $\mathcal{R}_{\gamma+\frac{1}{k}}(X)\geq 0$. Thus $X_{\infty}$ is the largest solution of $\mathcal{R}_{\gamma}(X)=0$ and $\sigma(\mathcal{R}_{\gamma})^{\prime}_{X_{\infty}}\subset\mathbb{C}_{-}\cup i\mathbb{R}$. If $\gamma=\|\mathbb{L}\|$ then $\sigma(\mathcal{R}_{\gamma})^{\prime}_{X_{\infty}}\cap i\mathbb{R}\neq\emptyset$ and [10, Theorem 3.2.3] yields that $0\in\sigma(\mathcal{R}_{\gamma})^{\prime}_{X_{\infty}}$.     $\Box$

Proof of Lemma 2.5:

The controllability Gramian is given by

[TABLE]

In the following consider an arbitrary matrix $X\leq 0$, $X\neq 0$, satisfying $(\mathcal{L}_{A}+\Pi_{N})(X)=Y\geq 0$. Then

[TABLE]

There exists a vector $u\in\mathbb{R}^{m}$ with $\|u\|_{2}=1$ and $u^{*}B^{T}XBu=-\|B^{T}XB\|_{2}\;.$ Moreover

[TABLE]

for $\alpha_{*}=\|B^{T}P^{\dagger}B\|_{2}$. To see this, note that the image of $B$ is contained in the image of $P$. Hence there exists a unitary $U$, such that

[TABLE]

The largest zero of $\chi(\alpha)=\det(\alpha P_{1}-B_{1}B_{1}^{T})$ is

[TABLE]

For $\alpha\geq\alpha_{*}$, we have $\alpha P-BB^{T}\geq 0$ which proves (19).

We set $\mu(X)=\langle Y,P\rangle=|\operatorname{trace}(B^{T}XB)|$. Let now $X$ satisfy (11). With the given data and $\eta>0$ this implies

[TABLE]

If we assume $\mu(X)>m^{2}\gamma^{2}\alpha_{*}$, then the right hand is negative for $\eta=\frac{\mu(X)}{m\gamma^{2}}$, which is a contradiction.

Hence, we have $|\operatorname{trace}(B^{T}XB)|\leq m^{2}\gamma^{2}\|B^{T}P^{\dagger}B\|_{2}$.     $\Box$

      Proof of Lemma 2.6: Note that $\mathcal{R}^{\det}_{\gamma_{1}}(X)\geq\mathcal{R}_{\gamma}(X)$ if $X\leq 0$ and $\gamma\leq\gamma_{1}$ and thus every solution of $\mathcal{R}_{\gamma_{1}}(X)>0$ also satisfies $\mathcal{R}^{\det}_{\gamma_{1}}(X)>0$. Hence $\mathcal{R}^{\det}_{\gamma_{1}}(X)=0$ possesses a stabilizing solution, and, consequently, also an anti-stabilizing solution $X_{-}$, which is the smallest solution of $\mathcal{R}^{\det}_{\gamma_{1}}(X)\geq 0$. Thus also $X_{-}\leq X$ for every solution $X$ of $\mathcal{R}_{\gamma_{1}}(X)\geq 0$.     $\Box$

      Proof of Lemma 3.2: We exploit the concavity of $\mathcal{R}_{\gamma}$ and the resolvent positivity of $(\mathcal{R}_{\gamma})^{\prime}_{0}$, see [10]. If $\|\mathbb{L}\|\leq\gamma$, then there exists $X\leq 0$ such that, by concavity,

[TABLE]

Assume that $\sigma\left((\mathcal{R}_{\gamma})^{\prime}_{0}\right)\not\subset\mathbb{C}_{-}$. Then by [10, Theorem 3.2.3] there exists $H\geq 0$, $\lambda\geq 0$, such that $(\mathcal{R}_{\gamma})^{\prime}_{0}(H)=\lambda H$. Taking the scalar product of inequality (20) with $H$, we get

[TABLE]

It follows that $\mathcal{R}_{\gamma}(0)H=0$, which implies $D^{T}CH=0$ and thus $A_{0}H=AH$. But then $(\mathcal{L}_{A}+\Pi_{N})^{*}(H)=\lambda H$ in contradiction to the stability of $(A,N)$.     $\Box$

A.3 Usage of LMI-solver

The LMI-solver was used as in the following MATLAB® listing.

⬇

1lmis([])

2X = lmivar(1,[n,1]);g = lmivar(1,[1,1]);

3lmiterm([1 1 1 X],N',N);

4lmiterm([1 1 1 X],A',1,'s');

5lmiterm([1 1 1 0],C'*C);

6lmiterm([1 1 2 X],1,B,'s');

7lmiterm([1 2 2 g],-1,1);

8lmisys = getlmis;

9c = mat2dec(lmisys,zeros(n),1);

10options = [tol,0,0,0,1];

11copt = mincx(lmisys,c,options);

12gamma = sqrt(copt)