Maximum Correntropy High-Order Extended Kalman Filter

I.   Introduction

Filter has been applied in various fields, including navigation, national defense construction, neural network training and so on^[1-6]. For general nominal linear systems, Kalman filter (KF) provides an optimal solution step-by-step^[7]. For nonlinear problems, extended Kalman filter is the most typical method, which approximates the nonlinear system by its first order linearization^[8]. Another nonlinear filters are based on sigma sampling, such as unscented Kalman filter (UKF), cubature Kalman filter (CKF), and so on^[9-12], which have better filtering performance compared with EKF.

KF and its extensions in general have been established to meet the conditions of the standard KF^[13-15]. Nevertheless, their performances may get worse when employed to non-Gaussian situations. This is because their objective function and criterion depend on large outliers, they are not suitable for non-Gaussian environments^[16].

For non-Gaussian systems, Chen proposed maximum correntropy KF (MCKF)^[17]. MCKF only needs a limited number of realizations for modeling error to establish a mean estimator of correntropy. Lately, the Maximum correentropy extended KF (MCEKF) that can solve nonlinear non-Gaussian noise systems have appeared based on EKF^[16]. However, similar to EKF, MCEKF can realize the approximation of the nonlinear function not higher than the second order.

In the paper, we develop a novel high-order extended KF, called the maximum correntropy high-order KF (H-MCEKF) by combining MCC with a fixed-point iterative algorithm. Similar to EKF, the H-MCEKF not only retains the propagation process of state, but that of covariance matrix. Therefore, the new algorithm has two significant advantages: real-time and recursion. Unlike the first-order Taylor expansion of MCEKF, HMCEKF uses polynomials to have a full description for nonlinear functions. It should be emphasized that it is not difficult to convert a nonlinear function into a high-order polynomial by employing Taylor expansion or multi-dimensional Taylor network^[17].

There are two contributions that can be shown in this paper. 1) The idea of converting a nonlinear function to a linear form is given: a) all high-order polynomials in the system are defined as implicit variables and treated as parameter variables; b) the original state model is equivalently formulated into a pseudo-linear form; c) an augmented linear state model is established by combing with pseudo-linear state model; d) the original measurement model is rewritten into linear form. 2) The solution of high-order polynomials are converted into KF.

This paper will be divided into several sections to describe and derive the design process of the new filter: Brief introduction in Section II; modeling the nonlinear non-Gaussian systems in Section III; linearizing the systems in Section IV; derivation process of H-MCEKF in Section V; simulations in Section VI and conclusions and future in Section VII finally.

II.   A Description of Correntropy

For one-dimensional random variable $X,Y$ with ${F_{XY}}(x,y)$, then the definition of correntropy is^[17]

where, E is expectation, and $\Upsilon (XY)$ is Mercer Kernel. For example, Gaussian kernel is as follows.

where, e is the difference between x and y, $\sigma$ is their covariance.

Perform Taylor expansion on Eq.(2),

then, the correntropy of Eq.(1) can be formulated

where, $E\left\{ {{{(X - Y)}^{2p}}} \right\} = \int {{{(x - y)}^{2p}}d{F_{XY}}(x,y)}$ is $2p$ order moment of the random variable $X,Y \in R$

In general, random variable pairs may be relatively easy to obtain, and ${F_{XY}}$ is difficult to obtain. In this case, we can get the correntropy by Eq.(5)

then, the correntrpy of random variable pairs $(X,Y)\,$ under limited data is

When $X,Y \in {R^n}$ and the components in $e =$$X - Y$ are independent of each other, the multi-dimensional correntropy can be obtained by Eq.(7)

where, $i$ denotes ${i{{\rm{th}}}}$ component of $e$; $j$ denotes ${j{{\rm{th}}}}$ sampler.

Remark 1　The larger the $\sigma$ is , the greater the influence of the second moment on the correntropy is.

III.   Description of the System

Given a class of state models and observation models with strong nonlinear characteristics

where, $x(u)$ and $y(u)$ are n-dimensional state variable and m-dimensional measurement variable, respectively; state modeling errors $\omega (u)$ and measurement modeling error $\nu (u)$ are both with non-Gaussian systems; $f(x(u))$ is nonlinear state function and $h(x(u + 1))$ represents nonlinear function of $x(u + 1)$.

IV.   Linearization of Nonlinear Models

1   Pseudo-linearized representation of nonlinear functions

For ease of understanding, we simplify the representation of the high-dimensional linearization process based on Eq.(8) and Eq.(9).

For facilitate understanding, we only describe and establish the filter for the two-dimensional systems and suppose.

where, ${a_{i,{l_1},{l_2}}}$ is the corresponding weight value for $\displaystyle\sum\limits_{ {l_1} + {l_2} = l \atop {l_1},{l_2} \leqslant l} {x_1^{{l_1}}(u)x_2^{{l_2}}(u)}$.

Definition 1　${x^{(l)}}(u): = \{ x_1^{{l_1}}(u)x_2^{{l_2}}(u),\;\;\;{l_1} + {l_2} = l;$$0 \leqslant {l_j}\leqslant l; l = 0,1, \ldots ,r\}$ is a set of hidden variables with sequential order.

Definition 2　$a_i^{(l)}$ is weight vector for ${x^{(l)}}(u)$ where

On the basis of Definition 1, 2 and Eq.(10), we have

Let

then,

Similarly, suppose the measurement function in Eq.(9) is shown as follows

Similar to Definition 1 and Definition 2, Eq.(10) and Eq.(11), Eq.(13) has the matrix form as following.

2   Linearized representation of nonlinear functions

In order to linearize the nonlinear functions, we build the dynamic model as follows

where, $A_l^{(q)}(u)$ can be solved by given prior information. Specially, we may suppose that

Combining Definition 1 and Definition 2, Eq.(11) and Eq.(12), state model Eq.(8) has further linear form

Let

Eq.(17) is equivalently rewritten as follows

where, $\underline W (u)$ denotes the modeling uncertainty of the augmented linearized state.

Similarly, the linear matrix form of the measurement model (9) is

On the basis of Eq.(19), we can obtain the linearization of Eq.(2).

where, $V(k)$ is the modeling error.

V.   Design of Maximum Correntropy High-Order Extended Kalman Filter

1   Characterization of non-Gaussian modeling error

For linear models Eq.(18) and Eq.(20), we have

where, $\underline {\hat X} (u + 1|u) = \underline A (u + 1,k)\underline {\hat X} (u|u)$ can be obtained by Eq.(18).

where, $- \underline {\tilde X} (u + 1|u){\rm{ = }} \underline X (u + 1|u + 1) - \underline {\hat X}(u + 1|u)$

and

where, $E [ \cdot ]$ can be obtained by finite sampling.

where

${\omega ^{(1,j)}}(u)$ is the $j{\rm{th}}$ sampling of modeling error ${\omega ^{(1)}}(u)$; and ${\omega ^{(l)}}(u) \sim N(0,{Q^{(l)}}),l = 2,3, \ldots ,r$.

Similarly, ${R_V}(u + 1)$ is the covariance matrix of $v(u+1)$ in original state model (9).

where $\bar v(u + 1) = \dfrac{1}{N}\displaystyle\sum\limits_{j = 1}^N {{v^{(j)}}(u + 1)}$, ${v^{(j)}}(u + 1)$ is the $j{{\rm{th}}}$ sampling of the non-Gaussian random noise vector $v(u + 1)$.

2   Statistical independence process of non-Gaussian modeling error vector ${\boldsymbol{u(k+1)}}$

$u$ in Eq.(23) is $L =\displaystyle \sum\limits_{l = 1}^r {{n_l}} + m$ dimension modeling error vector, and the relationship between components are not independent. In order to employ the correntropy shown in Eq.(7), $u(u+1)$ needs to be transformed independently. Then, Eq.(23) can be further formulated.

where, ${\underline B _X}(u + 1)$ and ${\underline B _Y}(u + 1)$ can be obtained by Cholesky decomposition of $\underline P (k + 1)$ and ${\underline R _V}(u + 1)$, respectively.

Combing Eq.(26) and ${\underline B ^{ - 1}}(u + 1)$, Eq.(21) can be simplified as

where

with

Therefore, the components of the obtained random variable $e(u + 1)$ become statistically independent after the non-Gaussian modeling error $u(u + 1)$ is equivalently transformed by ${B^{ - 1}}(u + 1)$.

3   Cost function for the optimal state estimation

We propose the following objective function for solving $\underline {\hat X} (u + 1|u + 1)$ by Eq.(21)

where, $d_i^{(j)}$ is the ${j{{\rm{th}}}}$ realization of the ${i{{\rm{th}}}}$ element of $D(k)$,

4   Fix-point iterative algorithm for state estimation under MCC

According to Eq.(30), we can get the optimal solution.

Further

Considering that $e_i^{(j)}(u + 1) = d_i^{(j)}(u + 1) - {s_i}(u + 1)\underline X (u +$$1)$, Eq.(31) can be written as a fixed point equation about $\underline X (u + 1)$

Let

then

where

5   Equivalent conversion from fixed point solution equation to Kalman filter

From Ref [14], Eq.(34) and Eq.(35), we arrive at

where

Further, we get

where

Thus, the equivalent conversion is completed from fixed-point equation to Kalman filter.

6   Online iterative solution of H-MCEKF

With the above derivations, we summarize the proposed H-MCEKF algorithm. Select a proper $\sigma$ and $\varepsilon$, where $\sigma$ controls the local scope of the Gaussian kernel function and $\varepsilon$ represents the estimation accuracy. Set an initial estimate $\underline {\hat X} {(k |k )_0} = \underline {\hat X} (k |k-1)$, $P(k|k)$, then

where

If Eq.(47) holds,

Remark 2 Without confusion, let $u + 1: = u$

Remark 3 Eq.(44) and Eq.(45) contribute to the adjustment of the noise uncertainty. As the system iterates, Eq.(44) and Eq.(45) will gradually converge to 1.

Remark 4 In the polynomial expansion method in the paper, as the dimensionality increases, higher-order terms will show sparseness and its proportion will be less and less. We can use the method of pruning to model it again, so that the dimensionality can be controlled, and at the same time, the dynamic statistical properties of high-order terms will be retained.

VI.   Simulations

This chapter verifies the performance of the proposed new filter through several simulations.

1   Example 1

Giving the following systems shown in Eq.(8) and Eq.(9)

where

The simulation is performed many times from 1-500, where $w(u) \sim 0.8N(0,0.01) + 0.2N(0,0.1)$ for the first 200 times and $w(u + 1)\sim 0.8N(0,0.01) + 0.2N(0,0.1) +$${[\begin{array}{*{20}{c}} {0.1}&{0.1} \end{array}]^{\rm{T}}}$ for the remaining; ${v}(u+1) \sim N(0,0.01)$; ${x^{(2)}}(u) = {\left[ {\begin{array}{*{20}{c}} {x_1^2(u)}&{x_2^2(u)} \end{array}} \right]^{\rm{T}}}$, ${x^{(3)}}(u) = {\left[ {\begin{array}{*{20}{c}} {x_1^3(u)}&{x_2^3(u)} \end{array}} \right]^{\rm{T}}}$, ${x^{(4)}}(u) = {\left[ {\begin{array}{*{20}{c}} {x_1^4(u)}&{x_2^4(u)} \end{array}} \right]^{\rm{T}}}$, and ${x^{(5)}}(u) = \left[ {c} {x_1^5(u)}\;\; {x_2^5(u)}\right]^{\rm{T}}$ are defined as implicit variables and their modeling errors satisfies that the mean is 0 and the variance is 0.01; $x(0) = {[1\;\;1]^{\rm{T}}}$, $\hat x(0|0) = {[1\;\;1]^{\rm{T}}} + N(0,0.01) \times randn(2,1)$ and $P(0|0) = 0.01 \times diag(1,1)$

Case 1 only considers the process noise is non-Gaussian and takes $\sigma = 10,\varepsilon = {10^{ - 1}}$ as an example to demonstrate the estimation results of MCEKF and H-MCEKF, as shown in Fig.1(a)-(d). Table 1 summarizes the estimation errors under different $\sigma$ and $\varepsilon$.

Figure  1.  The actual state $x_1$ and its estimate (a) with estimate error (c); The actual state $x_2$ and its estimate (b) with estimate error (d)

DownLoad: Full-Size Img PowerPoint

Table  1.  Estimation errors with different filters

		MSE of ${x_{\rm{1}}}$			MSE of ${x_{\rm{2}}}$			MSE
$\sigma$	$\varepsilon$	MCEKF	H-MCEKF	Improved	MCEKF	H-MCEKF	Improved	MCEKF	H-MCEKF	Improved
$\sigma {\rm{ = }}2$	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.0469	0.0467	0.42%	0.0226	0.0204	9.73%	0.0348	0.0335	3.74%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.0434	0.0421	3.00%	0.0344	0.0256	25.58%	0.0389	0.0338	13.11%
$\sigma {\rm{ = }}5$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.0205	0.0189	7.80%	0.0091	0.0089	2.20%	0.0148	0.0139	6.08%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.0107	0.0099	7.47%	0.0112	0.0083	25.89%	0.0110	0.0091	17.27%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.0201	0.0178	11.44%	0.0235	0.0195	17.02%	0.0218	0.0187	14.22%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.0203	0.0111	45.32%	0.0165	0.0148	10.30%	0.0184	0.0130	29.34%
$\sigma {\rm{ = }}10$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.0128	0.0105	17.96%	0.0107	0.0083	22.42%	0.0117	0.0094	19.65%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.0150	0.0140	6.67%	0.0147	0.0131	18.88%	0.0148	0.0136	8.11%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.0133	0.0092	30.82%	0.0122	0.0108	11.48%	0.0127	0.0100	21.26%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.0190	0.0169	11.05%	0.0155	0.0116	25.16%	0.0173	0.0142	17.92%
$\sigma {\rm{ = }}15$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.0081	0.0071	12.35%	0.0083	0.0061	26.50%	0.0082	0.0066	19.51%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.0092	0.0074	19.56%	0.0077	0.0055	28.57%	0.0084	0.0064	23.80%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.0080	0.0069	13.75%	0.0080	0.0057	28.75%	0.0080	0.0063	21.25%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.0086	0.0073	15.11%	0.0084	0.0053	36.90%	0.0085	0.0063	25.88%

 | Show Table

DownLoad: CSV

2   Example 2

Consider the nonlinear non-Gaussian system shown in Eq.(8) and Eq.(9).where

The process noises are uncorrelated Gaussian white noises but the observation process noises are non-Gaussian noise with a mixed-Gaussian distribution. The simulation is performed many times for 1−500, where $v(u+1) \sim 0.8N(0,0.01) + 0.2N(0,0.1)$ for the first 200 times and $v(u + 1)\sim 0.8N (0,0.01) + 0.2N(0,0.1) +$${[\begin{array}{*{20}{c}} {0.1}&{0.1} \end{array}]^{\rm{T}}}$ for the remaining; ${w}(u) \sim N(0,0.01)$; all high-order polynomials are regarded as implicit values，and their modeling errors satisfies that the mean is 0 and the variance is 0.01; $x(0) = {[1\;\;1]^{\rm{T}}}$, $\hat x(0|0) = {[1\;\;1]^{\rm{T}}} +$$N(0,0.01) \times randn(2,1)$ and $P(0|0) = 0.01 \times diag(1,1)$

Case 2 only considers the measurement noise is non- Gaussian and takes $\sigma = 10,\varepsilon = {10^{ - 1}}$ as an example to compare the filtering performance between MCEKF and H-MCEKF, as shown in Fig.2(a)−(d)). Table 2 describes the error comparison in detail under several $\sigma$ and $\varepsilon$.

Figure  2.  The actual state $x_1$ and its estimate (a) with estimate error (c); The actual state $x_2$ and its estimate (b) with estimate error (d)

DownLoad: Full-Size Img PowerPoint

Table  2.  Estimation errors with different filters

		MSE of ${x_{\rm{1}}}$			MSE of ${x_{\rm{2}}}$			MSE
$\sigma$	$\varepsilon$	MCEKF	H-MCEKF	Improved	MCEKF	H-MCEKF	Improved	MCEKF	H-MCEKF	Improved
$\sigma {\rm{ = }}2$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.5800	0.1351	76.70%	0.0547	0.0363	33.64%	0.3173	0.0857	72.99%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.5246	0.2246	57.18%	0.0508	0.0412	18.89%	0.2877	0.1329	53.80%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.3086	0.2535	17.85%	0.1115	0.0453	57.56%	0.2101	0.1494	28.89%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.4088	0.2095	48.75%	0.0798	0.0617	22.68%	0.2443	0.1356	44.49%
$\sigma {\rm{ = }}5$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.1981	0.1950	1.56%	0.0299	0.0287	4.01%	0.1140	0.1119	1.82%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.2336	0.1843	21.10%	0.0341	0.0271	0.53%	0.1338	0.1057	21.00%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.1997	0.1804	9.66%	0.0461	0.0343	25.59%	0.1170	0.1133	3.16%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.2268	0.1910	15.78%	0.0461	0.0245	46.85%	0.1256	0.1185	5.65%
$\sigma {\rm{ = }}10$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.5157	0.1795	65.19%	0.0355	0.0272	23.38%	0.2756	0.1034	62.48%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.5888	0.1336	77.31%	0.0274	0.0273	0.36%	0.3081	0.0805	73.87%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.5517	0.1667	69.78%	0.0311	0.0260	16.39%	0.2914	0.0964	66.92%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.5862	0.2497	57.40%	0.0428	0.0409	4.44%	0.3145	0.1453	53.80%
$\sigma {\rm{ = }}15$	$\varepsilon {\rm{ = }}{10^{ - 1}}$	0.5326	0.2681	49.66%	0.0496	0.0461	7.06%	0.2911	0.1571	46.03%
	$\varepsilon {\rm{ = }}{10^{ - 2}}$	0.5294	0.2383	54.98%	0.0368	0.0331	10.05%	0.2831	0.1357	52.06%
	$\varepsilon {\rm{ = }}{10^{ - 4}}$	0.5461	0.1566	71.32%	0.0373	0.0346	7.23%	0.2917	0.0956	67.22%
	$\varepsilon {\rm{ = }}{10^{ - 6}}$	0.5109	0.2092	59.05%	0.0326	0.0212	34.96%	0.2718	0.1152	57.61%

 | Show Table

DownLoad: CSV

Case 1 and Case 2 take $\sigma = 10,\varepsilon = {10^{ - 1}}$ as examples to express the filtering performance of H-MCEKF and MCEKF, respectively. Table 1 and Table 2 summarize the corresponding MSEs of original state variables under different kernel bandwidth $\sigma$ and threshold $\varepsilon$ for H-MCEKF and MCEKF. In the simulation, $\sigma$ is set to 2, 5, 10 and 15 and $\varepsilon$ is set to ${10^{ - 1}}$, ${10^{ - 2}}$, ${10^{ - 4}}$, ${10^{ - 6}}$. The results confirms that the proposed H-MCEKF can outperform the MCEKF under the assumption that the state noise and measurement noise are non-Gaussian. Besides, we can see compared with the kernel bandwidth $\varepsilon$,$\sigma$ has little effect.

VII.   Conclusions and Prospects

A novel maximum correntropy high-order extended Kalman filter (H-MCEKF) has been designed for nonlinear and non-Gaussian systems in this paper. 1) the nonlinear polynomials has been defined as implicit function variables, which transforms the state model into pseudo-linearization. 2) establishing the linear model between all implicit function variables; 3) then the state model have been equivalently rewritten to linear model by combing with original states and implicit variables, which is similar to the measurement model.

Function variables with additive form can directly employ the new filter for state estimation. For general nonlinear functions, multi-dimensional Taylor network can be employed to develop into additive polynomials. But this method is not perfect, and model uncertainty still exists. Therefore, there is still a lot of work to solve nonlinear non-Gaussian systems.