 # Numerical Method for a Class of Algebraic Riccati Equations

We study an iteration approach to solve the coupled algebraic Riccati equations when they appear in general two player closed-loop type Nash differential games over an infinite time horizon. Also, we propose an effective algorithm for finding positive definite solutions. In particular, we present various numerical examples connected with matrix Riccati equations according to different dimensions.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Before we present our research, it is worth pointing out the real life experience (community volunteer activity) relevant to our research work. Besides the interest in our own problem, this class of equations appears, for instance, in the optimal control problem for infinite Markov jump linear systems  or linear-quadratic two-person zero-sum differential games over an infinite horizon 

. An important recent development in the use of invariant subspace-based methods is the use of iterative refinement. As we know, this idea is closely related to the solution of the Riccati equation in Newton’s method. In particular, Newton’s method itself is an interesting yet useful technique, and now a fairly complete theory has been developed for this algorithm, including a computable estimate of the convergence area. We refer the interested reader to

[9, 6, 7, 4, 5] for details.

The algebraic Riccati equations play fundamental roles in engineering, management, economic, finance, linear-quadratic control and estimation systems as well as in many branches of applied mathematics. The purpose of this article is not to investigate the vast literature available for these equations, but to provide readers with references [9, 6, 7, 4, 5, 12, 8, 2, 10]

. Generally, algebraic Riccati equations that appear in game theory are non-symmetric, so this makes the effective numerical methods developed for symmetric equations based on control theory inapplicable. Our problem is to find the solution of the coupled algebraic Riccati equations generated by two-person linear-quadratic non-zero sum deterministic differential games over an infinite horizon with closed-loop Nash equilibria. These equations based on closed-loop type Nash equilibria are symmetric, but they are coupled together so that the methods established to solve control problems cannot be directly applied. In this article, we develop an iteration approach to tackle these coupled algebraic Riccati equations, and propose an effective algorithm for finding positive definite solutions.

The article is organized as follows. Section 2 introduces some basic notions and presents formulation and algorithm. In Section 3, we present several numerical examples according to different dimensions to illustrate the proposed approach. Section 4 concludes this research work.

## 2 Formulation and Algorithm

Throughout this paper, we let and be the set of all (real) matrices and symmetric (real) matrices. To simplify the notation, we denote (for )

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩A∈Rn×n,B=(B1,B2),B1∈Rn×m1,B2∈Rn×m2,Q1,Q2∈Sn,Ri=(Ri11Ri12Ri21Ri22),Ri11∈Sm1,Ri22∈Sm2,Ri12=R⊤i21∈Rm1×m2,

where the notation denotes transpose of matrix .

###### Definition 2.1.

The following is called a system of coupled algebraic Riccati equations (AREs):

 {X1A+A⊤X1+Θ⊤R1Θ+X1BΘ+Θ⊤B⊤X1+Q1=0,X2A+A⊤X2+Θ⊤R2Θ+X2BΘ+Θ⊤B⊤X2+Q2=0, (2.1)

where , and satisfies

 (B⊤1X1B⊤2X2)+(R111R112R221R222)(Θ1Θ2)=0, (2.2)

and

 R111≻0,R222≻0.

Problem. For , and , find a pair of solution to satisfy equations (2.1)-(2.2), where and .

Since we cannot solve the coupled algebraic Riccati equations (2.1)-(2.2) directly, we separate the coupled equations into two decoupled equations (2.4) and (2.7) for finding positive definite solutions using iterations in the following algorithm.

[t]0.9 Algorithm:

• Let be initial matrix solution of iteration and .

• Substituting into (2.2) yields .

• Set

 A1Δ=A+B2Θ∗2,Q1Δ=Q1+(Θ∗2)⊤R122Θ∗2,S11Δ=R112Θ∗2. (2.3)

Then the first equation in (2.1) and the first equation in (2.2) can be reduced to

 X1A1+A⊤1X1+Q1−(X1B1+S⊤11)R−1111(B⊤1X1+S11)=0, (2.4)

and

 B⊤1X1+S11+R111Θ1=0, (2.5)

respectively. Solving (2.4) yields its solution , and substituting into (2.5) leads to .

[t]0.9 Algorithm:

• Similarly, set

 A2Δ=A+B1Θ∗1,Q2Δ=Q2+(Θ∗1)⊤R211Θ∗1,S22Δ=R221Θ∗1. (2.6)

Then the second equation in (2.1) and the second equation in (2.2) can be reduced to

 X2A2+A⊤2X2+Q2−(X2B2+S⊤22)R−1222(B⊤2X2+S22)=0, (2.7)

and

 B⊤2X2+S22+R222Θ2=0, (2.8)

respectively. Solving (2.7) yields its solution , and substituting into (2.8) leads to .

• If , then set and go to Step 2. Otherwise, stop.

Since and for , we have . In addition, and . Therefore, applying the Schur method in , we can tackle equations (2.4) and (2.7

) directly. Although the Schur-based method has all the reliability and efficiency, it is relatively difficult to effectively parallelize and vectorize. Therefore, other methods have been re-examined in order to implement them on various advanced computing architectures. Due to its non-linearity in

, it is difficult to study its numerical analysis. Numerical theory will be the subject of our future research. In addition, we also can tackle equations (2.4) and (2.7) using a computational approach via a semi-definite programming associated with linear matrix inequalities, whose feasibility is equivalent to the solvability of equations (2.4) and (2.7).

## 3 Examples

In this section, we present some examples illustrating the algorithm in the previous section. We discuss about five cases: (1) , and ; (2) , and ; (3) , and ; (4) , and ; (5) , and . All computations are done on an iMac using Matlab software and double precision arithmetic.

###### Example 3.1.

Consider the following problem with , and . In this example,

 ⎧⎪⎨⎪⎩A=1,B1=B2=1,Q1=2,Q2=2.5,R1=(1000),R2=(0002).

Then the corresponding system of generalized AREs (2.1)-(2.2) reads

 {X1+X1+Θ⊤R1Θ+X1(1,1)Θ+Θ⊤(1,1)⊤X1+2=0,X1+Θ1=0,X2+X2+Θ⊤R2Θ+X2(1,1)Θ+Θ⊤(1,1)⊤X2+2.5=0,X2+2Θ2=0, (3.1)

which is equivalent to

 {2X1−X21−X1X2+2=0,2X2−12X22−2X1X2+2.5=0. (3.2)

This is an example with two scalar variables and two quadratic equations. Let be initial solution of iteration and . Using algorithm in Section 2, we can get

 {X1A1+A1X1+2−X21=0,X1+Θ1=0,X2A2+A2X2+2.5−12X22=0,X2+2Θ2=0.

Several iterations leads to a pair of solution satisfying (3.1) and its equivalent equations (3.2).

###### Example 3.2.

Consider the following problem with , and . In this example,

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩A=(−100−1),B=(B1,B2),B1=(10),B2=(01),Q1=(3002.75),Q2=(2002.5),R1=(1001),R2=(2002).

Then the corresponding system of generalized AREs (2.1)-(2.2) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩X1A+A⊤X1+Θ⊤R1Θ+X1BΘ+Θ⊤B⊤X1+Q1=0,X2A+A⊤X2+Θ⊤R2Θ+X2BΘ+Θ⊤B⊤X2+Q2=0,(B⊤1X1B⊤2X2)+(1002)(Θ1Θ2)=0, (3.3)

where , and . Let be initial matrix solution of iteration, where

 ˆX1=(0.5000.5),ˆX2=(0.5000.5), (3.4)

and . Using algorithm in Section 2, we can get

 {A1=A+B2Θ∗2,Q1=Q1,S11=0,A2=A+B1Θ∗1,Q2=Q2,S22=0.

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A1+A⊤1X1+Q1−X1B1(B⊤1X1)=0,X2A2+A⊤2X2+Q2−12X2B2(B⊤2X2)=0,(B⊤1X1B⊤2X2)+(1002)(Θ1Θ2)=0.

Several iterations leads to a pair of solution

 X∗1=(1001)andX∗2=(1001) (3.5)

satisfying (3.3).

###### Example 3.3.

Consider the following problem with , and . In this example,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩A=(−100−1),B=(B1,B2),B1=(1001),B2=(11),Q1=(4.252.252.253),Q2=(2.5002),R1=⎛⎜⎝10.500.510001⎞⎟⎠,R2=⎛⎜⎝200021012⎞⎟⎠.

Then the corresponding system of generalized AREs (2.1)-(2.2) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A+A⊤X1+Θ⊤R1Θ+X1BΘ+Θ⊤B⊤X1+Q1=0,X2A+A⊤X2+Θ⊤R2Θ+X2BΘ+Θ⊤B⊤X2+Q2=0,(B⊤1X1B⊤2X2)+⎛⎜⎝10.500.510012⎞⎟⎠(Θ1Θ2)=0, (3.6)

where , and . Let be initial matrix solution of iteration, where

 ˆX1=(0.5000.5),ˆX2=(0.5000.5), (3.7)

and . Using algorithm in Section 2, we can get

 {A1=A+B2Θ∗2,Q1=Q1+(Θ∗2)⊤Θ∗2,S11=0,A2=A+B1Θ∗1,Q2=Q2+(Θ∗1)⊤R211Θ∗1,S22=R221Θ∗1,

where and . Then (2.4)-(2.5) and (2.7)-(2.8) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A1+A⊤1X1+Q1−X1B1R−1111(B⊤1X1)=0,X2A2+A⊤2X2+Q2−12(X2B2+S⊤22)(B⊤2X2+S22)=0,(B⊤1X1B⊤2X2)+⎛⎜⎝10.500.510012⎞⎟⎠(Θ1Θ2)=0,

where . Several iterations leads to a pair of solution

 X∗1=(10.50.51)andX∗2=(1001) (3.8)

satisfying (3.6).

###### Example 3.4.

Consider the following problem with , and . In this example,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩A=(−100−1),B=(B1,B2),B1=(1001),B2=(10.50.51),Q1=⎛⎝412234234316⎞⎠,Q2=(8−12−12316),R1=⎛⎜ ⎜ ⎜⎝10.5000.51000010.5000.51⎞⎟ ⎟ ⎟⎠,R2=⎛⎜ ⎜ ⎜⎝2000021001210012⎞⎟ ⎟ ⎟⎠.

Then the corresponding system of generalized AREs (2.1)-(2.2) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A+A⊤X1+Θ⊤R1Θ+X1BΘ+Θ⊤B⊤X1+Q1=0,X2A+A⊤X2+Θ⊤R2Θ+X2BΘ+Θ⊤B⊤X2+Q2=0,(B⊤1X1B⊤2X2)+⎛⎜ ⎜ ⎜⎝10.5000.510001210012⎞⎟ ⎟ ⎟⎠(Θ1Θ2)=0, (3.9)

where , and . Let be initial matrix solution of iteration, where

 ˆX1=(0.5000.5),ˆX2=(0.5000.5), (3.10)

and . Using algorithm in Section 2, we can get

 {A1=A+B2Θ∗2,Q1=Q1+(Θ∗2)⊤R122Θ∗2,S11=0,A2=A+B1Θ∗1,Q2=Q2+(Θ∗1)⊤R211Θ∗1,S22=R221Θ∗1,

where , and . Then (2.4)-(2.5) and (2.7)-(2.8) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A1+A⊤1X1+Q1−X1B1R−1111(B⊤1X1)=0,X2A2+A⊤2X2+Q2−(X2B2+S⊤22)R−1222(B⊤2X2+S22)=0,(B⊤1X1B⊤2X2)+⎛⎜ ⎜ ⎜⎝10.5000.510001210012⎞⎟ ⎟ ⎟⎠(Θ1Θ2)=0,

where and . Several iterations leads to a pair of solution

 X∗1=(10.50.51)andX∗2=(2001) (3.11)

satisfying (3.9).

###### Example 3.5.

Consider the following problem with , and . In this example,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩A=⎛⎜⎝−50−10−100−10−5⎞⎟⎠,B=(B1,B2),B1=⎛⎜⎝100101⎞⎟⎠,B2=⎛⎜⎝10.50.5101⎞⎟⎠,Q1=⎛⎜⎝13.185219.11113.870419.111143.41673.97223.87043.972211.8241⎞⎟⎠,Q2=⎛⎜⎝21.8519−0.72223.9815−0.722226.166725.22223.981525.222231.6852⎞⎟⎠,R1=⎛⎜ ⎜ ⎜⎝10.5000.51000010.5000.51⎞⎟ ⎟ ⎟⎠,R2=⎛⎜ ⎜ ⎜⎝2000021001210012⎞⎟ ⎟ ⎟⎠.

Then the corresponding system of generalized AREs (2.1)-(2.2) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A+A⊤X1+Θ⊤R1Θ+X1BΘ+Θ⊤B⊤X1+Q1=0,X2A+A⊤X2+Θ⊤R2Θ+X2BΘ+Θ⊤B⊤X2+Q2=0,(B⊤1X1B⊤2X2)+⎛⎜ ⎜ ⎜⎝10.5000.510001210012⎞⎟ ⎟ ⎟⎠(Θ1Θ2)=0, (3.12)

where , and . Let be initial matrix solution of iteration, where

 ˆX1=⎛⎜⎝0.50000.50000.5⎞⎟⎠,ˆX2=⎛⎜⎝0.50000.50000.5⎞⎟⎠, (3.13)

and . Using algorithm in Section 2, we can get

 {A1=A+B2Θ∗2,Q1=Q1+(Θ∗2)⊤R122Θ∗2,S11=0,A2=A+B1Θ∗1,Q2=Q2+(Θ∗1)⊤R211Θ∗1,S22=R221Θ∗1,

where , and . Then (2.4)-(2.5) and (2.7)-(2.8) reads

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩X1A1+A⊤1X1+Q1−X1B1R−1111(B⊤1X1)=0,X2A2+A⊤2X2+Q2−(X2B2+S⊤22)R−1222(B⊤2X2+S22)=0,(B⊤1X1B⊤2X2)+⎛⎜ ⎜ ⎜⎝10.5000.510001210012⎞⎟ ⎟ ⎟⎠(Θ1Θ2)=0,

where and . Several iterations leads to a pair of solution

 X∗1=⎛⎜⎝110120001⎞⎟⎠andX∗2=⎛⎜⎝200011012⎞⎟⎠ (3.14)

satisfying (3.12).

## 4 Conclusion

We have learned many methods of algebraic Riccati equations and their performance in computer finite arithmetic environments. However, many important questions remain unanswered and are the subject of ongoing research topics.

## References

•  M. Ait Rami and X.Y. Zhou, Linear matrix inequalities, Riccati equations, and indefinite stochastic linear quadratic controls, IEEE Transactions on Automatic Control, 45 (2000), 1131-1143.
•  T.P. Azevedo-Perdicoúlis and G. Jank, Iterative solution of algebraic matrix Riccati equations in open loop Nash games, International Journal of Robust Nonlinear Control, 15 (2005), 55-62.
•  J. Baczynski and M.D. Fragoso, Maximal solution to algebraic Riccati equations linked to infinite Markov jump linear systems, Mathematics of Control, Signals, and Systems, 20 (2008), 157-172.
•  P. Benner and T. Stykel, Numerical solution of projected algebraic Riccati equations, SIAM Journal on Numerical Analysis, 52 (2014), 581-600.
•  A.S.A. Dilip and H.K. Pillai, Characterization of solutions of non-symmetric algebraic Riccati equations, Linear Algebra and its Applications, 507 (2016), 356-372.
•  C.H. Guo and N.J. Higham, Iterative solution of a non-symmetric algebraic Riccati equation, SIAM Journal on Matrix Analysis and Applications, 29 (2007), 396-412.
•  C.H. Guo and A.J. Laub, On the iterative solution of a class of non-symmetric algebraic Riccati equations, SIAM Journal on Matrix Analysis and Applications, 22 (2000), 376-391.
•  D. Kremer and R. Stefan, Non-symmetric Riccati theory and linear quadratic Nash games, Fifteenth International Symposium on Mathematical Theory of Networks and Systems, University of Notre Dame, August 12-16, 2002.
•  A.J. Laub, A Schur method for solving algebraic Riccati equation, IEEE Transactions on Automatic Control, 24 (1979), 913–921.
•  X. Li, J. Shi and J. Yong, Mean-field linear-quadratic stochastic differential games in an infinite horizon, arXiv:2007.06130, 2020.
•  J. Sun, J. Yong and S. Zhang, Linear quadratic stochastic two-person zero-sum differential games in an infinite horizon, ESAIM COCV, 22 (2016), 743–769.
•  J. Xue and R.C. Li, Highly accurate doubling algorithms for -matrix algebraic Riccati equations, Numerische Mathematik, 135 (2017), 733-767.