On the Convergence of the Modified Riccati Equation

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.

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 Section 2, 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 time 𝑘≥0: 𝑥𝑧_𝑘+1𝑘=𝐹𝑥_𝑘+𝑤_𝑘,=𝐻𝑥_𝑘+𝜐_𝑘,(1) where 𝑥_𝑘 is the 𝑛×1 state vector at time 𝑘, 𝑧_𝑘 is the 𝑚×1 measurement vector, 𝐹 is the 𝑛×𝑛 system transition matrix, and 𝐻 is the 𝑚×𝑛 output matrix. It is assumed that {𝑤_𝑘} and {𝜐_𝑘} are zero mean, independent, white, Gaussian noise processes with constant covariance matrices given by 𝐸[𝑤_𝑘𝑤^𝑇𝑘]=𝑄 and 𝐸[𝜐_𝑘𝜐^𝑇𝑘]=𝑅, that is, 𝑄  is the 𝑛×𝑛 plant noise covariance matrix 𝑅 is the 𝑚×𝑚 measurement noise covariance matrix.

At the initial time 𝑘=0, the state 𝑥₀ is independent of the processes {𝑤_𝑘} and {𝜐_𝑘} for any 𝑘 and 𝑥₀ is a Gaussian random variable with mean 𝑥⁻₀ and covariance 𝑃₀, that is, 𝐸𝑥₀=𝑥⁻₀,𝐸𝑥₀−𝑥⁻₀𝑥⁰−𝑥⁻₀^𝑇=𝑃⁰.(2) For 𝑘≥0, denoting 𝑍_𝑘={𝑧₀,𝑧₁,…,𝑧_𝑘}, the state prediction 𝑥_{𝑘/𝑘−1}=𝐸[𝑥_𝑘/𝑍_𝑘−1] and the prediction error covariance matrix 𝑃_{𝑘/𝑘−1}=𝐸[(𝑥_𝑘−𝑥_{𝑘/𝑘−1})(𝑥_𝑘−𝑥_{𝑘/𝑘−1})^𝑇/𝑍_𝑘−1], and using the discrete-time invariant Kalman filter equations as described in [11–13], we derive the following recursion for the symmetric 𝑛×𝑛 prediction error covariance matrix 𝑃_𝑘+1/𝑘, the Riccati equation: 𝑃_𝑘+1/𝑘=𝑄+𝐹𝑃_{𝑘/𝑘−1}𝐹^𝑇−𝐹𝑃_{𝑘/𝑘−1}𝐻^𝑇𝐻𝑃_{𝑘/𝑘−1}𝐻^𝑇+𝑅⁻¹𝐻𝑃_{𝑘/𝑘−1}𝐹^𝑇(3) with initial condition 𝑃_0/−1=𝑃₀.

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 𝑝_𝑑 [2], where 0≤𝑝_𝑑≤1.(4)

Using the Kalman filter equations we are able to derive [2, 3] a Kalman-like recursion for the symmetric 𝑛×𝑛 prediction error covariance matrix 𝑃_𝑘+1/𝑘, the modified Riccati equation: 𝑃_𝑘+1/𝑘=𝑄+𝐹𝑃_{𝑘/𝑘−1}𝐹^𝑇−𝑝_𝑑𝐹𝑃_{𝑘/𝑘−1}𝐻^𝑇𝐻𝑃_{𝑘/𝑘−1}𝐻^𝑇+𝑅⁻¹𝐻𝑃_{𝑘/𝑘−1}𝐹^𝑇.(5) Note that the prediction error covariance matrix is nonnegative definite (𝑃_𝑘+1/𝑘≥0).

Also in the modified Riccati equation (5) note that the 𝑚×𝑚 matrix 𝐻𝑃_{𝑘/𝑘−1}𝐻^𝑇+𝑅 is nonsingular, when 𝑅 is a positive definite matrix (𝑅>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)Setting 𝑝_𝑑=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 𝑝_𝑑. It is obvious that the modified Riccati equation (5) cannot be transformed into the classical Riccati equation (3).(ii)Setting 𝑝_𝑑=0 in (5), the classical Lyapunov equation is derived: 𝑃_𝑘+1/𝑘=𝑄+𝐹𝑃_{𝑘/𝑘−1}𝐹^𝑇(6) which arises from the Riccati equation (3) in the infinite measurement noise case (𝑅→∞).

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 𝐹 lie inside the unit circle, the modified Riccati equation always converges  and the limiting value 𝑃 of the prediction error covariance is the steady state solution of the discrete time modified Riccati equation.

The classical implementation of the modified Riccati Equation (cmRE) arises from (5), which consists of the direct implementation of the recursion of the following equation: 𝑃_𝑘+1/𝑘=𝐹𝑃_{𝑘/𝑘−1}−𝑝_𝑑𝐹𝑃_{𝑘/𝑘−1}𝐻^𝑇×𝐻𝑃_{𝑘/𝑘−1}𝐻^𝑇+𝑅⁻¹𝐻𝑃_{𝑘/𝑘−1}𝐹^𝑇+𝑄.(7) 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 𝐹.

Notice that if 𝑅>0, then the nonsingularity of 𝐻𝑃_{𝑘/𝑘−1}𝐻^𝑇+𝑅 in (7) is guaranteed.

It is known [2] that the steady state solution of the modified Riccati equation is independent of the initial condition 𝑃_0/−1=𝑃₀. So, for convenience, we are able to use zero initial condition 𝑃₀=0. Then we are able to use 𝑃_1/0=𝑄 as initial condition for the classical implementation.

Note that the convergence is achieved, when ‖𝑃_𝑘+1/𝑘−𝑃_{𝑘/𝑘−1}‖<𝜀, where 𝜀 is a small positive number and ‖𝑀‖ denotes the norm of the matrix 𝑀, which is equal to the largest singular value of 𝑀. Then the steady state solution 𝑃 satisfies the steady state modified Riccati equation: 𝑃=𝑄+𝐹𝑃𝐹^𝑇−𝑝_𝑑𝐹𝑃𝐻^𝑇𝐻𝑃𝐻^𝑇+𝑅⁻¹𝐻𝑃𝐹^𝑇.(8)

The proposed implementation of modified Riccati Equation (pmRE) consists of the direct implementation of the recursion of the following equation: 𝑃_𝑘+1/𝑘=+1𝑝−_𝑑𝑝_𝑑𝑄+𝐹𝑃_{𝑘/𝑘−1}𝐹^𝑇𝑄+𝐹𝑃^{−1𝑘/𝑘−1}+𝐻^𝑇𝑅⁻¹𝐻⁻¹𝐹^𝑇.(9) 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 𝑘≥1 the prediction error covariance matrix 𝑃_{𝑘/𝑘−1} is positive definite (𝑃_{𝑘/𝑘−1}>0). This is guaranteed, if 𝑄>0, due to the fact that 𝑃_1/0=𝑄, if we use zero initial condition 𝑃_0/−1=𝑃₀=0.

Notice that if 𝑅>0, then the existence of 𝐻^𝑇𝑅⁻¹𝐻 in (9) is guaranteed and the nonsingularity of 𝑃^{−1𝑘/𝑘−1}+𝐻^𝑇𝑅⁻¹𝐻 becomes obvious. Also note that 𝐻^𝑇𝑅⁻¹𝐻 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) where 𝐵_𝑎 is the per recursion calculation burden required for the online calculations of each algorithm, 𝑆_𝑎 is the number of recursions (steps) that each algorithm executes, and 𝑡_𝑜 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 variance 𝑃 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 Table 1 and needed for the implementation of the filtering algorithm [14].

Table 1 

Calculation burden of matrix operations.

Then, the (per recursion) computational requirements of the classical and the proposed algorithms for solving the modified Riccati equation are computed as 𝐵_cmRE=126𝑛³+3𝑛₂+𝑛²+𝑛𝑚+3𝑛𝑚²+1616𝑚³𝐵−3𝑚²−𝑚_pmRE,=1334𝑛³+3𝑛²+5𝑛.(11) The details of (11) are given in Tables 2 and 3.

Table 2 

Classical implementation of the modified Riccati equation (cmRE).

Table 3 

Proposed implementation of the modified Riccati equation (pmRE).

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 dimension 𝑛 and on the measurement dimension 𝑚, while the per recursion calculation burden of the proposed algorithm depends only on the state dimension 𝑛.(2)The proposed algorithm is faster than the classical one if the following relation holds: 𝐵_cmRE−𝐵_pmRE=3𝑛²𝑚+3𝑛𝑚²+1616𝑚³−1−3𝑚²−𝑚650𝑛³+3𝑛²+7𝑛>0.(12)Figure 1 depicts the relation between the dimensions 𝑛 and 𝑚 that may hold in order to decide which algorithm is faster. Then, relation (12) is approached by the relation: 𝑚>𝑛.(13) Thus it becomes obvious that we are able to establish the following method to distinguish the faster algorithm: “if 𝑚>𝑛, 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 𝐹 that lies strictly outside the unit circle), there exists a critical value of detection probability 𝑝_𝑑, 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 of 𝐹 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 𝐹 that lies strictly outside the unit circle), then there exists a critical value of detection probability 𝑝_𝑑, below which the modified Riccati equation diverges.(ii) The critical value of detection probability  𝑝_𝑑 increases in a logarithmic way as the maximum absolute eigenvalue of 𝐹 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 of 𝐹 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 of 𝐹, 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.