# An Improved Iterative Method for Solving the Discrete Algebraic Riccati Equation

Authors: Li Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The discrete algebraic Riccati equation has wide applications, especially in networked systems and optimal control systems. In this paper, according to the damped Newton method, two iterative algorithms with a stepsize parameter is proposed to solve the discrete algebraic Riccati equation, one of which is an extension of Algorithm (4.1) in Dai and Bai (2011). A numerical example demonstrates the convergence effect of the presented algorithm.


---

## Body

## 1. Introduction and Preliminaries

The discrete algebraic Riccati equation plays an important part in engineering, such as optimal control systems [1], modified filtering [2, 3], and networked systems[4–7]. Consider the following discrete-time linear system:(1)xk+1=Axk+Buk,where xk∈Rn is the state variable, uk∈Rr is the input variable, B∈Rn×r is the input matrix, and A∈Rn×n is the system matrix and is always invertible [8]. The optimal state feedback controller of (1) is(2)u∗k=−G+BTPB−1BTPAxk,which minimizes the quadratic performance index of (1) and is closely related to the discrete algebraic Riccati equation (DARE):(3)P=ATPA−ATPBG+BTPB−1BTPA+Q,where Q∈Rn×n is semipositive definite, G∈Rr×r is positive definite, and P∈Rn×n is the positive definite solution of the DARE (3). Let R=BG−1BT≥0. According to the matrix identity,(4)X−1+YZ−1=X−XYI+ZXY−1ZX,equation (3) can be transformed into(5)P=ATP−1+R−1A+Q.Due to the wide applications of the DARE, many works have been proposed to discuss the DARE. Various bounds and solutions about the DARE have been provided, such as upper and lower solution bounds [9–14], bounds about sum and product of eigenvalues [15, 16], determinant of the solution [17], and the existence of the solution [18–21]. However, in an optimal control system, we often need to compute the solution of the DARE to find the optimal state feedback controller which minimizes the quadratic performance index. It is very difficult to solve the DARE, especially when the dimensions of the coefficient matrices are high. So, many researchers provide a lot of iterative methods to solve this equation. Komaroff present a fixed-point iterative algorithm that needs to compute twice matrix inversion at each step [22]. By Newton’s method, Guo derived the maximal symmetric solution of the DARE in [23]. The structure-preserving doubling algorithms are discussed in [24–27]. The Schur method is adopted to solve algebraic Riccati equations [28]. Recently, Dai and Bai propose an iterative algorithm that partially avoids computing the matrix inversions by making use of the Schulz iteration [29].In Section2, we propose two iterative algorithms with a stepsize parameter to solve the DARE by the damped Newton method. One of the iterative algorithms is an extension of Algorithm 4.1 in [29]. Numerical example is given in Section 3 to demonstrate the convergence effect of our algorithms.We first introduce some symbol conventions.R denotes the real number field. Rn×m denotes the set of n×m real matrices. For X=xij∈Rn×n, let XT, X−1, X, and λminX denote the transpose, inverse, spectral norm, and the minimal eigenvalue of the matrix X, respectively. The inequality X>≥0 means X is a symmetric positive (semi-) definite matrix; and the inequality X>≥Y means X−Y is a symmetric positive (semi-) definite matrix. The identity matrix with appropriate dimensions is represented by I.Lemma 1.
(see [30]). If A,B∈Rn×n are symmetric positive definite matrices, then(6)A≥B,if and only ifB−1≥A−1.Lemma 2.
(see [31]). Let C and P be Hermitian matrices of the same order and let P>0. Then,(7)CPC+P−1≥2C.Lemma 3.
(see [32]). Let SandT be symmetric positive definite matrices. Then,(8)S≥T,ifS≥T≥0.


## 2. Improved Iterative Algorithms for Solving the DARE

To find the positive definite solution of the DARE (5), Dai and Bai, in [29], proposed an algorithm that partially avoids computing the matrix inversion as follows.Algorithm 1.
(see [29]). Take Y0=Q−1+R−1. For k=0,1,2,⋯ , compute(9)Pk+1=ATYkA+Q,Yk+1=Yk2I−Pk+1−1+RYk.
In this section, we propose two iterative algorithms to solve the DARE (5), which are motivated by the damped Newton method [33] and the methods in [34, 35]. Let us recall the damped Newton method to find the root of FZ=Z−1−B:(10)Zk+1=Zk+tZkI−BZk=1+tZk−tZkBZk,where t>0 is a stepsize parameter. If the initial matrix is near the solution of the problem, the unit stepsize t=1 can be accepted in the local Newton method. However, it is not suitable to choose t=1 if the initial matrix is far from the solution of the problem [33].
The DARE (5) can be translated into FP=0, where(11)FP=A−TP−QA−1−1−P−1+R.
LetZ=A−TP−QA−1 and B=P−1+R. Then, to find the root of FP is equivalent to find the root of FZ, we can solve the DARE (5) by constructing an iterative scheme. According to (10), we present the following iterative algorithms for the DARE (5).Algorithm 2.
Step 1: setP0=Q, Y0=Q−1+R−1 and t>0.Step 2: compute(12)Yk+1=1−tYk+t2Yk−YkPk−1+RYk,Pk+1=ATYk+1A+Q,k=0,1,2,⋯.Algorithm 3.
Step 1: setY0=Q−1+R−1, P0=Q and t>0.Step 2: compute(13)Pk+1=ATYkA+Q,Yk+1=1−tYk+t2Yk−YkPk−1+RYk,k=0,1,2,⋯.
About Algorithms2 and 3, we have the following results.Theorem 1.
LetP− be the positive definite solution of the DARE (5) and Q>0. The iterative sequences Pk and Yk are generated by Algorithm 2 with t∈0,1; then,(14)P0≤P1≤P2≤⋯,limk⟶∞Pk=P−;Y0≤Y1≤Y2≤⋯,limk⟶∞Yk=P−−1+R−1.Proof.
We first provePk and Yk are monotone increasing by induction. Since P− is positive definite solution of DARE (5), then(15)P−=ATP−−1+R−1A+Q.
Thus,P−≥Q.(i)
Since(16)Y1=1−tY0+t2Y0−Y0P0−1+RY0=1−tY0+t2Y0−Y0Q−1+RY0=1−tY0+t2Y0−Y0=Y0,then by Lemma1, we obtain(17)P1=ATY1A+Q=ATY0A+Q=ATQ−1+R−1A+Q,(18)≤ATP−−1+R−1A+Q=P−.From (17), we also obtain P1≥Q=P0; then,(19)P−≥P1≥P0,(20)Y0=Y1=Q−1+R−1≤P1−1+R−1≤P−−1+R−1.By Lemma2 and (20), we have(21)Y2=1−tY1+t2Y1−Y1P1−1+RY1≤1−tY1+tP1−1+R−1≤1−tP1−1+R−1+tP1−1+R−1=P1−1+R−1≤P−−1+R−1.By (20) and Lemma 1, we obtain(22)Y2=1−tY1+t2Y1−Y1P1−1+RY1≥1−tY1+t2Y1−Y1Q−1+RY1=1−tY1+t2Y1−Y1=Y1,thereby,(23)P2=ATY2A+Q≥ATY1A+Q=P1.By (21), we obtain(24)P2=ATY2A+Q≤ATP1−1+R−1A+Q≤ATP−−1+R−1A+Q=P−.Thus, from the above-mentioned proof, we have(25)P0≤P1≤P2≤P−,Y0=Y1≤Y2≤P1−1+R−1≤P−−1+R−1.(ii)
Assume that(26)Pi−1≤Pi≤⋯≤P−,Yi−1≤Yi≤Pi−1−1+R−1≤P−−1+R−1,i=1,2,⋯,k.
From (26), we get Yk≤Pk−1−1+R−1≤Pk−1+R−1; then,(27)Yk−1≥Pk−1+R.
Thus,(28)Yk+1=1−tYk+t2Yk−YkPk−1+RYk≥1−tYk+t2Yk−YkYk−1Yk≥1−tYk+tYk=Yk,(29)Yk+1=1−tYk+t2Yk−YkPk−1+RYk≤1−tYk+tPk−1+R−1≤1−tPk−1+R−1+tPk−1+R−1≤Pk−1+R−1≤P−−1+R−1.
By (28) and (29), we have(30)Pk+1=ATYk+1A+Q≥ATYkA+Q=Pk,(31)Pk+1=ATYk+1A+Q≤ATP−−1+R−1A+Q=P−.
So, we obtain(32)Pk≤Pk+1≤P−,Yk≤Yk+1≤Pk−1+R−1≤P−−1+R−1,k=0,1,2,⋯
Thus, the proof of induction is completed. Moreover, asPk and Yk are monotone increasing and they are bounded, then limk⟶∞Pk and limk⟶∞Yk exist. Taking limits in Algorithm 2 gives limk⟶∞Yk=P−−1+R−1 and limk⟶∞Pk=P−.Theorem 2.
LetP− be the positive definite solution of the DARE (5). After k steps of iteration for Algorithm 2, we have I−YkPk−1+R<ε; then,(33)ATPk−1+R−1A+Q−Pk≤εP−−Q.Proof.
According to (15), we have(34)ATP−−1+R−1A=P−−Q.
Then, by Algorithm2, Lemma 3, and (34), we obtain(35)ATPk−1+R−1A+Q−Pk=ATPk−1+R−1A−ATYkA=ATPk−1+R−1−YkA=ATI−YkPk−1+RPk−1+R−1A≤I−YkPk−1+R⋅ATPk−1+R−1A≤I−YkPk−1+R⋅ATP−−1+R−1A=I−YkPk−1+R⋅P−−Q≤εP−−Qbecause of Pk≤P−.
As the proof method is similar to Theorem1, we list the monotonicity and convergence of Algorithm 3 without proof.Theorem 3.
LetP− be the positive definite solution of the DARE (5) and Q>0. The iterative sequences Pk and Yk are generated by Algorithm 3 with t∈0,1 and start from Y0=Q−1+R−1 and P0=Q; then, Pk is monotone increasing and converges to P−, and Yk is monotone increasing and converges to P−−1+R−1.Remark 1.
For Algorithms2 and 3, we find the steps of iteration for Algorithm 2 are less than Algorithm 3 and the convergence speed of the Algorithm 2 is faster than Algorithm 3 from the numerical examples. Therefore, in the following example, we only discuss the superiority and effectiveness of Algorithm 2.


## 3. Numerical Examples

In this section, we present the following numerical example to show the effectiveness of our results. We also discuss the performance of Algorithm2 with different t values. The whole process is carried out on Matlab 7.1 and the precision is 10−8.Example 1.
Consider the discrete system (1) with(36)A=2.270.130.120.1−0.132.340.120.050.11−0.171.90.030.010.070.021.1,B=1.1500.01000.80000.040.900.02001.8,Q=0.1200.1002.2000.101.400000.7,G=I.
In [29], Dai and Bai choose the starting matrix Y0=Q−1+R−1. After 17 steps of iteration, the required precision is derived, and the residual ATP−1+R−1A+Q−P is 2.0754e−009.
For Algorithm2, we choose P0=Q,Y0=Q−1+R−1 and give the steps of iteration and the residual as Table 1 with a different parameter t when the process is stopped under the required precision. When t is near 1, we find that the steps of iteration are less than [29]. Especially, when t=1.2, it only needs 10 steps for Algorithm 2 to converge to the iterative solution:(37)P10=3.3299−0.31200.52020.1433−0.31209.6394−0.12920.19040.5202−0.12924.97310.08200.14330.19040.08200.9962with the residual ATP−1+R−1A+Q−P=5.6438e−009, and Algorithm 2 has faster convergence speed than Algorithm 1 from Figure 1. Moreover, from Table 1, we see that Algorithm 2 is more efficient when t>1. Although we only prove the convergence of Algorithm 2 when t∈0,1, in this paper, Algorithm 2 works well in practical computation when t>1.Table 1
Numerical results.
tIterationsResidualtIterationsResidual0.6359.0186e – 0091.2105.6438e − 0090.8245.4562e − 0091.3129.2094e − 0090.9204.7874e − 0091.5189.4721e − 0091172.0754e − 0091.8408.6619e − 0091.1141.0160e − 0092.01049.1141e − 009Figure 1
The relationship between iteration step and residual.



---

*Source: 3283157-2020-05-20.xml*