Numerical solution of Lyapunov equations related to Markov jump linear systems

TL;DR

This paper compares various numerical methods for solving Lyapunov equations in Markov jump linear systems, highlighting the effectiveness of fixed-point and trust-region approaches for large-scale problems.

Contribution

It introduces and evaluates multiple numerical algorithms, including fixed point, Krylov subspace, and optimization-based methods, for solving Lyapunov equations in Markov jump systems.

Findings

Fixed-point approach outperforms optimization methods.

Trust-region method is most effective for large-scale problems.

Application to networked control systems with wireless-induced Markov jumps.

Abstract

We suggest and compare different methods for the numerical solution of Lyapunov like equations with application to control of Markovian jump linear systems. First, we consider fixed point iterations and associated Krylov subspace formulations. Second, we reformulate the equation as an optimization problem and consider steepest descent, conjugate gradient, and a trust-region method. Numerical experiments illustrate that for large-scale problems the trust-region method is more effective than the steepest descent and the conjugate gradient methods. The fixed-point approach, however, is superior to the optimization methods. As an application we consider a networked control system, where the Markov jumps are induced by the wireless communication protocol.

Tables

Table 1. Table 1: Constructed example with n = 4 𝑛 4 n=4 , N = 2 𝑁 2 N=2 and known solution

	Time (seconds)	Iteration	Residual	Error
Algorithm 1	0.0041	3	1.1e$-$11	4.0e$-$13
Algorithm 2	2.5	1116	2.9e$-$05	7.5e$-$06
Algorithm 3	0.93	359	1.9e$-$07	2.5e$-$08
Algorithm 4	0.19	14	5.5e$-$08	2.6e$-$09

Table 2. Table 2: Random examples with fixed Γ ∈ ℝ 2 × 2 Γ superscript ℝ 2 2 \Gamma\in\mathbb{R}^{2\times 2} and positive definite A 1 , A 2 subscript 𝐴 1 subscript 𝐴 2 A_{1},A_{2}

	$n$	Time (seconds)	Iteration	Residual
Algorithm 1	100	0.11	5	7.1e$-$10
Algorithm 2	100	113	30000	1.0e02
Algorithm 3	100	131	30000	5.9e00
Algorithm 4	100	33	30	6.9e$-$08
Algorithm 1	300	1.3	4	1.1e$-$09
Algorithm 4	300	745	35	6.4e$-$08
Algorithm 1	1200	88.0	3.5	5.5e$-$08

Table 3. Table 3: Application of Algorithm 1 to random examples with nonsymmetric A j subscript 𝐴 𝑗 A_{j} and varying Γ Γ \Gamma

$n$	$N$	Time (seconds)	Iteration	Residual
200	20	11.7	6.5	1.1e$-$07
100	100	11.3	6	2.3e$-$08
200	100	64.0	7	1.5e$-$07
1000	3	241	6	3.7e$-$06

Table 4. Table 4: Transition matrix Θ Θ \Theta for τ = 3 𝜏 3 \tau=3 , ν = 2 𝜈 2 \nu=2 with stations A and B .

			e	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1
			t+1	A	A	A	A	B	B	B	B	A	A	A	A	B	B	B	B
			t	A	A	B	B	A	A	B	B	A	A	B	B	A	A	B	B
			t-1	A	B	A	B	A	B	A	B	A	B	A	B	A	B	A	B
e	t	t-1	t-2
0	A	A	A	0.19	0	0	0	0.78	0	0	0	0.01	0	0	0	0.02	0	0	0
0	A	A	B	0.32	0	0	0	0.65	0	0	0	0.01	0	0	0	0.02	0	0	0
0	A	B	A	0	0.42	0	0	0	0.55	0	0	0	0.01	0	0	0	0.02	0	0
0	A	B	B	0	0.42	0	0	0	0.55	0	0	0	0.01	0	0	0	0.02	0	0
0	B	A	A	0	0	0.55	0	0	0	0.42	0	0	0	0.02	0	0	0	0.01	0
0	B	A	B	0	0	0.55	0	0	0	0.42	0	0	0	0.02	0	0	0	0.01	0
0	B	B	A	0	0	0	0.65	0	0	0	0.32	0	0	0	0.02	0	0	0	0.01
0	B	B	B	0	0	0	0.78	0	0	0	0.19	0	0	0	0.02	0	0	0	0.01
1	A	A	A	0.25	0	0	0	0	0	0	0	0.75	0	0	0	0	0	0	0
1	A	A	B	0.25	0	0	0	0	0	0	0	0.75	0	0	0	0	0	0	0
1	A	B	A	0	0.25	0	0	0	0	0	0	0	0.75	0	0	0	0	0	0
1	A	B	B	0	0.25	0	0	0	0	0	0	0	0.75	0	0	0	0	0	0
1	B	A	A	0	0	0	0	0	0	0.25	0	0	0	0	0	0	0	0.75	0
1	B	A	B	0	0	0	0	0	0	0.25	0	0	0	0	0	0	0	0.75	0
1	B	B	A	0	0	0	0	0	0	0	0.25	0	0	0	0	0	0	0	0.75
1	B	B	B	0	0	0	0	0	0	0	0.25	0	0	0	0	0	0	0	0.75

Table 5. Table 5: Application of Algorithm 1 to networked control system

$\nu$	$n$	$N$	Time (seconds)	Iteration	Residual
3	6	54	0.04	3.5	2.3e$-$11
5	10	250	1.1	13.5	2.2e$-$12
10	20	2000	27.3	14	1.5e$-$10

Equations

$$\{\begin{array}{l}\dot{x}=A(r(t))x+B(r(t))u,\\ y=C(r(t))x,\end{array}$$

$$\mathrm{Pr}\left(r(t+\delta )=j|r(t)=i\right)=\{\begin{array}{l}{\pi }_{ij}\delta +o(\delta )\mathrm{if}i\ne j,\\ 1+{\pi }_{ii}\delta +o(\delta )\mathrm{if}i=j,\end{array}$$

$$\underset{t\to \infty }{\lim }E(\Vert x(t){\Vert }_2^2)=0$$

$$\Vert G{\Vert }_{H^2}:=\sqrt{\left(\sum _{s=1}^m\sum _{i=1}^N{\mu }_i\Vert y_{s,i}{\Vert }_2^2\right)},$$

$$\Vert G{\Vert }_{H^2}^2=\sum _{i=1}^N\mathrm{tr}(C_i{P}_i{C}_i^T)=\sum _{i=1}^N{\mu }_i\mathrm{tr}(B_i^T{Q}_i{B}_i),$$

$$A_i{P}_i+P_i{A}_i^T+\sum _{j=1}^N{\pi }_{ji}{P}_j+{\mu }_i{B}_i{B}_i^T=0,$$

$$A_i^T{Q}_i+Q_i{A}_i+\sum _{j=1}^N{\pi }_{ij}{Q}_j+C_i^T{C}_i=0$$

$$0$$

$$\mathcal{H}=\underset{N\text{-times}}{\underbrace{{\mathcal{H}}^n\times \cdots \times {\mathcal{H}}^n}}\text{ and }{\mathcal{H}}_{+}={\mathcal{H}}_{+}^n\times \cdots \times {\mathcal{H}}_{+}^n,$$

$${\mathcal{L}}_A:{\mathcal{H}}^n\to {\mathcal{H}}^n\text{ by }{\mathcal{L}}_A(H)=AH+H{A}^T.$$

$$\mathcal{L}(X_1,\dots ,X_N)=(Z_1,\dots ,Z_N)\text{ with }{Z}_i={\mathcal{L}}_{A_i}(X_i)+{\gamma }_{ii}{X}_i={\mathcal{L}}_{A_i+\frac{{\gamma }_{ii}}{2}I}(X_i)$$

$$\Pi (X_1,\dots ,X_N)=(Z_1,\dots ,Z_N)\text{ with }{Z}_i=\sum _{i\ne j}{\gamma }_{ij}{X}_j.$$

$$(\mathcal{L}+\Pi )(X)=-Y,$$

$$X^{(k+1)}$$

$$X_1^{(k+1)}$$

$$⋮$$

$$X_N^{(k+1)}$$

$$X_1^{(k+1)}$$

$$X_2^{(k+1)}$$

$$⋮$$

$$X_N^{(k+1)}$$

$$(I-T_J)(X)=-{\mathcal{L}}^{-1}(Y)=\tilde{Y}.$$

$$\begin{array}{rl}{\tilde{X}}_1& =-{\mathcal{L}}_{A_1+\frac{{\gamma }_{11}}{2}I}^{-1}\left({\gamma }_{12}{X}_2+{\gamma }_{13}{X}_3+\dots +{\gamma }_{1N}{X}_N\right)\\ {\tilde{X}}_2& =-{\mathcal{L}}_{A_2+\frac{{\gamma }_{22}}{2}I}^{-1}\left({\gamma }_{21}{\tilde{X}}_1+{\gamma }_{23}{X}_3+\dots +{\gamma }_{2N}{X}_N\right)\\ & ⋮\\ {\tilde{X}}_N& =-{\mathcal{L}}_{A_N+\frac{{\gamma }_{NN}}{2}I}^{-1}\left({\gamma }_{N1}{\tilde{X}}_1+{\gamma }_{N2}{\tilde{X}}_2+\dots +{\gamma }_{N,N-1}{\tilde{X}}_{N-1}\right).\end{array}$$

$$\begin{array}{rl}{\tilde{Y}}_1& =-{\mathcal{L}}_{A_1-\frac{{\gamma }_{11}}{2}I}^{-1}(Y_1)\\ {\tilde{Y}}_2& =-{\mathcal{L}}_{A_2-\frac{{\gamma }_{22}}{2}I}^{-1}(Y_2+{\gamma }_{21}{\tilde{Y}}_1)\\ & ⋮\\ {\tilde{Y}}_N& =-{\mathcal{L}}_{A_N-\frac{{\gamma }_{NN}}{2}I}^{-1}(Y_N+{\gamma }_{N1}{\tilde{Y}}_1+\dots +{\gamma }_{N,N-1}{\tilde{Y}}_{N-1}).\end{array}$$

$$\displaystyle\operatorname{vec}\Big{(}\sum_{j=1}^{N}\gamma_{ij}X_{j}\Big{)}$$

$$\mathrm{minimize}f(X):=\sum _{i=1}^N\Vert f_i(X_1,X_2,\dots ,X_N){\Vert }_F^2$$

$$\mathrm{subject}\mathrm{to}X\in M:={\mathbb{R}}^{n\times n}\times {\mathbb{R}}^{n\times n}\times \cdots \times {\mathbb{R}}^{n\times n}.$$

$$\mathrm{D}f(X)[(X^{\prime })]$$

$$\mathrm{grad}f(X)=(D_1,D_2,\dots ,D_N).$$

$$\mathrm{Hess}f(X)[({\xi }_1,{\xi }_2,\dots ,{\xi }_N)]:=\mathrm{D}\mathrm{grad}f(X)[({\xi }_1,{\xi }_2,\dots ,{\xi }_N)]=(D_1^{\prime },D_2^{\prime },\dots ,D_N^{\prime }),$$

$$p_{k+1}=p_k+t_k{d}_k.$$

$$⟨({\xi }_1,{\xi }_2,\dots ,{\xi }_N),({\eta }_1,{\eta }_2,\dots ,{\eta }_N)⟩:=\sum _{i=1}^N\mathrm{tr}({\xi }_i^T{\eta }_i),$$

$$f(p+{\beta }^{\gamma }\bar{\alpha }d)\le f(p)+\sigma ⟨\mathrm{grad}f(p),{\beta }^{\gamma }\bar{\alpha }d⟩$$

$$d_k=-\mathrm{grad}f(p_k)+{\beta }_k{d}_{k-1},$$

Full text

Numerical experiments illustrate that for large-scale problems the trust-region method is more effective than the steepest descent and the conjugate gradient methods. The fixed-point approach, however, is superior to the optimization methods. As an application we consider a networked control system, where the Markov jumps are induced by the wireless communication protocol.

1 Introduction

Markovian jump linear systems are expressed by

[TABLE]

where $x\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$, and $y\in\mathbb{R}^{p}$ are the state, input, and output, respectively. The parameter $r(t)$ denotes a continuous-time Markov process on a probability space, which takes values in a finite set $S:=\{1,2,\ldots,N\}$ with transition probabilities given by

[TABLE]

where $\delta>0$, and $\pi_{ij}$ denotes the transition probability rate from mode $i$ to mode $j$ when $i\neq j$. Furthermore, for all $i\in S$, $\pi_{ij}$ satisfies $\pi_{ij}\geq 0$ $(i\neq j)$ and $\pi_{ii}=-\sum_{j\in S,\,j\neq i}\pi_{ij}$. The Markov process $\{r(t),\,t>0\}$ is assumed to have an initial process $r(0)=(\mu_{1},\mu_{2},\ldots,\mu_{N})$. The matrices $A(r(t))$, $B(r(t))$, and $C(r(t))$ are contained in $\{A_{1},A_{2},\ldots,A_{N}\}$, $\{B_{1},B_{2},\ldots,B_{N}\}$, and $\{C_{1},C_{2},\ldots,C_{N}\}$, respectively, and if $r(t)=i\in S$, we have $A(r(t))=A_{i}$, $B(r(t))=B_{i}$, and $C(r(t))=C_{i}$. Applications to electric power systems have been considered in [16, 25], others are mentioned e.g. in the survey paper [23] or the monograph [6]. In subsection 5.3 we sketch an application to networked control system.

The system (1) is called mean-square stable if the solution $x(t)$ of $\dot{x}(t)=A(r(t))x(t)$ satisfies

[TABLE]

for any initial condition $x(0)=x_{0}$ and initial distribution for $r(0)=r_{0}$, where $\|\cdot\|_{2}$ denotes the Euclidean norm and $E(\cdot)$ the expected value. In [7], for mean-square stable Markovian jump linear systems (1), the $H^{2}$ norm has been defined as

[TABLE]

where $\|y_{s,i}\|_{2}:=\sqrt{\int_{0}^{\infty}E(y^{T}_{s,i}(t)y_{s,i}(t))dt}$ and $y_{s,i}$ is the output $\{y(t)|t>0\}$ when

• •

the input is given by $u(t)=e_{s}\delta(t)$, where $\delta(t)$ is the unit impulse, and $e_{s}$ is the $m$ dimensional unit vector formed by having $1$ at the $s$th position and zero elsewhere.

• •

$x(0)=0$ and $r(0)=i$.

For $N=1$, the definition reduces to the usual $H^{2}$ norm. If (1) is mean-square stable, then by [7]

[TABLE]

where $(P_{1},P_{2},\ldots,P_{N})$ and $(Q_{1},Q_{2},\ldots,Q_{N})$ are the unique solutions of the coupled Lyapunov equations

[TABLE]

for $i=1,2,\ldots,N$, respectively. Using the above $H^{2}$ norm concept, [7] and [24] have studied an $H^{2}$ optimal state-feedback control and an $H^{2}$ optimal model reduction method for Markovian jump linear systems (1), respectively. Thus it is important to study effective algorithms for solving the Lyapunov like equations (2) and (3). In a unified form we look for solutions $X=(X_{1},\ldots,X_{N})$ of coupled equations

[TABLE]

Here $A_{i},Y_{i}=Y_{i}^{T}\in\mathbb{R}^{n\times n}$ are constant matrices, and $\gamma_{ij}\geq 0$ $(i\neq j)$ and $\gamma_{ii}=-\sum_{j\in S,\,j\neq i}\gamma_{ij}$, or more generally, just $\gamma_{ii}<0$. Eqs. (2) and (3) are special cases of (4). If $\gamma_{ij}=0$ for all $i\neq j$, then (4) just consists of uncoupled standard Lyapunov equations. For the standard Lyapunov equation, several algorithms have been developed, see [10, 17], and references therein. Coupled Lyapunov equations of type (4) have been solved e.g. in [4, 5, 11, 27, 29, 15], partly by fixed-point iterations, partly by optimization methods. However, all methods were only applied to examples of small dimensions. In trying to extend methods of model order reduction to Markov jump linear systems, we found the necessity to develop a more efficient solver. In this paper we compare different classes of algorithms for their applicabilty at least to medium-sized problems with $nN\approx 1000$ or larger.

Our first and main approach follows the ideas in [10] and uses a fixed point iteration. Here the efficiency can be improved by considering Krylov-subspace methods and appropriate vectorization. The idea is rather simple, but works much better than other more intricate methods.

In our second approach, we reformulate the equation as an optimization problem on the product space of $N$ Euclidean matrix spaces and derive the gradient and the Hessian of the objective function. The gradient is used to develop a steepest descent and a conjugate gradient method; the Hessian is applied to establish a trust-region method.

Numerical experiments illustrate that all our methods work, but only the fixed-point iteration lends itself for larger problems.

2 Preliminaries

Let $\mathcal{H}^{n}$ denote the space of symmetric (i.e. real Hermitian) $n\times n$ matrices and $\mathcal{H}_{+}^{n}$ the cone of nonnegative definite matrices.

Then we set

[TABLE]

such that $\mathcal{H}$ is an ordered real vector space with ordering cone $\mathcal{H}_{+}$. For $A\in\mathbb{R}^{n\times n}$, we define the Lyapunov operator

[TABLE]

On $\mathcal{H}$ we consider the blockwise Lyapunov operator $\mathcal{L}:\mathcal{H}\to\mathcal{H}$, defined by

[TABLE]

and the positive operator $\Pi:\mathcal{H}\to\mathcal{H}$, $\Pi(\mathcal{H}_{+})\subset\mathcal{H}_{+}$, defined by

[TABLE]

From [6, Thm. 3.15] and also [21, 9] we cite a well-known stability result.

Theorem 2.1

The following statements are equivalent.

• (a)

System (1) is asymptotically mean-square stable.

• (b)

\sigma(\mathcal{L}+\Pi)\subset\mathbb{C}_{-}=\left\{\lambda\in\mathbb{C}\;\big{|}\;\Re\lambda<0\right\}.

• (c)

$\sigma(\mathcal{L})\subset\mathbb{C}_{-}$* and $\rho(\mathcal{L}^{-1}\Pi)<1$.*

• (d)

$\exists X\in\mathcal{H}_{+}$: $(\mathcal{L}+\Pi)(X)<0$.

• (e)

$\forall Y\in\mathcal{H}_{+}$:$\exists X\in\mathcal{H}_{+}$: $(\mathcal{L}+\Pi)(X)=-Y$.

The coupled Lyapunov equations (4) can be written in the form

[TABLE]

where $Y\in\mathcal{H}$. Under the assumptions of asymptotic mean-square stability there exists a unique solution $X\in\mathcal{H}$. Moreover, if $Y\in\mathcal{H}_{+}$, then $X\in\mathcal{H}_{+}$. Since $X$ contains $\frac{1}{2}{Nn(n+1)}$ scalar unknowns, a direct solution e.g. via Kronecker-product representation and Gaussian elimination has complexity $O(N^{3}n^{6})$.

However, it is well-known, that a single Lyapunov equation for an unknown $n\times n$-matrix can be solved with $O(n^{3})$ operations by the Bartels-Stewart-algorithm [2]. Since $\mathcal{L}^{-1}(X)$ for $X\in\mathcal{H}$ is obtained by solving $N$ independent Lyapunov equations, the cost of evaluating $\mathcal{L}^{-1}(X)$ is only $O(Nn^{3})$. Based on this observation we suggest some fixed point iterations.

3 Fixed point formulations and Krylov subspace methods

3.1 Jacobi and Gauss-Seidel schemes

If $\rho(\mathcal{L}^{-1}\Pi)<1$ then the solution to (5) is obtained as the limit of the iterative scheme

[TABLE]

Blockwise, we get the Jacobi-type fixed point iteration

[TABLE]

In the $i$-th row, it is natural to replace $X_{j}^{(k)}$ by $X_{j}^{(k+1)}$ for $j<i$, because this update is already available. Thus we get the Gauss-Seidel-type scheme

[TABLE]

3.2 Preconditioned Krylov subspace iterations

Both schemes can be written in the form $X^{(k+1)}=T(X^{(k)})$ and interpreted as preconditioners for Krylov subspace iterations as has been explained e.g. in [10, Sec 4.2]. The basic idea is to find an optimal approximation to the solution within the Krylov subspace $X^{(0)}+\operatorname{span}\{X^{(1)},\ldots,X^{(k)}\}$. In the Jacobi formulation, this means that we set $T_{J}=-\mathcal{L}^{-1}\Pi$ and apply some standard Krylov subspace method to the equation

[TABLE]

That is, we replace $-Y$ by $\tilde{Y}=-\mathcal{L}^{-1}(Y)$ and then solve with the linear mapping $I-\mathcal{L}^{-1}\Pi$. In the Gauss-Seidel formulation, the linear mapping $T_{\text{GS}}:X\to\tilde{X}$ is described by the scheme

[TABLE]

The update of the right hand side $Y\mapsto\tilde{Y}$ is obtained via

[TABLE]

Thus we can apply a Krylov subspace method to the equation $(I-T_{\text{GS}})(X)=\tilde{Y}$.

3.3 Avoiding loops by vectorization

Terms of the form $\sum_{j=1}^{N}\gamma_{ij}X_{j}$ can be vectorized efficiently. If $\operatorname{vec}X_{i}$ as before denotes the vector obtained by stacking all columns of $X_{i}$ one above the other, and (by a slight abuse of notation) $\operatorname{vec}X=[\operatorname{vec}X_{1},\ldots,\operatorname{vec}X_{N}]$, then

[TABLE]

where $\Gamma_{i}$ is the $i$-th row of $\Gamma$. Thus, we can write the sum in (8) as a matrix-vector product, which is processed faster than a loop over $j$.

3.4 An implementation

Some of the specific vectorization ideas can be seen more clearly in the MATLAB®-listing which we add for convenience.

⬇

1ction X=MJLSlyap(A,Y,G,tol)

2% Solve algebraic Lyapunov equation for Markov jump linear system

3% The n^2 x N matrices A and Y contain A1(:),...,AN(:) and Y1(:),...,YN(:)

4% as columns. The matrix G=Gamma is N x N

5% Krylov subspace approach with Gauss-Seidel-type preconditioning

6

7N=size(G,2);n2=size(A,1);n=sqrt(n2);

8d=diag(G);G=G'-diag(d);d=d/2;

9A=reshape(A+reshape(eye(n),n2,1)*d',n,n,N);

10Y=reshape(Y,n,n,N);

11for i=1:N

12 Yi=Y(:,:,i)+reshape(reshape(Y(:,:,1:i),n2,i)*G(1:i,i),n,n);

13 Y(:,:,i)=lyap(A(:,:,i),Yi);

14end

15X=reshape(bicgstab(@T,Y(:),tol),n,n,N);

16 function Z = T(X)

17 X=reshape(X,n,n,N);Z=zeros(size(X));

18 for j=1:N

19 Xj=lyap(A(:,:,j),reshape(reshape(X,n2,N)*G(:,j),n,n));

20 Z(:,:,j)=X(:,:,j)-Xj;

21 X(:,:,j)=Xj;

22 end

23 Z=Z(:);

24 end

25end

4 Optimization-based approach

This section formulates an optimization problem for solving the Lyapunov like equations (4) and suggests a number of methods to solve this problem.

4.1 Reformulation as optimization problem

To solve (4), we consider the following optimization problem.

Problem 4.1

[TABLE]

If $f(X)\approx 0$, we obtain $f_{i}(X)\approx 0$ $(i=1,2,\ldots,N)$, i.e., $X$ is an approximate solution of (4).

4.2 Gradient and Hessian of the objective function $f$

To develop optimization algorithms for Problem 4.1, this section derives the gradient of the objective function $f$, and then gives the Hessian of $f$.

The Fréchet derivative of $f$ at $X=(X_{1},X_{2},\ldots,X_{N})$ in the direction $X^{\prime}=(X^{\prime}_{1},X^{\prime}_{2},\ldots,X^{\prime}_{N})$ can be calculated as

[TABLE]

for \displaystyle D_{i}:=2\Big{(}A_{i}^{T}f_{i}+f_{i}A_{i}+\sum_{j=1}^{N}\gamma_{ji}f_{j}\Big{)}.

Since the gradient ${\rm grad}f(X)$ satisfies $\displaystyle{\rm D}f(X)[(X^{\prime})]=\sum_{i=1}^{N}{\rm tr}(X^{\prime T}_{i}({\rm grad}f(X))_{X_{i}})$, equation (9) yields

[TABLE]

Furthermore, from (4) and (10), the Hessian of $f$ is given by

[TABLE]

where \displaystyle D^{\prime}_{i}=2\Big{(}A_{i}^{T}f^{\prime}_{i}+f^{\prime}_{i}A_{i}+\sum_{j=1}^{N}\gamma_{ji}f^{\prime}_{j}\Big{)} with $f^{\prime}_{i}=A_{i}\xi_{i}+\xi_{i}A_{i}^{T}+\sum_{j=1}^{N}\gamma_{ij}\xi_{j}$.

4.3 Optimization algorithms for Problem 4.1

In the vector space $M$, optimization methods based on line search can be developed. In such methods, given the current point $p_{k}\in M$, the search direction $d_{k}\in T_{p_{k}}M\cong M$, and the step size $t_{k}>0$, the next point $p_{k+1}\in M$ is computed as

[TABLE]

We consider three optimization algorithms for Problem 4.1. To this end, for any $(\xi_{1},\xi_{2},\ldots,\xi_{N}),(\eta_{1},\eta_{2},\ldots,\eta_{N})\in T_{p}M\cong M$, we define the inner product as

[TABLE]

and the induced norm by $\|(\xi_{1},\xi_{2},\ldots,\xi_{N})\|:=\sqrt{\langle(\xi_{1},\xi_{2},\ldots,\xi_{N}),(\xi_{1},\xi_{2},\ldots,\xi_{N})\rangle}.$

4.3.1 Steepest descent method for Problem 4.1

In the steepest descent method, the negative gradient of the objective function $f$ at a current iterate $p_{k}\in M$ can be chosen as a search direction $d_{k}\in T_{x_{k}}M$ at $p_{k}$, i.e., $d_{k}:=-{\rm grad}\,f(p_{k})$. As a step size $t_{k}$, the following Armijo step size is popular [18]: Given a point $p\in M$, a tangent vector $d\in T_{p}M$, and scalars $\bar{\alpha}>0,\beta,\sigma\in(0,1)$, the Armijo step size $t^{A}:=\beta^{\gamma}\bar{\alpha}$ is defined in such a way that $\gamma$ is the smallest nonnegative integer satisfying

[TABLE]

Algorithm 2 describes the steepest descent method for solving Problem 4.1.

4.3.2 Conjugate gradient method for Problem 4.1

In the conjugate gradient method, the search direction $d_{k}$ at the current point $p_{k}$ is computed as

[TABLE]

where $\beta_{k}>0$. The Dai and Yuan type parameter $\beta_{k}$ is given by

[TABLE]

In more detail, see [8, 20]. To guarantee the convergence property, we use the step size $t_{k}$ satisfying

[TABLE]

where $0<c_{1}<c_{2}<1$. The condition (15) is called the Wolfe condition. Algorithm 3 describes the conjugate gradient method for solving Problem 4.1.

4.3.3 Trust-region method for Problem 4.1

At each iterate $p$ in the trust-region method on the vector space $M$, we evaluate the quadratic model $\hat{m}_{p}$ of the objective function $f$ within a trust region:

[TABLE]

for $d\in T_{p}M$, where ${\rm grad}\,f(p)$ and ${\rm Hess}\,f(p)[d]$ are given by (10) and (11), respectively. A trust-region with a radius $\Delta>0$ at $p\in M$ is defined as a ball with center [math] in $T_{p}M$. Thus the trust-region subproblem at $p\in M$ with a radius $\Delta$ is defined as a problem of minimizing $\hat{m}_{p}(d)$ subject to $d\in T_{p}M$, $\|d\|\leq\Delta$. This subproblem can be solved by the truncated conjugate gradient method [1]. Then we compute the ratio of the decreases in the objective function $f$ and the model $\hat{m}_{p}$ attained by the resulting $d_{*}$ to decide whether $d_{*}$ should be accepted and whether the trust-region with the radius $\Delta$ is appropriate. Algorithm 4 describes the process. The constants $\frac{1}{4}$ and $\frac{3}{4}$ in the condition expressions in Algorithm 4 are commonly used in the trust-region method for a general unconstrained optimization problem. These values ensure the convergence properties of the algorithm [1].

4.4 Stopping criterion for Algorithms 2, 3, and 4

In practice, we need a stopping criterion for Algorithms 2, 3, and 4. In this paper, we stop the algorithm when $\|{\rm grad}\,f(p_{k})\|<10^{-5}$ or the iteration number reached 30000.

5 Numerical experiments and an application

We first discuss some artificial and random examples to compare the performance of Algorithms 1–4. It turns out that Algorithm 1 is our method of choice. Then we provide some details on an application of Markov jump linear systems to networked control systems. Our method works fine for this example as well. All computations were carried out using MATLAB® R2016b on an Intel® Core(TM) i7-7500U CPU @ 2.70GHz 2.90GHz and 16.0 GB RAM. We report on computation times, numbers of required iterations, residuals and – where available – the absolute errors. Note that for Algorithm 1 each evaluation of $T_{\text{GS}}$ is half an iteration step. The residual is given as the square root of the objective function $f$.

5.1 A system with known solution

Our first example is constructed with known solution. Let $\begin{pmatrix}\gamma_{11}&\gamma_{12}\\ \gamma_{21}&\gamma_{22}\end{pmatrix}=\begin{pmatrix}-1&1\\ 2&-2\end{pmatrix}$ and

[TABLE]

Then $Y_{j}=-(\mathcal{L}+\Pi)(X)$ is given by

[TABLE]

For Algorithm 1, 2, 3, and 4, Table 1 shows computational times, iteration numbers, the residual and the error $\|\Delta\|_{F}$, i.e. the deviation of the computed solution from the known solution, in the Frobenius norm. According to Table 1, Algorithm 1 is superior to the others.

5.2 Random systems

We consider random mean square asymptotically systems. In our first experiments we always set $\Gamma=\begin{pmatrix}-0.3&0.3\\ 0.5&-0.5\end{pmatrix}$ and generate $A_{1},A_{2},Y_{1},Y_{2}\in\mathbb{R}^{n\times n}$ such that $-A_{1},-A_{2},Y_{1},Y_{2}$ are symmetric positive definite. Under this assumption we have $(\mathcal{L}+\Pi)([I,I])=-2([A_{1},A_{2}])<0$ implying asymptotic mean-square stability.

Table 2 shows the computational times, iteration numbers, and the residuals $\sqrt{f(X_{1},X_{2})}$ in Algorithms 1, 2, 3, and 4 for different $n$. For $n=100$ Algorithms 2 and 3 turn out to be inadequate and for larger dimensions Algorithm 4 becomes too slow.

In the next experiments we generate random nonsymmetric matrices $A_{j}\in\mathbb{R}^{n\times n}$, $j=1,\ldots,N$ with $\sigma(A_{j})\subset\mathbb{C}_{-}$ and $\Gamma=(\gamma_{ij})\in\mathbb{R}^{N\times N}$ with $\gamma_{ii}<0$ and $\gamma_{ij}\geq 0$ for $i,j\neq 0$. Mean square stability of the system is enforced by appropriate scaling of the off-diagonal entries $\gamma_{ij}$. In these experiments, whose results are presented in Table 3, only Algorithm 1 is considered.

5.3 A communication network example

It has been argued in the literature (e.g. [19, 14]) that Markov jump linear systems can be used to model communication phenomena in cyber-physical systems. In the following we sketch such an example, where we have a fixed number $\nu$ of entities or agents to be controlled. Each entity is seen as a transmitting station that transmits its observed values via a medium which follows the Carrier Sense Multiple Access/Collision Avoidance (CSMA/CS) principle, as it is used e.g. in a WLAN transmission. A Markov jump linear system can then be used to determine which station is allowed to send its data.

In a CSMA/CA network after every transmission a backoff value is assigned to the station which did just transmit. The backoff value is a uniformly distributed random variable in a given interval [0,ContentionWindowMax]. While the physical medium of a CSMA/CS network is idle, the backoff value of each station is reduced continuously. When the backoff value reaches [math] the respective station is again allowed to transmit its data. See figure 1 for a visualization.

We then keep track of the $\tau$ last transmissions and whether the last transmission was faulty or not. This information can then be encoded as states in a transition rate matrix for a Markov chain. Assuming $2$ entities and a transition rate matrix which keeps track of the last $3$ transmissions this leads to a transition rate matrix of size $16\times 16$, which we present here in Table 4.

In the matrix we see the probabilities to get from the current situation, presented in the row of the matrix, to the future situation presented in the columns. We denote by $e=0$ that the transmission is correct and $e=1$ that the transmission is incorrect. The values for $t$ denote which station sends at the moment, $t-1$ which station did send before and $t-2$ which station did send before $t-1$. The station sending next is then denoted by $t+1$. The probability for an error in the next transmission is set to $0.03$, if the current transmission is correct, and $0.75$ if we have an error in the current transmission, as e.g. suggested in [13], see also [12].

Now we want to approximate the probability for a fixed station to be the next to send. Let

[TABLE]

where for each station, the set $J$ contains the index of its last occurrence in the memory. Then the approximated probability for a fixed station to be the next to send (neglecting the probability of an error), is either $\bar{w}$, if it did not send in the last $\tau$ transmissions, or $\frac{\bar{w}}{\tau+1-\left(i-1\right)}$, if its last transmission was at $t+1-i$. The transition matrix $\Theta$ in Table 4 contains these entries. As an example consider the transition 0AAA $\leadsto$ 0BAA. Here $J=\{1\}$ because the first position of AAA contains the station A. Hence $\bar{w}=(2-1+\frac{1}{3+1-(1-1)})^{-1}=\frac{4}{5}$. Thus the transition probability is $0.8\cdot 0.97\approx 0.78=\Theta_{1,5}$. For a different case consider 0BBA $\leadsto$ 0ABB. Here $J=\{1,3\}$ because the first occurences of B and A in BBA are at positions $1$ and $3$. Here $\bar{w}=(2-2+\frac{1}{3+1-(1-1)}+\frac{1}{3+1-(3-1)})^{-1}=\frac{4}{3}$ and the transition probability is $\frac{4}{3}\cdot\frac{1}{2}\cdot 0.97\approx 0.65=\Theta_{7,4}$. The formula was obtained by a discretization of the interval [0,ContentionWindow] and an approximation for the exact probability for each station to send next. Details of the technical derivation can be found in [26]. In the case of an error, the transmission memory is not changed until the error state changes. So for arbitrary X,Y,Z, only 1XYZ$\leadsto$1XYZ or 1XYZ$\leadsto$0XYZ is possible, with the corresponding probabilities. Note that we have thus approximated the transition probabilities for a discrete time setting. The transformation to the continuous time situation is obtained by $\Gamma=a(\Theta-I)$ where $a$ is the average time spent in one transmission mode, see [28].

Together with this model of communication, we consider the simple system described in figure 2.

We denote with $e\in\left\{0,1\right\}$ whether or not an error ocurred in the last transmission and with $j\in\left\{1,2,3\right\}$ the station which did send last and its corresponding diagonal block. Note that both $e$ and $j$ are determined by $r\left(t\right)$ but are added here for easier understanding of the model. The state vector of the whole system is $\begin{pmatrix}{s_{1}\left(t\right)}&{v_{1}\left(t\right)}&{s_{2}\left(t\right)}&{v_{2}\left(t\right)}&{s_{3}\left(t\right)}&{v_{3}\left(t\right)}\end{pmatrix}^{T}$. The linearized dynamics are given by the equations

[TABLE]

Here $I_{2}$ denotes the identity matrix and $0_{2}$ the [math] matrix in $\mathbb{R}^{2\times 2}$, while $m,g,R$ denote the mass of each cart, the gravitational acceleration, and the friction coefficient, respectively. In this application only the matrices $C_{i}$ are depending on the Markov process. There is still a need for fast solving of Lyapunov like equations in this applications since the size of the transition rate matrix, and therefore the number of coupled Lyapunov equations, scales with a factor $\nu^{\tau}$, where $\nu$ is the number of entities and $\tau$ is the number of tracked transmissions.

We compute the observability Gramian (3). Let $\nu=3$, $\tau=3$, $m=1$, $g=9.81$, and $R=0.1$. Table 5 shows the computational time, iteration number, and the residual in Algorithm 1. The results show that the algorithm works well in practical applications with relevant dimensions.

6 Conclusion

We have compared optimization based methods and a preconditioned Krylov subspace iteration for the solution of Lyapunov equations related to Markov jump linear systems. From our numerical experiments we deduce that only the Krylov subspace iteration lends itself for large systems. As an application we have sketched a Markov jump linear system model for a networked control system with WLAN based communication. In ongoing research, we plan to elaborate further on this model. We expect that efficient Lyapunov solvers will be an essential tool. A next step in this direction could be the development of low-rank methods like they were considered e.g. in [3, 22] for other types of Lyapunov equations.