The LMI Approach for Stabilizing of Linear Stochastic Systems

Abstract

Stochastic linear systems subjected both to Markov jumps and to multiplicative white noise are considered. In order to stabilize such type of stochastic systems, the so-called set of generalized discrete-time algebraic Riccati equations has to be solved. The LMI approach for computing the stabilizing symmetric solution (which is in fact the equilibrium point) of this system is studied. We construct a new modification of the standard LMI approach, and we show how to apply the new modification. Computer realizations of all modifications are compared. Numerical experiments are given where the LMI modifications are numerically compared. Based on the experiments the main conclusion is that the new LMI modification is faster than the standard LMI approach.

1. Introduction

In this paper we investigate general stochastic algebraic Riccati equations which are related to LQ control models for stochastic linear systems with multiplicative white noise and markovian jumping. We consider a stochastic system described by 𝑑𝑦(𝑡)=𝐴₀𝜂(𝑡)𝑦(𝑡)+𝐵₀𝜂(𝑡)𝑢(𝑡)+𝑑𝑡_𝑘𝑗=1𝐴_𝑗𝜂(𝑡)𝑦(𝑡)+𝐵_𝑗𝜂(𝑡)𝑢(𝑡)𝑑𝑊_𝑗(𝑡)𝑦(0)=𝑦₀,𝑧(𝑡)=𝐶𝜂(𝑡)𝑥(𝑡)+𝐷𝜂(𝑡)𝑢(𝑡), (1) where𝑊(𝑡)=(𝑊₁(𝑡),…,𝑊_𝑘(𝑡))^𝑇is a 𝑘-dimensional standard Brownian motion with𝑡∈[0,+∞)and𝑊(0)=0, defined by a filtered probability space(Ω,ℱ,ℱ_𝑡,𝑃). In addition, 𝜂(𝑡), 𝑡≥0is a right continuous homogeneous Markov chain with the state space the set𝒟={1,…,𝑁}and the probability transition matrix𝑃(𝑡)=(𝑝_𝑖𝑗(𝑡))=𝑒^𝑄𝑡, 𝑡≥0, and𝑄=(𝑞_𝑖𝑗)with∑^𝑘_𝑗=1𝑞_𝑖𝑗=0, 𝑖∈𝒟, and 𝑞_𝑖𝑗≥0if 𝑖=^╱𝑗. It is assumed that the 𝑊(𝑡) and𝜂(𝑡)are independent stochastic processes and𝑃{𝜂(0)=𝑖}>0for all𝑖∈𝒟. The state vector 𝑥 is a𝑛×1real vector, 𝑢 denotes the vector of 𝑚 control variables, and 𝑧 is the regulated output vector with 𝑝 components. The matrix coefficients 𝐴_𝑗,𝐵_𝑗,0≤𝑗≤𝑘, 𝐶(𝑖), 𝐷(𝑖), 𝑖∈𝒟are constant matrices of appropriate dimensions with real elements.

The stochastic systems with multiplicative white noise naturally arise in control problems of linear uncertain systems with stochastic uncertainty. It is important for applications to find a stabilizing controller for the above stochastic system (for more details see [1, 2]). For this purpose it is enough to compute the stabilizing solution to the following stochastic generalized Riccati algebraic equations: 𝒯_𝑖(𝐗)∶=𝐴₀(𝑖)^𝑇𝑋(𝑖)+𝑋(𝑖)𝐴₀(𝑖)+^𝑘_𝑗=1𝐴_𝑗(𝑖)^𝑇𝑋(𝑖)𝐴_𝑗(𝑖)+𝑄(𝑖)+^𝑁_𝑗=1𝑞_𝑖𝑗𝑋𝑗−𝑆(𝑖,𝑋)^𝑇[𝑅(𝑖,𝑋)]⁻¹𝑆()()𝑅()(𝑖)𝑖,𝑋=0with𝑅𝑖,𝑋>0,𝑖∈𝒟,𝑖,𝑋=𝐷^𝑇𝐷(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑋(𝑖)𝐵_𝑗(𝑖),𝑄(𝑖)=𝐶^𝑇(𝑖)𝐶(𝑖),𝑆(𝑖,𝑋)=𝐵₀(𝑖)^𝑇𝑋(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑋(𝑖)𝐴_𝑗(𝑖)+𝐷(𝑖)^𝑇𝐶(𝑖). (2)(3)

The concepts of stabilizing solution of the Riccati-type equation (2) and of stabilizability for the triple(𝐀,𝐁,𝑄), where as usual𝐀=(𝐴₀,𝐴₁,…,𝐴_𝑘),𝐁=(𝐵₀,𝐵₁,…,𝑏_𝑘), are defined in a standard way (see [1]).

Applying Theorem 4.9 of Dragan and Morozan [2] we deduce that system (2) has a unique stabilizing solution(𝑋(1),𝑋(2),…,𝑋(𝑁))with [𝐷(𝑖)^𝑇𝐷(𝑖)+∑^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑋(𝑖)𝐵_𝑗(𝑖)]>0. The control 𝑢(𝑡)=−𝑅𝑖,𝑋⁻¹𝑆𝑖,𝑋𝑥(𝑡) (4) stabilizes system (1). An e1ective iterative convergent algorithm to compute these stabilizing solutions is presented in [3].

Let us consider the special case of (1) where𝒟={1}. The stochastic linear quadratic model studied by Yao et al. in [4] is obtained. In this model the special functional𝐸∫^∞₀[𝑥(𝑡)^𝑇𝑄𝑥(𝑡)+𝑢(𝑡)^𝑇𝑅𝑢(𝑡)]𝑑𝑡is minimized (see equation (8) from [4]). The matrices𝑄and𝑅are so-called cost weighting matrices. It is allowed that the cost weighting matrices for the state and the control are singular, that is, 𝑄=𝐶^𝑇𝐶and𝑅=𝐷^𝑇𝐷are singular matrices. Such type models belong to a wide class of indefinite SLQ models. Applications of the indefinite LQ problems can be found in [5] for pollution control, [4, 6–9] in which problem appears in field of the mathematical finance. In this case the stochastic LQ problem can be solved via the following stochastic algebraic Riccati equation (the symmetric matrix𝑃being the unknown one): 𝒯(𝑃)∶=𝐴^𝑇𝑃+𝑃𝐴+𝑄+^𝑘_𝑗=1𝐶_𝑗^𝑇𝑃𝐶_𝑗−𝐵^𝑇𝑃+^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐶_𝑗^𝑇×𝑅+^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐷_𝑗⁻¹𝐵^𝑇𝑃+^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐶_𝑗=0with𝑅+^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐷_𝑗>0. (5)

The coefficient matrices of (5)𝐴, 𝐵, 𝑄, 𝑅, 𝐶_𝑗, and𝐷_𝑗for𝑗=1,…,𝑘are given ones of sizes𝑛×𝑛,𝑛×𝑝,𝑛×𝑛,𝑚×𝑚,𝑛×𝑛, and𝑛×𝑚, respectively, and𝑄denotes the state matrix, and𝑅is the control matrix. For a deterministic case (𝐶_𝑗=𝐷_𝑗=0,𝑗=1…,𝑘) it is assumed that the matrix𝑄 is positive semidefinite and𝑅is a positive definite one.

Here the computation of stabilizing solution of system (1) is explicitly expressed in terms of the solutions of some linear matrix inequalities (LMI). The paper is devoted to the LMI approach and its modifications. The LMI approach is very important for the practice and real-world problems. Very often the LMI approach is an only method for solving a given class of problems. The application of the LMI approach to the solution of the optimal control problems is studied in [10–14]. We introduce a modification set of nonlinear equations equivalent to (1) which lead us to the new convex optimization problems. The LMI approach applied to the new optimization problem gives a fast way to find the stabilizing solution to (1). We will compare the numerical effectiveness of the introduced LMI solvers. Numerical simulations are used to demonstrate the performance of the considered solvers.

The notations used in this note are standard. Here,𝑋≥0 (or𝑋>0); it is denoted that𝑋=𝑋^𝑇is positive semidefinite (or positive definite) and‖⋅‖₂denotes the spectral matrix norm.

2. The LMI Approach to the Generalized Stochastic Riccati Equations

Rami and Zhou [15] have investigated the stochastic algebraic Riccati equation (5) in case𝑘=1, where the numerical method to compute the maximal solution to (5) with𝑘=1is derived. The method is based on the solution of a convex optimization problem over linear matrix inequalities. The LMI approach is considered as powerful tool in optimization. Further on, the authors in [4] have developed a computational approach to such SLQ models(𝑘>1)using an LMI formulation. The LMI optimization is a successful method for solving (5) in this case. Although we cannot solve the stochastic algebraic Riccati equation we can still find the optimal control law via the LMI approach. It is well known [16] that if the stochastic system is stabilizable, in the mean square sense, the LMI optimization method always yields the maximal positive semidefinite solution to (5). The problem of stability and optimality of SLQ model (5) with𝑘=1(the case of one-dimensional Brownian motion) is treated by [7]. Moreover, this study is continued in terms of the multidimensional model (𝑘>1) by [4]. The optimization problem associated with (5) is max𝐼_𝑛×,𝑃𝑄+𝐴^𝑇+𝑃+𝑃𝐴_𝑘𝑗∑₌₁𝐶_j^𝑇𝑃𝐶_𝑗𝐵^𝑇𝑃+∑^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐶_𝑗^𝑇𝐵^𝑇𝑃+∑^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐶_𝑗𝑅+∑^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐷_𝑗≥0𝑅+^𝑘_𝑗=1𝐷_𝑗^𝑇𝑃𝐷_𝑗>0𝑃=𝑃^𝑇, (6) with respect to the variable𝑃which is a symmetric matrix.

Under notations 𝑊(𝑖,𝑋)=𝐴₀(𝑖)+0.5𝑞_𝑖𝑖𝐼^𝑇𝑋(𝑖)+𝑋(𝑖)𝐴₀(𝑖)+0.5𝑞_𝑖𝑖𝐼=𝐴₀(𝑖)^𝑇𝑋(𝑖)+𝑋(𝑖)𝐴₀(𝑖),Π₁(𝑖,𝑋)=^𝑘_𝑗=1𝐴_𝑗(𝑖)^𝑇𝑋(𝑖)𝐴_𝑗(𝑖)+_𝑗╱=𝑖𝑞_𝑖𝑗𝑋𝑗, (7) the optimization problem associated with the set of equations (2) is max^𝑁_𝑖=1𝐼_𝑛(𝑖×𝑊(),𝑋)subjectto𝑖∈𝒟𝑖,𝑋+Π₁(𝑖,𝑋)+𝐶(𝑖)^𝑇𝐶(𝑖)𝑆(𝑖,𝑋)^𝑇𝑆()𝑅()𝑅()𝑋(𝑖)(𝑖)𝑖,𝑋𝑖,𝑋≥0𝑖,𝑋>0=𝑋^𝑇. (8)

In this paper we investigate the numerical solvability of the semidefinite programming problem (8) for different types of matrices𝐷(𝑖)^𝑇𝐷(𝑖), 𝑖∈𝒟. However, the numerical experiments for finding the maximal solution of (2) show that the LMI method (8) is slowly working for different types of matrices𝐷(𝑖)^𝑇𝐷(𝑖)in the case 𝑘=1. Here we introduce a new modification to accelerate the LMI method for solving the optimization problem (8). In many applications of control system theory the following fact is exploited. If any matrix𝑅=𝑅^𝑇is singular or zero and the matrix𝐵has the full rank, then there exists a symmetric matrix 𝑍 such that𝑅+𝐵^𝑇𝑍𝐵is a positive definite one. In our investigation the matrices𝐵(𝑖)=^{𝐵1𝐵(𝑖)2}_⋮𝐵^(𝑖)_𝑘(𝑖)have the full rank. Thus, we replace𝑋(𝑖)=𝑍(𝑖)+𝑌(𝑖)in𝑅(𝑖,𝑋). It is obtained 𝑅(𝑖,𝑍+𝑌)=𝐷(𝑖)^𝑇𝐷(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇(𝑍(𝑖)+𝑌(𝑖))𝐵_𝑗(𝑖)=𝑅(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑌(𝑖)𝐵_𝑗(𝑖), (9) where𝑌(𝑖)is the new unknown matrix. We apply this conclusion to set of equations (2). We construct the matrices𝑍(1),…,𝑍(𝑁)such that all matrices 𝑅(𝑖)=𝐷(𝑖)^𝑇𝐷(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑍(𝑖)𝐵_𝑗(𝑖) (10) are positive definite ones for𝑖∈𝒟. Then, the following set of Riccati equations is obtained regarding 𝑌(1),…,𝑌(𝑁): 𝒲_𝑖(𝐘)∶=𝑊(𝑖,𝑌)+𝑄(𝑖)−𝑆(𝑖,𝑌)^𝑇𝑅(𝑖,𝑌)⁻¹×𝑆(𝑖,𝑌)=0,𝑖∈𝒟, (11) where 𝑄(𝑖)=𝐶(𝑖)^𝑇𝐶(𝑖)+𝐴₀(𝑖)^𝑇𝑍(𝑖)+𝑍(𝑖)𝐴₀(𝑖)+^𝑘_𝑗=1𝐴_𝑗(𝑖)^𝑇𝑍(𝑖)𝐴_𝑗(𝑖)+_𝑗╱=𝑖𝑞_𝑖𝑗𝑍𝑗𝑅(𝑖,𝑌)=𝑅(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑌(𝑖)𝐵_𝑗(𝑖)𝑆(𝑖,𝑌)=𝐵₀(𝑖)^𝑇𝑌(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑌(𝑖)𝐴_𝑗(𝑖)+𝐿(𝑖),𝐿(𝑖)=𝐵₀(𝑖)^𝑇𝑍(𝑖)+^𝑘_𝑗=1𝐵_𝑗(𝑖)^𝑇𝑍(𝑖)𝐴_𝑗(𝑖)+𝐷(𝑖)^𝑇𝐶(𝑖). (12) Thus, the new optimization problem over LMI conditions related to (11) is derived: max^𝑁_𝑖=1𝐼_𝑛(𝑖×𝑊(),𝑌)subjectto𝑖∈𝒟𝑖,𝑌+Π₁(𝑖,𝑌)+𝑄(𝑖)𝑆(𝑖,𝑌)^𝑇𝑆()𝑅()𝑅()𝑌(𝑖)(𝑖)𝑖,𝑌𝑖,𝑌≥0𝑖,𝑌>0=𝑌^𝑇. (13)

Our experience on computations with LMI approach shows that there are examples where 𝑌(1),…,𝑌(𝑁)are not positive semidefinite, but the matrices𝑅(𝑖,𝑌)are still positive definite. In this reason we can consider a practical implementation of (13) where the constrain𝑅(𝑖,𝑌)>0is omitted. This leads us to the following problem: max^𝑁_𝑖=1𝐼_𝑛(𝑖×𝑊(),𝑌)subjectto𝑖∈𝒟𝑖,𝑌+Π₁(𝑖,𝑌)+𝑄(𝑖)𝑆(𝑖,𝑌)^𝑇𝑆()𝑅()𝑌(𝑖)(𝑖)𝑖,𝑌𝑖,𝑌≥0=𝑌^𝑇. (14)

3. Numerical Experiments

We carry out numerical simulations to present the numerical behaviour of introduced methods. In our experiments we apply three semidefinite programming problems (8), (13), and (14) for solving the stochastic algebraic Riccati equations (2). Our experiments are executed in MATLAB on an 2,16 GHz Intel(R) Dual CPU computer. The solutions of above optimization problems are obtained under the MATLAB lmi solvers which are executed with relative accuracy tol=1.0𝑒−9.

For all examples we take𝑄(1)=𝑄(2)=𝑄(3)to be diagonal matrices with entries(1,1,…,1,0,0). We have executed a set of examples with different values of𝑛and constant weighting matrices𝑅(1)=𝑅(2)=𝑅(3)=zeros(3,3). We compare all iterations introducing the following parameters: “m It”—the biggest number of iterations, “av It”— the average number of iterations. To determine the numbers “m It” and “av It” we count those examples of each size for which the corresponding iteration converges.

We consider a family of examples in case𝑁=3, 𝑘=2, and 𝑛=8,…,12,15,20, where the coefficient real matrices are given as follows:𝐴₀(𝑖), 𝐴₁(𝑖), 𝐴₂(𝑖), 𝐵₀(𝑖), 𝐵₁(𝑖), 𝐵₂(𝑖), 𝐿(𝑖), 𝑖=1,2,3 were constructed using the MATLAB notations 𝐴₀(1)=randn(𝑛,𝑛)();𝐴10−eye𝑛,𝑛₀(2)=randn(𝑛,𝑛)();𝐴10−eye𝑛,𝑛₀(3)=randn(𝑛,𝑛)();𝐴10−eye𝑛,𝑛₁(1)=randn(𝑛,𝑛);𝐴10₁(2)=randn(𝑛,𝑛);𝐴10₁(3)=randn(𝑛,𝑛);𝐴10₂(1)=randn(𝑛,𝑛);𝐴10₂(2)=randn(𝑛,𝑛);𝐴10₂(3)=randn(𝑛,𝑛);𝐵10₀(1)();𝐵=2∗rand𝑛,3₀(2)();𝐵=2∗rand𝑛,3₀(3)();𝐵=2∗rand𝑛,3₁(1)=randn(𝑛,3);𝐵100₁(2)=randn(𝑛,3);𝐵100₁(3)=randn(𝑛,3);𝐵100₂(1)=randn(𝑛,3);𝐵100₂(2)=randn(𝑛,3);𝐵100₂(3)=randn(𝑛,3);𝐿(1)();𝐿(2)();𝐿(3)().100=zeros𝑛,3=zeros𝑛,3=zeros𝑛,3 (15)

In our definitions the functions randn(p,k) and rand(p,k) return a 𝑝-by-𝑘 matrix of pseudorandom scalar values (for more information see the MATLAB description). The following transition probability matrix (see [3]) 𝑞_𝑖𝑗=0.3330.60.60.3330.30.10.3330.10.3 (16) is applied for all examples.

For our purpose we have executed hundred examples of each value of 𝑛 for the test. The maximal number of iterations “m It” and average number of iterations “av It” of each size for all examples needed for achieving the relative accuracy are reported in the table. There are three columns where the maximal errors 𝐸_𝑇=max_1,…,100max_𝑖=1,2,3‖𝒯_𝑖(𝐗)‖₂ and 𝐸_𝑊=max_1,…,100max_𝑖=1,2,3‖𝒲_𝑖(𝐘)‖₂ for each 𝑛 for the test are presented. Here𝐗=(𝑋(1),𝑋(2),𝑋(3))is a computed solution to (2) via (8), and𝐘=(𝑌(1),𝑌(2),𝑌(3))is a computed solution to (11) via (13) and (14). In addition, the time of execution for each method in cases𝑛=15and𝑛=20is reported. Results from experiments are given in Table 1.

Table 1 

Results from 100 runs for each value of 𝑛. 𝑍(1)=𝑍(2)=𝑍(3)=0.6𝐼_𝑛.

4. Conclusion

We have made numerical experiments for computing this solution, and we have compared the numerical results. Our numerical experiments confirm the effectiveness of the proposed new transformations which lead us to the equivalent semidefinite programming problem. We have compared the results from the experiments in regard to number of iterations and CPU time for executing the above optimization problems for𝑛=15, 𝑛=20. The solution of the optimization problems achieves the same accuracy for different number of iterations. The executed examples have demonstrated that the LMI problem performances (8) and (13) require the same average numbers of iterations (see the corresponding columns “av It” for all tests). In addition, the LMI performance (14) requires less number of iterations than the remaining approaches, and it is faster than the others—see the CPU time for execution from Table 1. This property of (14) follows from the structure of problem (14)—there is one inequality less than (13). However, this property is not valid in all cases.

Acknowledgment

This research was financially supported in part by the Sofia University “St. Kl. Ohridski” under a project of the 2013 year.