# The extended Hamiltonian algorithm for the solution of the algebraic Riccati equation.

1. IntroductionThe algebraic Riccati equations (AREs) have been widely used in control system syntheses [1, 2], especially in optimal control [3], robust control [4], signal processing [5], and the LMI-based design [6]. And, in the past several decades, there has been an increasing interest in the solution problems of the AREs [7-9]. The lower matrix bounds for the AREs were discussed in [10, 11]. The structure-preserving doubling algorithm (SDA) was given under assumptions which are weaker than stabilizability and detectability, as well as practical issues involved in the application of the SDA to continuous-time AREs [12]. It is the purpose of this paper to investigate the unique (under suitable hypotheses) symmetric positive definite solution of an algebraic Riccati equation. Some effective methods, such as the Schur method (SM) [13] and the new subspace iteration method (NSIM) [14], were given for the numerical solution of the algebraic Riccati equation. Recently, the Euclidean gradient algorithm (EGA) and the Riemannian gradient algorithm (RGA) for the numerical solution of the algebraic Riccati equation with a cost function of the Riemannian distance on the curved Riemannian manifold were proposed in [15].

It is well known that the Euclidean gradient algorithm (EGA) and the Riemannian gradient algorithm (RGA) are first-order learning algorithms; hence the convergence speeds of the EGA and RGA are very slow. The inclusion of a momentum term had been found to increase the rate of convergence dramatically [16, 17]. Based on this, the extended Hamiltonian algorithm (EHA) on manifolds was developed in [18], while its numerical implementation was discussed in [19]. Furthermore, the EHA is a second-order learning algorithm and converges faster than the first-order learning algorithm for optimization problem if certain conditions are satisfied.

In this paper, we will apply the EHA to give the numerical solution of the AREs. Furthermore, we give two simulations to show the efficiency of our algorithm.

Briefly, the rest of the paper is organized as follows. Section 2 is concluded with some fundamental knowledge on manifold that will be used throughout the paper. Section 3 introduces the extended Hamiltonian algorithm on manifold of positive definite symmetric matrices and Section 4 illustrates the convergence speeds of the EHA compared with other algorithms using two numerical examples. Section 5 concludes the paper.

2. Preliminaries

In this section we recall some differential geometric facts of the space of positive definite symmetric matrices that will be used in the present analysis. More details can be found in [20-23].

2.1. The Frobenius Inner Product and the Euclidean Distance. Let M(n) be the set of nx n real matrices and GL (n) be its subset containing all the nonsingular matrices. GL(n) is a Lie group, that is, a group which is also a differentiable manifold and for which the operations of group multiplication and inverse are smooth. The tangent space at the identity is called the corresponding Lie algebra denoted by gl(n). It is the space of all linear transformations in R", that is, M(n).

On M(n), we use the Euclidean inner product, known as the Frobenius inner product and defined by

[<P, Q>.sub.F] = tr ([P.sup.T]Q), (1)

where tr stands for the trace of the matrix and superscript T denotes the transpose. The associated norm [[parallel]P[parallel].sub.F] = [[tr([P.sup.T] P)].sup.1/2] is used to define the Euclidean distance on M(n) by

[d.sub.F] (P, Q) = [[parallel]P - Q[parallel].sub.F]. (2)

2.2. Geometry of Positive Definite Symmetric Matrices. We denote by

Sym (n) = {s [member of] M (n) | s = [s.sup.T]} (3)

the vector space of symmetric matrices and denote by

PD (n) = {x [member of] Sym (n) | x > 0} (4)

the set of all n x n positive definite symmetric matrices, which is a differentiable manifold of dimension n(n + 1)/2. Here x > 0 means the quadratic form [z.sup.T]xz > 0 for all z [member of] [R.sup.n] {0}. Furthermore, PD(n) is an open convex cone; namely, if [x.sub.1] and [x.sub.2] are in PD(n), so is [x.sub.1] + t[x.sub.2] for any t > 0.

The exponential of any symmetric matrix is a positive definite symmetric matrix and the (principal) logarithm of any positive definite symmetric matrix is a symmetric matrix. Thus the exponential map from Sym(n) to PD(n) is one-to-one and onto. It follows that Sym(n) provides a parameterization of PD(n) via the exponential map.

2.3. Metric and Distances on PD(n). Let g denote the Riemannian metric on manifold PD(n). For any x [member of] PD(n) the tangent space [T.sub.x]PD(n) is identified with Sym(n). On [T.sub.x]PD(n) we define the inner product and the corresponding norm as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

which depend on the point x. The positive definiteness of this metric is a consequence of the positive definiteness of the Frobenius inner product for

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (6)

Compared to manifold PD(n) with the Frobenius inner product (1) resulting in a flat Riemannian manifold, manifold PD(n) with the Riemannian metric (5) becomes a curved Riemannian manifold.

To measure closeness of two positive definite symmetric matrices [x.sub.1] and [x.sub.2] one can use the Euclidean distance 2) of the ambient space M(n); that is,

[d.sub.F] ([x.sub.1], [x.sub.1]) = [[parallel][x.sub.1] - [x.sub.2][parallel].sub.F]. (7)

It will be more appropriate to use the Riemannian distance induced from (5), which is intrinsic to PD(n) and defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (8)

where [[lambda].sub.i], i = 1, 2, ..., n, are the positive eigenvalues of [x.sup.-1.sub.1][x.sub.2]. The positivity of the [[lambda].sub.i] is a consequence of the similarity between the (in general nonsymmetric) matrix [x.sup.-1.sub.1] [x.sub.2] and the positive definite symmetric matrix [x.sup.-1/2.sub.1] [x.sub.2] [x.sup.-1/2.sub.1]. It is straightforward to see that the Riemannian distance (8) is invariant under inversion, that is,

[d.sub.R] ([x.sup.-1.sub.1], [x.sup.-1.sub.2]) = [d.sub.R] ([x.sub.1], [x.sub.2]) (9)

and is also invariant under congruent transformations; that is,

[d.sub.R] ([x.sub.1], [x.sub.2]) = [d.sub.R] ([P.sup.T] [x.sub.1] P, [P.sup.T][x.sub.2]P), (10)

where P [member of] GL(n).

Now, we discuss the geodesics on manifold PD(n) under the different metrics. The geodesic passing through x in the direction of v [member of] Sym(n) with respect to the Euclidean metric (2) is given by

[gamma](t) = x + tv, 0 [less than or equal to] t < [delta] (11)

for some positive number [delta] > 0. The restriction that t is not too large is necessary in order to ensure that the line stays within PD(n). Similarly, the geodesic passing through x in the direction of v [member of] Sym(n) with respect to the Riemannian metric (5) is expressed by

[gamma] (t) = [x.sup.1/2] exp ([x.sup.-1/2] v[x.sup.-1/2]t) [x.sup.1/2], (12)

which does not suffer the previous restriction. Let the exponential map Exp: Sym(n) [right arrow] PD(n); then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)

Notice that the geodesic (12) is defined for all values of t, which is called geodesic completeness. The Hopf-Rinow theorem states that a geodesic complete manifold is also a complete metric space. And, manifold PD(n) with the Riemannian metric (5) has nonpositive sectional curvature.

3. The Extended Hamiltonian Algorithm for the Solution of AREs

3.1. Problem Formulation. The infinite time state regulator, namely, the state equation of linear time-invariant system, is in the form

[??](t) = Ay (t) + Bw (t), y ([t.sub.0]) = [y.sub.0], (14)

where y(t) [member of] [R.sup.n] and w(t) [member of] [R.sup.m] are the state and the control vector, respectively, A and B are constant matrices of appropriate dimensions, and the terminal time [t.sub.f] = [infinity]. The problem would be to find an optimal control w(t) satisfying (14) while minimizing the quadratic performance index as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (15)

where the designed parameters satisfy that both Q and R are positive definite symmetric matrices. It has been proved that the optimal control [w.sup.*] (t) of the cost function (15) is uniquely given by [24]

[w.sup.*] (t) = -[R.sup.-1] [B.sup.T] xy (t) (16)

and the quadratic optimal performance index is gotten by

[J.sup.*] (y ([t.sub.0])) = 1/2 [y.sup.T] ([t.sub.0]) xy ([t.sub.0]), (17)

where x is a positive definite symmetric matrix satisfying the AREs; that is,

xA + [A.sup.T] x - xB[R.sup.-1] [B.sup.T]x + Q = 0. (18)

It is assumed that the pair (A, B) is completely controllable and that the pair (A, D) is completely observable, where [DD.sup.T] = Q. Under these assumptions, (18) is known to have a unique positive definite solution x [25]. There are, of course, other solutions to (18), but for the algorithm presented in this paper the emphasis will be on computing the positive definite one. Moreover, the optimal trajectory y(t) satisfies the following optimal asymptotically stable closedloop system:

[??](t) = (A - [BR.sup.-1] [B.sup.T] x) y (t) = (A - BK) y (t), (19)

where the state-feedback gain K is given by

K = [R.sup.-1] [B.sup.T] x. (20)

In order to solve (18), our purpose is to seek a matri x on manifold PD(??)such that the matrix [xBR.sup.-1] [B.sup.T] x - xA - [A.sup.T] x is as close as possible to the given matrix Q.

3.2. Distances and Gradients. In fact, the idea of formulating an equation on a manifold as a minimal-distance problem was suggested previously in [26]. Hence we use the two distance functions (7) and 8), and the difference between Q and [xBR.sup.-1] [B.sup.T] x - xA - [A.sup.T]x can be expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (21)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)

Let J : PD(n) [right arrow] R be a scalar function of n x n matrix. The Riemannian gradient of J at x, denoted by [[nabla].sub.x]J, is a map from PD(n) to Sym(n) that satisfies [<[[nabla].sub.x]J, s>.sub.x] = [<[[partial derivative].sub.x]J, s>.sub.F], for every s [member of] Sym(n), where [[partial derivative].sub.x]J is the Euclidean gradient of J at x. It follows from the inner product (5) that [27]

[[nabla].sub.x]J = x([[partial derivative].sub.x]J)x. (23)

Next, we will state some known facts, which are very helpful in computing [[partial derivative].sub.x][J.sub.F] and [[nabla].sub.x][J.sub.R].

Lemma 1 (see [28]). Let A, B denote two constant matrices and let X, Y denote matrix functions. One has the following properties:

(1) dtr (X) = tr(dX),

(2) d(A) = 0,

(3) d(AXB) = A(dX)B,

(4) d(X [+ or -] Y) = dX [+ or -] dY,

(5) d(XY) = (dX)Y + X(dY),

where d denotes the differential of a matrix (scalar) function.

Lemma 2 (see [28]). Let f(x) be a scalar function of n x n matrix x, if

df (x) = tr (Wdx) (24)

holds; then

[[partial derivative].sub.x]f(x) = [W.sup.T], (25)

where [[partial derivative].sub.x]f(x) denotes the gradient of f(x) and W denotes the derivative of f(x).

Theorem 3. Let [J.sub.F](x) be a scalar function of n x n matrix, if

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

where the state-feedback gain K is given by

holds; then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (27)

Proof. Since Q + xA+ [A.sup.T]x - [xBR.sup.-1][B.sup.T]x is symmetric, we have from Lemma 1 that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)

Using the above Lemma 2, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (29)

Therefore, it is shown that (27) is valid.

Theorem 4. Let [J.sub.R](x) be a scalar function of n x n matrix x, if

[J.sub.R] (x) = 1/2 tr {[log.sup.2] [[Q.sup.-1/2] ([xBR.sup.-1] [B.sup.T] x - xA - [A.sup.T] x) [Q.sup.-1/2]]} (30)

holds; then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (31)

Proof. Since [xBR.sup.-1][B.sup.T] - xA - [A.sup.T]x is symmetric, it follows immediately from Lemma 1 that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (32)

Using the above Lemma 2, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (33)

Combining (23) and (33), (31) is established. This completes the proof of Theorem 4.

3.3. The Extended Hamiltonian Algorithm on PD(n). The Euclidean gradient algorithm (EGA) and the Riemannian gradient algorithm (RGA) are first-order learning algorithms and inherit some drawbacks that are well known about gradient-based optimization like, for instance, the slowlearning phenomenon in presence of plateau in the error surface. In order to overcome this difficulty, Fiori generalized EGA and RGA and proposed the extended Hamiltonian algorithm on manifold. In particular, on manifold PD(n), the extended Hamiltonian algorithm can be expressed by [18]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (34)

where x [member of] PD(n), v denotes the instantaneous learning velocity, V : PD(n) [right arrow] R denotes a cost function, [[nabla].sub.x]V denotes the Riemannian gradient of V, and the constant [mu] > 0 denotes a viscosity term. The function [K.sub.x] : [T.sub.x]PD(n) x [T.sub.x] PD(n) [right arrow] R denotes the kinetic energy of the particle in a point x and has the corresponding expression [K.sub.x] = (1/2) tr[[([vx.sup.-1]).sup.2]].

System (34) maybe implemented by [19]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (35)

where [eta] is a small positive number known as the learning rate. The effectiveness of the algorithm is ensured if and only if

[square root of 2[[lambda].sub.max]] < [mu] < 1/[eta], (36)

where [[lambda].sub.max] denotes the largest eigenvalue of the Hessian matrix of the cost function V(x). See 19] for more details.

Since the Riemannian distance is the shortest one on manifold PD(n), we take the Riemannian distance (22) as the cost function V(x) of (18); that is,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (37)

In order to apply (35) to solve (18), we plug (31) into (34) and get the final extended Hamiltonian algorithm on manifold PD(n) as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (38)

which is implemented by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (39)

where [eta] > 0 denotes the learning rate and the constant [mu] denotes a viscosity term.

Now, based on the above discussion, we give the iterative formula of the EHA for the numerical solution of (18).

Algorithm 5 (EHA). For any x belonging to the considered manifold PD(n), the extended Hamiltonian algorithm is given by the following.

(1) Set [x.sub.0] as an initial input matrix x and [v.sub.0] as an initial direction of v, and then choose a desired tolerance [epsilon] > 0.

(2) Compute V([x.sub.k]).

(3) If V([x.sub.k]) < [epsilon] then stop.

(4) Update x and v by (39) and return to step (2).

4. Simulations

In this section, we give two computer simulations to demonstrate the effectiveness and performance of the proposed algorithm.

Example 1. First consider a second-order algebraic Riccati equation. In the present experiment, [eta] = 0.01, [mu] = 3, and [epsilon] = [10.sup.-8]. Moreover, A, B, Q, R, and D used in this simulation are as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (40)

which satisfy the unique definite solution condition.

To study the performance differences among the EHA, the EGA, the RGA, and the NSIM, these algorithms are, respectively, applied to solve (18). In this paper, the optimal step size [eta], which corresponds to the best convergence speed in each algorithm, is obtained by numerical experiments. For instance, it can be found that the iterations will reduce gradually as the step size changes from 0.01 to 0.1, while the algorithm will be divergent once the step size is bigger than 0.1. Thus this optimal step size [eta] is convenient to be obtained experimentally. Based on (36), the constant [mu] > 0 depends on the selection of initial value. That is to say, different initial values will be corresponding to different [mu]; hence the best [mu] is also obtained by the experiment.

As Figure 1 shows, the behavior of the cost function is shown. Clearly, in the early stages of learning, the EHA decreases much faster than the EGA, the RGA, and the NSIM with the same learning step size. The result shows that the EHA has the fastest convergence speed among four algorithms and needs 28 iterations to obtain the numerical solution of the ARE as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (41)

However, RGA realizes the given error through 36 iterations and the slowest one among them is the EGA with 88 steps. And Figure 2 shows that the kinetic energy function finally converges to 0.

Example 2. A third-order ARE is considered following the similar procedure as Example 1. In this experiment, [eta] = 0.01, [mu] = 2.5, and [epsilon] = [10.sup.-8], and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (42)

As Figure 3 shows, the convergence speed of the EHA is still the fastest one among four algorithms. Finally, we get the numerical solution of the ARE as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (43)

Moreover, the kinetic energy function also converges to 0 in the Figure 4.

5. Summary

In this paper, a second-order learning method called the extended Hamiltonian algorithm is recalled from literature and applied to solve the numerical solution for the algebraic Riccati equation. And, the convergence speed of the provided algorithm is compared with the EGA, the RGA, and the NSIM using two simulation examples. It is shown that the convergence speed of the extended Hamiltonian algorithm is the fastest one among these algorithms. http://dx.doi.org/10.1155/2014/693659

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This subject is supported by the National Natural Science Foundations of China (nos. 61179031 and 10932002).

References

[1] K. Zhou, J. Doyle, and K. Glover, Robust and Optimal Control, Prentice Hall, Upper Saddle River, NJ, USA, 1996.

[2] Z. Lin, "Global control of linear systems with saturating actuators," Automatica, vol. 34, no. 7, pp. 897-905, 1998.

[3] D.-Z. Zheng, "Optimization of linear-quadratic regulator systems in the presence of parameter Perturbations," IEEE Transactions on Automatic Control, vol. 31, no. 7, pp. 667-670, 1986.

[4] M. L. Ni and H. X. Wu, "A Riccati equation approach to the design of linear robust controllers," Automatica, vol. 29, no. 6, pp. 1603-1605, 1993.

[5] T. Kaneko, S. Fiori, and T. Tanaka, "Empirical arithmetic averaging over the compact Stiefel manifold," IEEE Transactions on Signal Processing, vol. 61, no. 4, pp. 883-894, 2013.

[6] H.-P. Liu, F.-C. Sun, K.-Z. He, and Z.-Q. Sun, "Simultaneous stabilization for singularly perturbed systems via iterative linear matrix inequalities," Acta Automatica Sinica, vol. 30, no. 1, pp. 1-7, 2004.

[7] V. Sima, Algorithms for Linear-Quadratic Optimization, vol. 200 of Monographs and Textbooks in Pure and Applied Mathematics, Marcel Dekker, New York, NY, USA, 1996.

[8] P. Benner, J.-R. Li, and T. Penzl, "Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems," Numerical Linear Algebra with Applications, vol. 15, no. 9, pp. 755-777, 2008.

[9] H. G. Xu and L. Z. Lu, "Properties of a quadratic matrix equation and the solution of the continuous-time algebraic Riccati equation," Linear Algebra and Its Applications, vol. 222, pp. 127-145, 1995.

[10] H. H. Choi and T.-Y. Kuc, "Lower matrix bounds for the continuous algebraic Riccati and Lyapunov matrix equations," Automatica, vol. 38, no. 8, pp. 1147-1152, 2002.

[11] R. Davies, P. Shi, and R. Wiltshire, "New lower solution bounds of the continuous algebraic Riccati matrix equation," Linear Algebra and Its Applications, vol. 427, no. 2-3, pp. 242-255, 2007.

[12] E. K.-W. Chu, H.-Y. Fan, and W.-W. Lin, "A structure-preserving doubling algorithm for continuous-time algebraic Riccati equations," Linear Algebra and Its Applications, vol. 396, pp. 55-80, 2005.

[13] A. J. Laub, "A Schur method for solving algebraic Riccati equations," IEEE Transactions on Automatic Control, vol. 24, no. 6, pp. 913-921, 1979.

[14] Y. Lin and V. Simoncini, "A new subspace iteration method for the algebraic Riccati equation," preprint, http://arxiv.org/ abs/1307.3843.

[15] X. Duan, Information geometry of matrix Lie groups and applications [Ph.D. thesis], Beijing Institute of Technology, 2013, (Chinese).

[16] D. Rumelhart and J. McClelland, Parallel Distributed Processing, MIT Press, 1986.

[17] N. Qian, "On the momentum term in gradient descent learning algorithms," Neural Networks, vol. 12, no. 1, pp. 145-151, 1999.

[18] S. Fiori, "Extended Hamiltonian learning on Riemannian manifolds: theoretical aspects," IEEE Transactions on Neural Networks, vol. 22, no. 5, pp. 687-700, 2011.

[19] S. Fiori, "Extended Hamiltonian learning on Riemannian manifolds: numerical aspects," IEEE Transactions on Neural Networks and Learning Systems, vol. 23, no. 1, pp. 7-21, 2012.

[20] M. Moakher, "A differential geometric approach to the geometric mean of symmetric positive-definite matrices," SIAM Journal on Matrix Analysis and Applications, vol. 26, no. 3, pp. 735-747, 2005.

[21] M. Moakher, "On the averaging of symmetric positive-definite tensors," Journal of Elasticity, vol. 82, no. 3, pp. 273-296, 2006.

[22] A. Schwartzman, Random ellipsoids and false discovery rates: statistics for diffusion tensor imaging data [Ph.D. thesis], Stanford University, 2006.

[23] J. Jost, Riemannian Geometry and Geometric Analysis, Universitext, Springer, Berlin, Germany, Third edition, 2002.

[24] R. E. Kalman, "Contributions to the theory of optimal control," Boletin de la Sociedad Matematica Mexicana, vol. 5, no. 2, pp. 102-119, 1960.

[25] R. Kalman, "When is a linear control system optimal?" Journal of Basic Engineering, vol. 86, no. 1, pp. 51-60, 1964.

[26] S. Fiori, "Solving minimal-distance problems over the manifold of real-symplectic matrices," SIAM Journal on Matrix Analysis and Applications, vol. 32, no. 3, pp. 938-968, 2011.

[27] R. Gregorio and P. Oliveira, "A proximal technique for computing the Karcher mean of symmetric positive definite matrices," Optimization Online, 2013.

[28] X. Zhang, Matrix Analysis and Application, Springer, New York, NY, USA, 2004.

Zhikun Luo, Huafei Sun, and Xiaomin Duan

School of Mathematics, Beijing Institute of Technology, Beijing 100081, China

Correspondence should be addressed to Huafei Sun; huafeisun@bit.edu.cn

Received 17 March 2014; Revised 27 April 2014; Accepted 8 May 2014; Published 20 May 2014

Academic Editor: Naseer Shahzad

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Luo, Zhikun; Sun, Huafei; Duan, Xiaomin |

Publication: | Journal of Applied Mathematics |

Article Type: | Report |

Date: | Jan 1, 2014 |

Words: | 3877 |

Previous Article: | A numerical convolution representation of potential for a disk surface density in 3D. |

Next Article: | Data filtering based recursive least squares algorithm for two-input single-output systems with moving average noises. |

Topics: |