# The Extended Hamiltonian Algorithm for the Solution of the Algebraic Riccati Equation

Authors: Zhikun Luo; Huafei Sun; Xiaomin Duan
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We use a second-order learning algorithm for numerically solving a class
of the algebraic Riccati equations. Specifically, the extended Hamiltonian algorithm based on manifold
of positive definite symmetric matrices is provided. Furthermore, this algorithm is compared with the
Euclidean gradient algorithm, the Riemannian gradient algorithm, and the new subspace iteration method.
Simulation examples show that the convergence speed of the extended Hamiltonian algorithm is the
fastest one among these algorithms.


---

## Body

## 1. Introduction

The 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. Section2 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

LetM
(
n
) be the set of n
×
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 g
l
(
n
). It is the space of all linear transformations in R
n, that is, M
(
n
).OnM
(
n
), we use the Euclidean inner product, known as the Frobenius inner product and defined by
(1)
〈
P
,
Q
〉
F
=
tr
⁡
(
P
T
Q
)
,
where tr
⁡ stands for the trace of the matrix and superscript T denotes the transpose. The associated norm ∥
P
∥
F
=
[
tr
⁡
(
P
T
P
)
]
1
/
2 is used to define the Euclidean distance on M
(
n
) by
(2)
d
F
(
P
,
Q
)
=
∥
P
-
Q
∥
F
.


### 2.2. Geometry of Positive Definite Symmetric Matrices

We denote by(3)
Sym
(
n
)
=
{
s
∈
M
(
n
)
∣
s
=
s
T
}
the vector space of symmetric matrices and denote by
(4)
PD
(
n
)
=
{
x
∈
Sym
(
n
)
∣
x
>
0
}
the set of all n
×
n positive definite symmetric matrices, which is a differentiable manifold of dimension n
(
n
+
1
)
/
2. Here x
>
0 means the quadratic form z
T
x
z
>
0 for all z
∈
R
n
∖
{
0
}. Furthermore, PD
(
n
) is an open convex cone; namely, if x
1 and x
2 are in PD
(
n
), so is x
1
+
t
x
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 fromSym
(
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
)

Letg denote the Riemannian metric on manifold PD
(
n
). For any x
∈
PD
(
n
) the tangent space T
x
PD
(
n
) is identified with Sym
(
n
). On T
x
PD
(
n
) we define the inner product and the corresponding norm as follows:
(5)
g
x
(
u
,
v
)
=
〈
u
,
v
〉
x
=
tr
⁡
(
x
-
1
u
x
-
1
v
)
,
∥
v
∥
x
=
[
tr
⁡
(
(
x
-
1
v
)
2
)
]
1
/
2
,
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
(6)
g
x
(
v
,
v
)
=
〈
v
,
v
〉
x
=
tr
⁡
(
x
-
1
/
2
v
x
-
1
/
2
x
-
1
/
2
v
x
-
1
/
2
)
=
〈
x
-
1
/
2
v
x
-
1
/
2
,
x
-
1
/
2
v
x
-
1
/
2
〉
.Compared to manifoldPD
(
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 matricesx
1 and x
2 one can use the Euclidean distance (2) of the ambient space M
(
n
); that is,
(7)
d
F
(
x
1
,
x
2
)
=
∥
x
1
-
x
2
∥
F
.
It will be more appropriate to use the Riemannian distance induced from (5), which is intrinsic to PD
(
n
) and defined by
(8)
d
R
(
x
1
,
x
2
)
=
∥
log
⁡
(
x
1
-
1
/
2
x
2
x
1
-
1
/
2
)
∥
F
=
(
∑
i
=
1
n
‍
ln
⁡
2
λ
i
)
1
/
2
,
where λ
i,   i
=
1,2
,
…
,
n, are the positive eigenvalues of x
1
-
1
x
2. The positivity of the λ
i is a consequence of the similarity between the (in general nonsymmetric) matrix x
1
-
1
x
2 and the positive definite symmetric matrix x
1
-
1
/
2
x
2
x
1
-
1
/
2. It is straightforward to see that the Riemannian distance (8) is invariant under inversion, that is,
(9)
d
R
(
x
1
-
1
,
x
2
-
1
)
=
d
R
(
x
1
,
x
2
)
and is also invariant under congruent transformations; that is,
(10)
d
R
(
x
1
,
x
2
)
=
d
R
(
P
T
x
1
P
,
P
T
x
2
P
)
,
where P
∈
GL
(
n
).Now, we discuss the geodesics on manifoldPD
(
n
) under the different metrics. The geodesic passing through x in the direction of v
∈
Sym
(
n
) with respect to the Euclidean metric (2) is given by
(11)
γ
(
t
)
=
x
+
t
v
,
0
≤
t
<
δ
for some positive number δ
>
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
∈
Sym
(
n
) with respect to the Riemannian metric (5) is expressed by
(12)
γ
(
t
)
=
x
1
/
2
exp
⁡
(
x
-
1
/
2
v
x
-
1
/
2
t
)
x
1
/
2
,
which does not suffer the previous restriction. Let the exponential map Exp: Sym
(
n
)
→
PD
(
n
); then
(13)
x
2
=
Ex
p
x
1
(
v
)
=
x
1
1
/
2
exp
⁡
(
x
1
-
1
/
2
v
x
1
-
1
/
2
)
x
1
1
/
2
.
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.


## 2.1. The Frobenius Inner Product and the Euclidean Distance

LetM
(
n
) be the set of n
×
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 g
l
(
n
). It is the space of all linear transformations in R
n, that is, M
(
n
).OnM
(
n
), we use the Euclidean inner product, known as the Frobenius inner product and defined by
(1)
〈
P
,
Q
〉
F
=
tr
⁡
(
P
T
Q
)
,
where tr
⁡ stands for the trace of the matrix and superscript T denotes the transpose. The associated norm ∥
P
∥
F
=
[
tr
⁡
(
P
T
P
)
]
1
/
2 is used to define the Euclidean distance on M
(
n
) by
(2)
d
F
(
P
,
Q
)
=
∥
P
-
Q
∥
F
.


## 2.2. Geometry of Positive Definite Symmetric Matrices

We denote by(3)
Sym
(
n
)
=
{
s
∈
M
(
n
)
∣
s
=
s
T
}
the vector space of symmetric matrices and denote by
(4)
PD
(
n
)
=
{
x
∈
Sym
(
n
)
∣
x
>
0
}
the set of all n
×
n positive definite symmetric matrices, which is a differentiable manifold of dimension n
(
n
+
1
)
/
2. Here x
>
0 means the quadratic form z
T
x
z
>
0 for all z
∈
R
n
∖
{
0
}. Furthermore, PD
(
n
) is an open convex cone; namely, if x
1 and x
2 are in PD
(
n
), so is x
1
+
t
x
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 fromSym
(
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
)

Letg denote the Riemannian metric on manifold PD
(
n
). For any x
∈
PD
(
n
) the tangent space T
x
PD
(
n
) is identified with Sym
(
n
). On T
x
PD
(
n
) we define the inner product and the corresponding norm as follows:
(5)
g
x
(
u
,
v
)
=
〈
u
,
v
〉
x
=
tr
⁡
(
x
-
1
u
x
-
1
v
)
,
∥
v
∥
x
=
[
tr
⁡
(
(
x
-
1
v
)
2
)
]
1
/
2
,
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
(6)
g
x
(
v
,
v
)
=
〈
v
,
v
〉
x
=
tr
⁡
(
x
-
1
/
2
v
x
-
1
/
2
x
-
1
/
2
v
x
-
1
/
2
)
=
〈
x
-
1
/
2
v
x
-
1
/
2
,
x
-
1
/
2
v
x
-
1
/
2
〉
.Compared to manifoldPD
(
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 matricesx
1 and x
2 one can use the Euclidean distance (2) of the ambient space M
(
n
); that is,
(7)
d
F
(
x
1
,
x
2
)
=
∥
x
1
-
x
2
∥
F
.
It will be more appropriate to use the Riemannian distance induced from (5), which is intrinsic to PD
(
n
) and defined by
(8)
d
R
(
x
1
,
x
2
)
=
∥
log
⁡
(
x
1
-
1
/
2
x
2
x
1
-
1
/
2
)
∥
F
=
(
∑
i
=
1
n
‍
ln
⁡
2
λ
i
)
1
/
2
,
where λ
i,   i
=
1,2
,
…
,
n, are the positive eigenvalues of x
1
-
1
x
2. The positivity of the λ
i is a consequence of the similarity between the (in general nonsymmetric) matrix x
1
-
1
x
2 and the positive definite symmetric matrix x
1
-
1
/
2
x
2
x
1
-
1
/
2. It is straightforward to see that the Riemannian distance (8) is invariant under inversion, that is,
(9)
d
R
(
x
1
-
1
,
x
2
-
1
)
=
d
R
(
x
1
,
x
2
)
and is also invariant under congruent transformations; that is,
(10)
d
R
(
x
1
,
x
2
)
=
d
R
(
P
T
x
1
P
,
P
T
x
2
P
)
,
where P
∈
GL
(
n
).Now, we discuss the geodesics on manifoldPD
(
n
) under the different metrics. The geodesic passing through x in the direction of v
∈
Sym
(
n
) with respect to the Euclidean metric (2) is given by
(11)
γ
(
t
)
=
x
+
t
v
,
0
≤
t
<
δ
for some positive number δ
>
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
∈
Sym
(
n
) with respect to the Riemannian metric (5) is expressed by
(12)
γ
(
t
)
=
x
1
/
2
exp
⁡
(
x
-
1
/
2
v
x
-
1
/
2
t
)
x
1
/
2
,
which does not suffer the previous restriction. Let the exponential map Exp: Sym
(
n
)
→
PD
(
n
); then
(13)
x
2
=
Ex
p
x
1
(
v
)
=
x
1
1
/
2
exp
⁡
(
x
1
-
1
/
2
v
x
1
-
1
/
2
)
x
1
1
/
2
.
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(14)
y
˙
(
t
)
=
A
y
(
t
)
+
B
w
(
t
)
,
y
(
t
0
)
=
y
0
,
where y
(
t
)
∈
R
n and w
(
t
)
∈
R
m are the state and the control vector, respectively, A and B are constant matrices of appropriate dimensions, and the terminal time t
f
=
∞. The problem would be to find an optimal control w
(
t
) satisfying (14) while minimizing the quadratic performance index as follows:
(15)
J
=
1
2
∫
0
∞
‍
[
y
T
(
t
)
Q
y
(
t
)
+
w
T
(
t
)
R
w
(
t
)
]
d
t
,
where the designed parameters satisfy that both Q and R are positive definite symmetric matrices. It has been proved that the optimal control w
*
(
t
) of the cost function (15) is uniquely given by [24]
(16)
w
*
(
t
)
=
-
R
-
1
B
T
x
y
(
t
)
and the quadratic optimal performance index is gotten by
(17)
J
*
(
y
(
t
0
)
)
=
1
2
y
T
(
t
0
)
x
y
(
t
0
)
,
where x is a positive definite symmetric matrix satisfying the AREs; that is,
(18)
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
+
Q
=
0
.It is assumed that the pair(
A
,
B
) is completely controllable and that the pair (
A
,
D
) is completely observable, where D
D
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 closed-loop system:
(19)
y
˙
(
t
)
=
(
A
-
B
R
-
1
B
T
x
)
y
(
t
)
=
(
A
-
B
K
)
y
(
t
)
,
where the state-feedback gain K is given by
(20)
K
=
R
-
1
B
T
x
.In order to solve (18), our purpose is to seek a matrix x on manifold PD
(
n
) such that the matrix x
B
R
-
1
B
T
x
-
x
A
-
A
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 x
B
R
-
1
B
T
x
-
x
A
-
A
T
x can be expressed as follows:
(21)
J
F
=
1
2
d
F
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
,
(22)
J
R
=
1
2
d
R
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
∥
log
⁡
(
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
)
∥
F
2
=
1
2
tr
⁡
{
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
.LetJ
:
PD
(
n
)
→
R be a scalar function of n
×
n matrix. The Riemannian gradient of J at x, denoted by ∇
x
J, is a map from PD
(
n
) to Sym
(
n
) that satisfies 〈
∇
x
J
,
s
〉
x
=
〈
∂
x
J
,
s
〉
F, for every s
∈
Sym
(
n
), where ∂
x
J is the Euclidean gradient of J at x. It follows from the inner product (5) that [27]
(23)
∇
x
J
=
x
(
∂
x
J
)
x
.Next, we will state some known facts, which are very helpful in computing∂
x
J
F and ∇
x
J
R.Lemma 1 (see [28]).
LetA
,
B denote two constant matrices and let X
,
Y denote matrix functions. One has the following properties: (1)
d
tr
(
X
)
=
tr
⁡
(
d
X
),
(2)
d
(
A
)
=
0,
(3)
d
(
A
X
B
)
=
A
(
d
X
)
B,
(4)
d
(
X
±
Y
)
=
d
X
±
d
Y,
(5)
d
(
X
Y
)
=
(
d
X
)
Y
+
X
(
d
Y
),
whered denotes the differential of a matrix (scalar) function.Lemma 2 (see [28]).
Letf
(
x
) be a scalar function of n
×
n matrix x, if
(24)
d
f
(
x
)
=
tr
⁡
(
W
d
x
)
holds; then
(25)
∂
x
f
(
x
)
=
W
T
,
where ∂
x
f
(
x
) denotes the gradient of f
(
x
) and W denotes the derivative of f
(
x
).Theorem 3.
LetJ
F
(
x
) be a scalar function of n
×
n matrix x, if
(26)
J
F
(
x
)
=
1
2
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
holds; then
(27)
∂
x
J
F
(
x
)
=
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
.Proof.
SinceQ
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x is symmetric, we have from Lemma 1 that
(28)
d
J
F
(
x
)
=
1
2
d
{
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
}
=
1
2
tr
⁡
[
d
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
=
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
d
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
]
=
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
(
d
x
A
+
A
T
d
x
-
d
x
B
R
-
1
B
T
x
-
x
B
R
-
1
B
T
d
x
)
]
=
tr
⁡
{
[
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
(
A
T
-
x
B
R
-
1
B
T
)
]
d
x
[
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
}
.
Using the above Lemma 2, we have
(29)
∂
x
J
F
(
x
)
=
[
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
]
T
=
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
.
Therefore, it is shown that (27) is valid.Theorem 4.
LetJ
R
(
x
) be a scalar function of n
×
n matrix x, if
(30)
J
R
(
x
)
=
1
2
tr
⁡
{
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
holds; then
(31)
∇
x
J
R
(
x
)
=
x
{
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
}
x
.Proof.
Sincex
B
R
-
1
B
T
-
x
A
-
A
T
x is symmetric, it follows immediately from Lemma 1 that
(32)
d
J
R
(
x
)
=
1
2
d
{
tr
⁡
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
1
2
tr
⁡
{
d
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
d
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
-
1
×
Q
-
1
/
2
d
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
(
d
x
B
R
-
1
B
T
x
-
x
B
R
-
1
B
T
d
x
-
d
x
A
-
A
T
d
x
)
}
=
tr
⁡
{
[
(
B
R
-
1
B
T
x
-
A
)
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
-
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
(
x
B
R
-
1
B
T
+
A
T
)
]
d
x
}
.
Using the above Lemma 2, we have
(33)
∂
x
J
R
(
x
)
=
{
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
(
x
B
R
-
1
B
T
+
A
T
)
]
(
B
R
-
1
B
T
x
-
A
)
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
-
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
(
x
B
R
-
1
B
T
+
A
T
)
}
T
=
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
.
Combining (23) and (33), (31) is established. This completes the proof of Theorem 4.


### 3.3. The Extended Hamiltonian Algorithm onP
D
(
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 slow-learning 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 manifoldPD
(
n
), the extended Hamiltonian algorithm can be expressed by [18]
(34)
x
˙
=
v
,
v
˙
=
v
x
-
1
v
-
∇
x
V
-
μ
v
,
where x
∈
PD
(
n
), v denotes the instantaneous learning velocity, V
:
PD
(
n
)
→
R denotes a cost function, ∇
x
V denotes the Riemannian gradient of V, and the constant μ
>
0 denotes a viscosity term. The function K
x
:
T
x
PD
(
n
)
×
T
x
PD
(
n
)
→
R denotes the kinetic energy of the particle in a point x and has the corresponding expression K
x
=
(
1
/
2
)
tr
⁡
[
(
v
x
-
1
)
2
].System (34) may be implemented by [19]
(35)
x
k
+
1
=
x
k
1
/
2
exp
⁡
(
η
x
k
-
1
/
2
v
k
x
k
-
1
/
2
)
x
k
1
/
2
,
v
k
+
1
=
η
[
v
k
x
k
-
1
v
k
-
∇
x
k
V
]
+
(
1
-
η
μ
)
v
k
,
where η is a small positive number known as the learning rate. The effectiveness of the algorithm is ensured if and only if
(36)
2
λ
max
⁡
<
μ
<
1
η
,
where λ
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 manifoldPD
(
n
), we take the Riemannian distance (22) as the cost function V
(
x
) of (18); that is,
(37)
V
(
x
)
=
1
2
d
R
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
∥
log
⁡
(
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
)
∥
F
2
.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:
(38)
x
˙
=
v
,
a
=
def
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
,
b
=
def
x
[
a
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
a
]
x
,
v
˙
=
v
x
-
1
v
-
b
-
μ
v
,
which is implemented by
(39)
x
k
+
1
=
x
k
1
/
2
exp
⁡
(
η
x
k
-
1
/
2
v
k
x
k
-
1
/
2
)
x
k
1
/
2
,
a
k
=
def
(
x
k
B
R
-
1
B
T
x
k
-
x
k
A
-
A
T
x
k
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
k
B
R
-
1
B
T
x
k
-
x
k
A
-
A
T
x
k
)
Q
-
1
/
2
]
,
b
k
=
def
x
k
[
a
k
(
x
k
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
k
+
A
)
a
k
]
x
k
,
v
k
+
1
=
η
(
v
k
x
k
-
1
v
k
-
b
k
)
+
(
1
-
η
μ
)
v
k
,
where η
>
0 denotes the learning rate and the constant μ 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 anyx belonging to the considered manifold PD
(
n
), the extended Hamiltonian algorithm is given by the following.(1)
Setx
0 as an initial input matrix x and v
0 as an initial direction of v, and then choose a desired tolerance ε
>
0.
(2)
ComputeV
(
x
k
).
(3)
IfV
(
x
k
)
<
ε then stop.
(4)
Updatex and v by (39) and return to step (2).


## 3.1. Problem Formulation

The infinite time state regulator, namely, the state equation of linear time-invariant system, is in the form(14)
y
˙
(
t
)
=
A
y
(
t
)
+
B
w
(
t
)
,
y
(
t
0
)
=
y
0
,
where y
(
t
)
∈
R
n and w
(
t
)
∈
R
m are the state and the control vector, respectively, A and B are constant matrices of appropriate dimensions, and the terminal time t
f
=
∞. The problem would be to find an optimal control w
(
t
) satisfying (14) while minimizing the quadratic performance index as follows:
(15)
J
=
1
2
∫
0
∞
‍
[
y
T
(
t
)
Q
y
(
t
)
+
w
T
(
t
)
R
w
(
t
)
]
d
t
,
where the designed parameters satisfy that both Q and R are positive definite symmetric matrices. It has been proved that the optimal control w
*
(
t
) of the cost function (15) is uniquely given by [24]
(16)
w
*
(
t
)
=
-
R
-
1
B
T
x
y
(
t
)
and the quadratic optimal performance index is gotten by
(17)
J
*
(
y
(
t
0
)
)
=
1
2
y
T
(
t
0
)
x
y
(
t
0
)
,
where x is a positive definite symmetric matrix satisfying the AREs; that is,
(18)
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
+
Q
=
0
.It is assumed that the pair(
A
,
B
) is completely controllable and that the pair (
A
,
D
) is completely observable, where D
D
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 closed-loop system:
(19)
y
˙
(
t
)
=
(
A
-
B
R
-
1
B
T
x
)
y
(
t
)
=
(
A
-
B
K
)
y
(
t
)
,
where the state-feedback gain K is given by
(20)
K
=
R
-
1
B
T
x
.In order to solve (18), our purpose is to seek a matrix x on manifold PD
(
n
) such that the matrix x
B
R
-
1
B
T
x
-
x
A
-
A
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 x
B
R
-
1
B
T
x
-
x
A
-
A
T
x can be expressed as follows:
(21)
J
F
=
1
2
d
F
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
,
(22)
J
R
=
1
2
d
R
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
∥
log
⁡
(
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
)
∥
F
2
=
1
2
tr
⁡
{
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
.LetJ
:
PD
(
n
)
→
R be a scalar function of n
×
n matrix. The Riemannian gradient of J at x, denoted by ∇
x
J, is a map from PD
(
n
) to Sym
(
n
) that satisfies 〈
∇
x
J
,
s
〉
x
=
〈
∂
x
J
,
s
〉
F, for every s
∈
Sym
(
n
), where ∂
x
J is the Euclidean gradient of J at x. It follows from the inner product (5) that [27]
(23)
∇
x
J
=
x
(
∂
x
J
)
x
.Next, we will state some known facts, which are very helpful in computing∂
x
J
F and ∇
x
J
R.Lemma 1 (see [28]).
LetA
,
B denote two constant matrices and let X
,
Y denote matrix functions. One has the following properties: (1)
d
tr
(
X
)
=
tr
⁡
(
d
X
),
(2)
d
(
A
)
=
0,
(3)
d
(
A
X
B
)
=
A
(
d
X
)
B,
(4)
d
(
X
±
Y
)
=
d
X
±
d
Y,
(5)
d
(
X
Y
)
=
(
d
X
)
Y
+
X
(
d
Y
),
whered denotes the differential of a matrix (scalar) function.Lemma 2 (see [28]).
Letf
(
x
) be a scalar function of n
×
n matrix x, if
(24)
d
f
(
x
)
=
tr
⁡
(
W
d
x
)
holds; then
(25)
∂
x
f
(
x
)
=
W
T
,
where ∂
x
f
(
x
) denotes the gradient of f
(
x
) and W denotes the derivative of f
(
x
).Theorem 3.
LetJ
F
(
x
) be a scalar function of n
×
n matrix x, if
(26)
J
F
(
x
)
=
1
2
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
holds; then
(27)
∂
x
J
F
(
x
)
=
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
.Proof.
SinceQ
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x is symmetric, we have from Lemma 1 that
(28)
d
J
F
(
x
)
=
1
2
d
{
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
}
=
1
2
tr
⁡
[
d
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
2
]
=
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
d
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
]
=
tr
⁡
[
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
(
d
x
A
+
A
T
d
x
-
d
x
B
R
-
1
B
T
x
-
x
B
R
-
1
B
T
d
x
)
]
=
tr
⁡
{
[
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
×
(
A
T
-
x
B
R
-
1
B
T
)
]
d
x
[
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
}
.
Using the above Lemma 2, we have
(29)
∂
x
J
F
(
x
)
=
[
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
]
T
=
(
A
-
B
R
-
1
B
T
x
)
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
+
(
Q
+
x
A
+
A
T
x
-
x
B
R
-
1
B
T
x
)
(
A
T
-
x
B
R
-
1
B
T
)
.
Therefore, it is shown that (27) is valid.Theorem 4.
LetJ
R
(
x
) be a scalar function of n
×
n matrix x, if
(30)
J
R
(
x
)
=
1
2
tr
⁡
{
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
holds; then
(31)
∇
x
J
R
(
x
)
=
x
{
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
}
x
.Proof.
Sincex
B
R
-
1
B
T
-
x
A
-
A
T
x is symmetric, it follows immediately from Lemma 1 that
(32)
d
J
R
(
x
)
=
1
2
d
{
tr
⁡
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
1
2
tr
⁡
{
d
log
⁡
2
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
d
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
-
1
×
Q
-
1
/
2
d
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
}
=
tr
⁡
{
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
(
d
x
B
R
-
1
B
T
x
-
x
B
R
-
1
B
T
d
x
-
d
x
A
-
A
T
d
x
)
}
=
tr
⁡
{
[
(
B
R
-
1
B
T
x
-
A
)
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
-
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
(
x
B
R
-
1
B
T
+
A
T
)
]
d
x
}
.
Using the above Lemma 2, we have
(33)
∂
x
J
R
(
x
)
=
{
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
(
x
B
R
-
1
B
T
+
A
T
)
]
(
B
R
-
1
B
T
x
-
A
)
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
-
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
(
x
B
R
-
1
B
T
+
A
T
)
}
T
=
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
×
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
.
Combining (23) and (33), (31) is established. This completes the proof of Theorem 4.


## 3.3. The Extended Hamiltonian Algorithm onP
D
(
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 slow-learning 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 manifoldPD
(
n
), the extended Hamiltonian algorithm can be expressed by [18]
(34)
x
˙
=
v
,
v
˙
=
v
x
-
1
v
-
∇
x
V
-
μ
v
,
where x
∈
PD
(
n
), v denotes the instantaneous learning velocity, V
:
PD
(
n
)
→
R denotes a cost function, ∇
x
V denotes the Riemannian gradient of V, and the constant μ
>
0 denotes a viscosity term. The function K
x
:
T
x
PD
(
n
)
×
T
x
PD
(
n
)
→
R denotes the kinetic energy of the particle in a point x and has the corresponding expression K
x
=
(
1
/
2
)
tr
⁡
[
(
v
x
-
1
)
2
].System (34) may be implemented by [19]
(35)
x
k
+
1
=
x
k
1
/
2
exp
⁡
(
η
x
k
-
1
/
2
v
k
x
k
-
1
/
2
)
x
k
1
/
2
,
v
k
+
1
=
η
[
v
k
x
k
-
1
v
k
-
∇
x
k
V
]
+
(
1
-
η
μ
)
v
k
,
where η is a small positive number known as the learning rate. The effectiveness of the algorithm is ensured if and only if
(36)
2
λ
max
⁡
<
μ
<
1
η
,
where λ
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 manifoldPD
(
n
), we take the Riemannian distance (22) as the cost function V
(
x
) of (18); that is,
(37)
V
(
x
)
=
1
2
d
R
2
(
Q
,
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
=
1
2
∥
log
⁡
(
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
)
∥
F
2
.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:
(38)
x
˙
=
v
,
a
=
def
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
B
R
-
1
B
T
x
-
x
A
-
A
T
x
)
Q
-
1
/
2
]
,
b
=
def
x
[
a
(
x
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
+
A
)
a
]
x
,
v
˙
=
v
x
-
1
v
-
b
-
μ
v
,
which is implemented by
(39)
x
k
+
1
=
x
k
1
/
2
exp
⁡
(
η
x
k
-
1
/
2
v
k
x
k
-
1
/
2
)
x
k
1
/
2
,
a
k
=
def
(
x
k
B
R
-
1
B
T
x
k
-
x
k
A
-
A
T
x
k
)
-
1
×
log
⁡
[
Q
-
1
/
2
(
x
k
B
R
-
1
B
T
x
k
-
x
k
A
-
A
T
x
k
)
Q
-
1
/
2
]
,
b
k
=
def
x
k
[
a
k
(
x
k
B
R
-
1
B
T
-
A
T
)
-
(
B
R
-
1
B
T
x
k
+
A
)
a
k
]
x
k
,
v
k
+
1
=
η
(
v
k
x
k
-
1
v
k
-
b
k
)
+
(
1
-
η
μ
)
v
k
,
where η
>
0 denotes the learning rate and the constant μ 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 anyx belonging to the considered manifold PD
(
n
), the extended Hamiltonian algorithm is given by the following.(1)
Setx
0 as an initial input matrix x and v
0 as an initial direction of v, and then choose a desired tolerance ε
>
0.
(2)
ComputeV
(
x
k
).
(3)
IfV
(
x
k
)
<
ε then stop.
(4)
Updatex 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,η
=
0.01, μ
=
3, and ε
=
1
0
-
8. Moreover, A
,
B
,
Q
,
R, and D used in this simulation are as follows:
(40)
A
=
(
0
0
1
0
)
,
B
=
(
1
0
)
,
Q
=
(
2
0
0
2
)
,
R
=
1
2
,
D
=
(
2
0
0
2
)
,
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 η, 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 η is convenient to be obtained experimentally. Based on (36), the constant μ
>
0 depends on the selection of initial value. That is to say, different initial values will be corresponding to different μ; hence the best μ is also obtained by the experiment.
As Figure1 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:
(41)
(
1.4142136
0.9999999
1.0000000
2.8284271
)
.
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.Figure 1
Comparison of convergence speeds of the EGA, the RGA, the NSIM, and the EHA.Figure 2
Kinetic energy in the EHA.Example 2.
A third-order ARE is considered following the similar procedure as Example1. In this experiment, η
=
0.01, μ
=
2.5, and ε
=
1
0
-
8, and
(42)
A
=
(
0
1
0
1
0
1
0
1
0
)
,
B
=
(
1
0
0
)
,
Q
=
(
1
0
0
0
1
0
0
0
1
)
,
R
=
1
,
D
=
(
1
0
0
0
1
0
0
0
1
)
.
As Figure3 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:
(43)
(
3.6502815
6.1622776
5.0644951
6.1622776
13.7792717
12.3245553
5.0644951
12.3245553
12.3650582
)
.Figure 3
Comparison of convergence speeds of the EGA, the RGA, the NSIM, and the EHA.Moreover, the kinetic energy function also converges to 0 in the Figure4.Figure 4
Kinetic energy in the EHA.


## 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.



---

*Source: 693659-2014-05-20.xml*