# On the Convergence of the Modified Riccati Equation

Authors: Nicholas Assimakis; Maria Adam
Journal: ISRN Signal Processing
 (2012)

---

## Abstract

The modified Riccati equation arises in the implementation of Kalman filter in target tracking under measurement uncertainty and it cannot be transformed into  an equation of the form of the Riccati equation. An iterative solution algorithm of the modified Riccati equation is proposed. A method is established to decide when the proposed algorithm is faster than the classical one. Both algorithms have the same behavior: if the system is stable, then there exists a steady-state solution, while if the system is unstable, then there exists a critical value of the measurement detection probability, below which both iterative algorithms diverge. It is established that this critical value increases in a logarithmic way as the system becomes more unstable.


---

## Body

## 1. Introduction

The discrete time modified Riccati equation emanating from Kalman filter was originally formulated in [1]. It plays an important role in target tracking [1–10]. Theoretical properties of the modified Riccati equation have been derived in [2, 3]. It is well known [2] that the modified Riccati equation cannot be transformed into an equation of the form of the Riccati equation. The discrete time Riccati equation arises in linear estimation, namely, in the implementation of the discrete time Kalman filter [11]. The modified Riccati equation is solvable under certain conditions [2, 9] and has existence and uniqueness properties similar to the Riccati equation [2].In Section2, the modified Riccati equation associated with target tracking under measurement uncertainty is presented; the case without clutter but with detection probability of less than one is considered. In Section 3, an iterative solution algorithm of the modified Riccati equation is proposed and compared to the classical one. A method is established to distinguish the faster algorithm. If the system is stable, then both algorithms do converge to the steady-state solution. If the system is unstable, then there exists a critical value of the measurement detection probability, below which both algorithms diverge. In Section 4, it is established that this critical value increases in a logarithmic way as the system becomes more unstable.


## 2. The Modified Riccati Equation

Consider the following state space equations at timek≥0:
(1)xk+1=Fxk+wk,zk=Hxk+υk,
where xk is the n×1 state vector at time k, zk is the m×1 measurement vector, F is the n×n system transition matrix, and H is the m×n output matrix. It is assumed that {wk} and {υk} are zero mean, independent, white, Gaussian noise processes with constant covariance matrices given by E[wkwkT]=Q and E[υkυkT]=R, that is, Q  is the n×n plant noise covariance matrix R is the m×m measurement noise covariance matrix.At the initial timek=0, the state x0 is independent of the processes {wk} and {υk} for any k and x0 is a Gaussian random variable with mean x-0 and covariance P0, that is,
(2)E[x0]=x-0,E[(x0-x-0)(x0-x-0)T]=P0.
For k≥0, denoting Zk={z0,z1,…,zk}, the state prediction xk/k-1=E[xk/Zk-1] and the prediction error covariance matrix Pk/k-1=E[(xk-xk/k-1)(xk-xk/k-1)T/Zk-1], and using the discrete-time invariant Kalman filter equations as described in [11–13], we derive the following recursion for the symmetric n×n prediction error covariance matrix Pk+1/k, the Riccati equation:
(3)Pk+1/k=Q+FPk/k-1FT-FPk/k-1HT[HPk/k-1HT+R]-1HPk/k-1FT
with initial condition P0/-1=P0.The state space equations in (1) can be used in target tracking to describe a linear target motion and measurement model. In addition, consider that a measurement is received with detection probability pd [2], where
(4)0≤pd≤1.Using the Kalman filter equations we are able to derive [2, 3] a Kalman-like recursion for the symmetric n×n prediction error covariance matrix Pk+1/k, the modified Riccati equation:
(5)Pk+1/k=Q+FPk/k-1FT-pdFPk/k-1HT[HPk/k-1HT+R]-1HPk/k-1FT.
Note that the prediction error covariance matrix is nonnegative definite (Pk+1/k≥0).Also in the modified Riccati equation (5) note that the m×m matrix HPk/k-1HT+R is nonsingular, when R is a positive definite matrix (R>0), which has the significance that no measurement is exact; this is reasonable in physical problems.It is remarkable that the following special cases are implied by (5).(i)
Settingpd=1, the classical Riccati equation (3) is derived. The difference between the modified Riccati equation and the Riccati equation is the term of detection probability pd. It is obvious that the modified Riccati equation (5) cannot be transformed into the classical Riccati equation (3).(ii)
Settingpd=0 in (5), the classical Lyapunov equation is derived:
(6)Pk+1/k=Q+FPk/k-1FT
which arises from the Riccati equation (3) in the infinite measurement noise case (R→∞).


## 3. Iterative Solutions of the Modified Riccati Equation

Concerning the modified Riccati equation, it is known [9] that for stable systems, which means that all eigenvalues of F lie inside the unit circle, the modified Riccati equation always converges  and the limiting value P of the prediction error covariance is the steady state solution of the discrete time modified Riccati equation.Theclassical implementation of the modified Riccati Equation (cmRE) arises from (5), which consists of the direct implementation of the recursion of the following equation:
(7)Pk+1/k=F{11Pk/k-1-pdFPk/k-1HT×[HPk/k-1HT+R]-1HPk/k-1}FT+Q.
It is obvious that this equation is equivalent to the modified Riccati equation (5) achieving a reduction in computational burden by using as a common factor F.Notice that ifR>0, then the nonsingularity of HPk/k-1HT+R in (7) is guaranteed.It is known [2] that the steady state solution of the modified Riccati equation is independent of the initial condition  P0/-1=P0. So, for convenience, we are able to use zero initial condition P0=0. Then we are able to use P1/0=Q as initial condition for the classical implementation.Note that the convergence is achieved, when∥Pk+1/k-Pk/k-1∥<ε, where ε is a small positive number and ∥M∥ denotes the norm of the matrix M, which is equal to the largest singular value of M. Then the steady state solution P satisfies the steady state modified Riccati equation:
(8)P=Q+FPFT-pdFPHT[HPHT+R]-1HPFT.Theproposed implementation of modified Riccati Equation (pmRE) consists of the direct implementation of the recursion of the following equation:
(9)Pk+1/k=(1-pd)(Q+FPk/k-1FT)+pd(Q+F[Pk/k-1-1+HTR-1H]-1)FT.
This equation is equivalent to the modified Riccati equation (5) and can be derived from (5) using the matrix inversion lemma, under the condition that for k≥1 the prediction error covariance matrix Pk/k-1 is positive definite (Pk/k-1>0). This is guaranteed, if  Q>0, due to the fact that P1/0=Q, if we use zero initial condition P0/-1=P0=0.Notice that ifR>0, then the existence of HTR-1H in (9) is guaranteed and the nonsingularity of Pk/k-1-1+HTR-1H becomes obvious. Also note that HTR-1H in (9) is computed once (initialization process).Both the classical and the proposed algorithms for solving the modified Riccati equation are recursive ones. Thus, the total computational time required for the implementation of each algorithm is(10)ta=BaSato,
where Ba is the per recursion calculation burden required for the online calculations of each algorithm, Sa is the number of recursions (steps) that each algorithm executes, and to is the time required to perform a scalar operation.Note that the two algorithms are equivalent to each other with respect to their behaviour: they calculate theoretically such steady-state prediction error varianceP that satisfies (8). Then, it is reasonable to assume that both algorithms compute the limiting solution of the modified Riccati equation executing the same number of recursions, depending on the desired accuracy. Thus, in order to compare the algorithms with respect to their computational time, we have to compare their per recursion calculation burden required for the online calculations; the calculation burden of the offline calculations (initialization process) is not taken into account.The computational analysis is based on the calculation burden of the matrix operations, which are summarized in Table1 and needed for the implementation of the filtering algorithm [14].Table 1
Calculation burden of matrix operations.
Matrix operation
Calculation burden
A
(
n
×
n
)
+
B
(
n
×
n
)
=
S symmetric
(
1
/
2
)
(
n
2
+
n
)
c
·
S
(
n
×
n
) symmetric
(
1
/
2
)
(
n
2
+
n
)
A
(
n
×
m
)
·
B
(
m
×
k
)
=
C
(
n
×
k
)
2
n
m
k
-
n
k
A
(
n
×
m
)
·
B
(
m
×
n
)
=
S symmetric
n
2
m
+
n
m
-
(
1
/
2
)
(
n
2
+
n
)
[
A
(
n
×
n
)
]
-
1
=
B
(
n
×
n
)
(
1
/
6
)
(
16
n
3
-
3
n
2
-
n
)Then, the (per recursion) computational requirements of the classical and the proposed algorithms for solving the modified Riccati equation are computed as(11)BcmRE=12(6n3+n2+n)+3n2m+3nm2+16(16m3-3m2-m),BpmRE=13(34n3+3n2+5n).
The details of (11) are given in Tables 2 and 3.Table 2
Classical implementation of the modified Riccati equation (cmRE).
Matrix operation
Matrix dimensions
Calculation burden
H
·
P
k
/
k
-
1
(
m
×
n
)
·
(
n
×
n
)
2
n
2
m
-
n
m
H
P
k
/
k
-
1
·
H
T
(
m
×
n
)
·
(
n
×
m
)
*
n
m
2
+
n
m
-
(
1
/
2
)
(
m
2
+
m
)
O
k
≡
H
P
k
/
k
-
1
H
T
+
R
(
m
×
m
)
+
(
m
×
m
)
*
(
1
/
2
)
(
m
2
+
m
)
O
k
-
1
(
m
×
m
)
(
1
/
6
)
(
16
m
3
-
3
m
2
-
m
)
K
k
≡
O
k
-
1
·
H
P
k
/
k
-
1
(
m
×
m
)
·
(
m
×
n
)
2
n
m
2
-
n
m
W
≡
P
k
/
k
-
1
H
T
·
K
k
(
n
×
m
)
·
(
m
×
n
)
*
n
2
m
+
n
m
-
(
1
/
2
)
(
n
2
+
n
)
p
d
·
W
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
Π
1
≡
P
k
/
k
-
1
-
p
d
·
W
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
F
·
Π
1
(
n
×
n
)
·
(
n
×
n
)
2
n
3
-
n
2
Π
2
≡
F
Π
1
·
F
T
(
n
×
n
)
·
(
n
×
n
)
*
n
3
+
(
1
/
2
)
(
n
2
-
n
)
P
k
+
1
/
k
=
Q
+
Π
2
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
B
cmRE
=
(
1
/
2
)
(
6
n
3
+
n
2
+
n
)
+
3
n
2
m
+
3
n
m
2
+
(
1
/
6
)
(
16
m
3
-
3
m
2
-
m
)
*Symmetric matrix.Table 3
Proposed implementation of the modified Riccati equation (pmRE).
Matrix operation
Matrix dimensions
Calculation burden
P
k
/
k
-
1
-
1
(
n
×
n
)
(
1
/
6
)
(
16
n
3
-
3
n
2
-
n
)
M
≡
P
k
/
k
-
1
-
1
+
H
T
R
-
1
H
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
M
-
1
(
n
×
n
)
(
1
/
6
)
(
16
n
3
-
3
n
2
-
n
)
F
·
M
-
1
(
n
×
n
)
·
(
n
×
n
)
2
n
3
-
n
2
F
M
-
1
·
F
T
(
n
×
n
)
·
(
n
×
n
)
*
n
3
+
(
1
/
2
)
(
n
2
-
n
)
Ω
≡
Q
+
F
M
-
1
F
T
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
p
d
·
Ω
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
F
·
P
k
/
k
-
1
(
n
×
n
)
·
(
n
×
n
)
2
n
3
-
n
2
F
P
k
/
k
-
1
·
F
T
(
n
×
n
)
·
(
n
×
n
)
*
n
3
+
(
1
/
2
)
(
n
2
-
n
)
Ψ
≡
Q
+
F
P
k
/
k
-
1
F
T
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
(
1
-
p
d
)
·
Ψ
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
P
k
+
1
/
k
=
(
1
-
p
d
)
Ψ
+
p
d
Ω
(
n
×
n
)
+
(
n
×
n
)
*
(
1
/
2
)
(
n
2
+
n
)
B
pmRE
=
(
1
/
3
)
(
34
n
3
+
3
n
2
+
5
n
)
*Symmetric matrix.From the above computational requirements where we derive the following conclusions.(1)
The per recursion calculation burden of the classical algorithm depends on the state dimensionn and on the measurement dimension m, while the per recursion calculation burden of the proposed algorithm depends only on the state dimension n.(2)
The proposed algorithm is faster than the classical one if the following relation holds:(12)BcmRE-BpmRE=3n2m+3nm2+16(16m3-3m2-m)-16(50n3+3n2+7n)>0.Figure 1 depicts the relation between the dimensions n and m that may hold in order to decide which algorithm is faster. Then, relation (12) is approached by the relation:
(13)m>n.
Thus it becomes obvious that we are able to establish the following method to distinguish the faster algorithm:“ifm>n, then the proposed algorithm is faster than  the classical algorithm, else the classical algorithm is faster than the proposed algorithm”.Figure 1
Proposed algorithm may be faster than the classical one.


## 4. Convergence of the Modified Riccati Equation

Concerning the modified Riccati equation, it is known [9] that, for unstable systems (there is at least one eigenvalue of F that lies strictly outside the unit circle), there exists a critical value of detection probability pd, below which the modified Riccati equation diverges.Simulation results were taken concerning the modified Riccati equation. Both the classical and the proposed algorithms were implemented in order to solve the modified Riccati equation. Various stable as well as unstable models were considered. The following results were confirmed.(i)
Both algorithms have the same behavior. If the system is stable (all eigenvalues ofF lie inside the unit circle), then the modified Riccati equation always converges: there always exists a steady-state solution. If the system is unstable (there is at least one eigenvalue of F that lies strictly outside the unit circle), then there exists a critical value of detection probability pd, below which the modified Riccati equation diverges.(ii)
The critical value of detection probabilitypd increases in a logarithmic way as the maximum absolute eigenvalue of F increases, that is, the system becomes more unstable. Figure 2 depicts the relation between the system stability and the modified Riccati equation convergence.(iii)
In the special case where the maximum absolute eigenvalue ofF lies in the unit circle, the critical value of detection probability takes its minimum value of the order of 0.04.(iv)
The maximum absolute eigenvalue ofF, below which the modified Riccati equation always diverges, is of the order of 10.Figure 2
Modified Riccati equation convergence with respect to system stability.



---

*Source: 625897-2012-07-08.xml*