Design Method for a Higher Order Extended Kalman Filter Based on Maximum Correlation Entropy and a Taylor Network System

Abstract

This paper proposes one new design method for a higher order extended Kalman filter based on combining maximum correlation entropy with a Taylor network system to create a nonlinear random dynamic system with modeling errors and unknown statistical properties. Firstly, the transfer function and measurement function are transformed into a nonlinear random dynamic model with a polynomial form via system identification through the multidimensional Taylor network. Secondly, the higher order polynomials in the transformed state model and measurement model are defined as implicit variables of the system. At the same time, the state model and the measurement model are equivalent to the pseudolinear model based on the combination of the original variable and the hidden variable. Thirdly, higher order hidden variables are treated as additive parameters of the system; then, we establish an extended dimensional linear state model and a measurement model combining state and parameters via the previously used random dynamic model. Finally, as we only know the results of the limited sampling of the random modeling error, we use the combination of the maximum correlation estimator and the Kalman filter to establish a new higher order extended Kalman filter. The effectiveness of the new filter is verified by digital simulation.

1. Introduction

The application of filters occupies an important position in various fields at the national and international levels. The progress and development of filters play important roles in national economic construction—especially national defense construction—such as in real-time estimation and target tracking. In 1960, Kalman proposed a method of filtering under the minimum mean squared error criterion for linear systems, and it soon began to be widely used [1]. In order to solve nonlinear problems, extended Kalman filters (EKFs) [2], unscented Kalman filters (UKFs) [3], and cubature Kalman filters (CKFs) have since emerged. However, the above-mentioned filtering methods require the modeling error to be Gaussian white noise. As such, their performances are likely to worsen when applied to non-Gaussian situations, especially when the systems are disturbed by impulsive noise. Impulsive noise arises from heavy-tailed distributions [4] (such as some mixed Gaussian distributions), and is common in many real scenarios of automatic control and target tracking (for instance, the measurement noise in the radar system is often not Gaussian noise, but heavy-tailed non-Gaussian noise [5]). In 1993, Gordon and Salmond proposed particle filtering when the density function is known [6]; this achieves an approximation of the distribution function by sampling a large number of particles therein. However, this method is very complicated; it requires a large number of particles, and it will cause particle degradation after re-sampling. In general, the density function is difficult to obtain. For this reason, for the linear system, Chen designed the corresponding Kalman filter under the maximum correlation entropy criterion based on the limited realization of random variables [7]; this is called the maximum correlation entropy Kalman filter (MCKF) [8]. On this basis, the maximum correntropy extended Kalman filter (MCEKF) and the maximum correntropy unscented Kalman filter (MCUKF), which can solve nonlinear non-Gaussian systems, have since emerged [9]. However, in MCEKFs, all higher order terms in the Taylor expansion are discarded. Therefore, a large truncation error will be generated, and the filtering performance will decrease or even diverge as the nonlinearity of the system increases. In addition, each step of the state estimation needs to recalculate the Taylor expansion coefficient, which will undoubtedly increase the complexity of the calculation. MCUKFs use UT transformation and sigma point sampling [10]; this is called deterministic sampling. There is only one sampling point for a dimensional system. Neither low-dimensional nor high-dimensional systems have a strong claim to superiority. A large number of experiments have shown that both EKFs and UKFs can be approximated by a second-order polynomial at most [11], which will produce a large rounding error. Hence, both will eventually face the problems of degraded filtering performance and divergence as their nonlinearity increases [12].

This project proposes a higher order extended Kalman filter method based on maximum correlation entropy, under the assumption that both state and measurement equations can be modeled and based on a strong nonlinear function. The main contributions of this paper are as follows: (1) using multidimensional Taylor nets to convert the general expression of nonlinear functions into higher order polynomials; (2) defining each order of polynomial in the system as hidden variables of the corresponding order, and treating them as time-variable parameters; (3) establishing the dynamic relationship between the time-variable parameters and combining them with the original variables to further establish the expanded dimension state model; (4) based on the expanded linear state variables, equivalently rewriting the measurement model into the corresponding linear form; and (5) according to the established linear state and measurement model of the new extended dimension system, establishing a higher order extended Kalman filter method based on maximum correlation entropy.

The remaining parts of this paper are organized as follows: the first chapter is the preface of our knowledge, which introduces the definition of “entropy”; the Section 2 presents a method for identifying nonlinear functions based on multidimensional Taylor networks; the Section 3 presents a higher order extended Kalman filter method; the Section 4 presents the detailed design process of the maximum correlation entropy higher order extended Kalman filter; the Section 5 concerns simulation verification; and the Section 6 and Section 7 presents a summary and outlook.

2. Description of Correntropy

Correntropy is a generalized similarity measure between two random variables [13]. Given two one-dimensional random variables $\phi ,\zeta \in R^1$, their joint distribution function is $F_{\varphi \xi }(\phi ,\zeta )$; then, the correlation entropy is defined as follows:

$$V(\varphi ,\xi )=\epsilon \left[\alpha (\phi ,\zeta )\right]={\displaystyle \int \alpha (\phi ,\zeta )d}{F}_{\varphi \xi }(\phi ,\zeta )$$

(1)

where $\epsilon$ is the expectation operator and $\alpha (\cdot ,\cdot )$ is the translation-invariant Mercer kernel. In this article, it is not particularly emphasized that this kernel function is a Gaussian kernel, which is defined as follows:

$$\alpha \left(\phi ,\zeta \right)=G_{\tau }(e)=exp\left(-\frac{e^2}{2{\tau }^2}\right)$$

(2)

where $e=\phi -\zeta$, $\tau >0$ represents the kernel’s bandwidth.

By expanding Equation (2) with a Taylor series, we can obtain the following:

$$\alpha (\phi ,\zeta )=G_{\tau }(e)=exp\left(-\frac{e^2}{2{\tau }^2}\right)={\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}\epsilon \left\{{(\phi -\zeta )}^{2k}\right\}$$

(3)

and then the correlation entropy of Equation (1) has the following expression:

$$\begin{array}{l}V(\varphi ,\xi )=\epsilon [\alpha (\phi ,\zeta )]={\displaystyle \int {\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}}{(\phi -\zeta )}^{2k}d{F}_{\varphi \xi }(\phi ,\zeta )\\ \hspace{1em}\hspace{1em} ={\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}{\displaystyle \int {(\phi -\zeta )}^{2k}}d{F}_{\varphi \xi }(\phi ,\zeta )={\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}\epsilon \left\{{(\varphi -\xi )}^{2k}\right\}\end{array}$$

(4)

where $\epsilon \left\{{\left(\varphi -\xi \right)}^{2k}\right\}={\int }^{ }{\left(\phi -\zeta \right)}^{2k}d{F}_{\varphi \xi }\left(\phi ,\zeta \right)$ is the $2k$ truncation statistic of the random variable $\varphi ,\xi \in R$.

However, in most practical cases, joint distribution $F_{\varphi \xi }$ is usually unknown, and there are often finite implementations $\left({\phi }^{\left(j\right)},{\zeta }^{\left(j\right)}\right),j=1,2,\cdots ,N$ of $\left(\varphi ,\xi \right)$ for random variables. In these cases, the sample mean estimator can be used to estimate the heterogeneity:

$$\epsilon \left\{{(\varphi -\xi )}^{2k}\right\}=\frac{1}{N}\left({\displaystyle \sum _{j=0}^N{\left({\phi }^{(j)}-{\zeta }^{(j)}\right)}^{2k}}\right)$$

(5)

Then, the entropy expression of the random variable pair $\left(\phi , \zeta \right)$ is driven by finite data:

$$\begin{array}{l}\widehat{V}(\varphi ,\xi )=\epsilon [\alpha (\phi ,\zeta )]\\ \hspace{1em}\hspace{1em} ={\displaystyle \int {\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}}{(\phi -\zeta )}^{2k}d{F}_{\varphi \xi }(\phi ,\zeta )\\ \hspace{1em}\hspace{1em} ={\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}\frac{1}{N}\left({\displaystyle \sum _{j=0}^N{\left({\phi }^{(j)}-{\zeta }^{(j)}\right)}^{2k}}\right)\\ \hspace{1em}\hspace{1em} =\frac{1}{N}{\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}\left({\displaystyle \sum _{j=0}^N{\left({\phi }^{(j)}-{\zeta }^{(j)}\right)}^{2k}}\right)\\ \hspace{1em}\hspace{1em} =\frac{1}{N}{\displaystyle \sum _{k=0}^{\infty }\frac{{(-1)}^k}{2^k{\tau }^{2k}k!}}\left({\displaystyle \sum _{j=0}^N{\left(e^{(j)}\right)}^{2k}}\right)=\frac{1}{N}{\displaystyle \sum _{j=0}^N{G}_{\tau }}\left(e^{(j)}\right)\end{array}$$

(6)

When $\varphi ,\xi \in R^n$, and the components of vector $e=\varphi -\xi$ are independent of one another, multidimensional correlation entropy is based on N sampling.

3. Non-Linear Model Identification Based on Multidimensional Taylor Networks

Given that the state model and observation model are complex dynamic systems with nonlinear characteristics [14]:

$$\alpha (\tau +1)=\sigma (\alpha (\tau ))+\gamma (\tau )$$

(7)

$$\mathsf{\Gamma }(\tau +1)=\delta (\alpha (\tau +1))+\theta (\tau +1)$$

(8)

where $\alpha (\tau )\in R^h$ is an h-dimensional state vector; $\mathrm{T}(\tau +1)\in R^d$ represents the d-dimensional measure vector; and ${\sigma }_i(\alpha (\tau )), i=1,2,\cdots ,h$ and ${\delta }_j(\alpha (\tau +1)), j=1,2,\cdots ,d$ represent the state function and the measurement function, respectively. The modeling errors for non-Gaussian systems are $\gamma (k)$ and $\theta (k)$, while $\vartheta =diag\{{\vartheta }_1,{\vartheta }_2,\cdots ,{\vartheta }_{\mathrm{h}}\}$ and $\eta =diag\left(\begin{array}{cccc}{\eta }_1& {\eta }_2& \cdots & {\eta }_d\end{array}\right)$ are the process noise variance and the measurement noise variance, respectively.

Lemma 1.

Any continuous function defined in a closed interval can be approximated accurately with a polynomial function [15].

Lemma 2.

For continuous functions, $\sigma (\alpha (k))$, defined in a closed interval, can be approximated by the following [16]:

$${\displaystyle \sum _{i=1}^{N(h,l)}{\psi }_i}(k){\displaystyle \prod _{t=1}^l{\alpha }_t^{{\lambda }_{i,t}}}(k)$$

(9)

where $N\left(h,l\right)$ denotes the total number of terms in the expansion and ${\lambda }_{i,t}$ denotes the power of the variable ${\alpha }_t$ in the product of the ith variable.

3.1. Multidimensional Taylor Network Structure

The multidimensional Taylor network model can replace the traditional neural network with the dynamic model and control the system under certain conditions; it is characterized by a nonlinear autoregressive moving-average model composed of polynomials. The multidimensional Taylor network (MTN) uses a forward single intermediate layer structure, including an input layer, an intermediate layer, and an output layer. Supposing that the input layer comprises n nodes—$\alpha (\tau )={\left[\begin{array}{cccc}{\alpha }_1(\tau )& {\alpha }_2(\tau )& \cdots & {\alpha }_h(\tau )\end{array}\right]}^T\in R^h$—the output layer is $\alpha (\tau +1)$, the middle layer is the network processing layer, and each input variable realizes the weighted summation of each power product term in this layer. The middle layer is composed of various power product terms and the corresponding connection weight vector ${\psi }_j(\tau )$:

$${\psi }_j(\tau )={\left[{\psi }_{j,1}(\tau ),{\psi }_{j,2}(\tau ),\cdots ,{\psi }_{j,N(h,l)}(\tau )\right]}^T$$

which represents the output weight vector connecting the intermediate layer and the output node of the network.

According to the multivariate Taylor equation, if a function is differentiable to the $h+1$th order at a certain point, then the function expands to a form where the power series of the variable is not greater than m times. The model can be expressed as a dynamic equation, as follows:

$${\alpha }_j(\tau +1)=\sigma (\alpha (\tau ))={\displaystyle \sum _{i=1}^{N(h,l)}{\psi }_{j,i}}(\tau ){\displaystyle \prod _{t=1}^l{\alpha }_t^{{\lambda }_{i,t}}}(\tau )+\mathsf{\Delta }\sigma (\tau )$$

(10)

where $\sigma (\cdot )$ is a function of nonlinearity described by a multidimensional Taylor network model, ${\psi }_i$ represents the weight before the product item of the ith variable, $N\left(h,l\right)$ denotes the total number of terms in the expansion, ${\lambda }_{i,t}$ denotes the power of the variable ${\alpha }_t$ in the product of the ith variable, and $\mathsf{\Delta }\sigma (\tau )$ is the error—also known as the remainder—produced by the identification of a function by a multidimensional Taylor network.

3.2. Parameter Identification Method Based on Kalman Filtering

Model Establishment of a Kalman Filter

A Kalman filter can be regarded as an optimized autoregressive data processing method that describes the entire system through a state equation and an observation equation.

State equation:

$${\psi }_{j,i}(\tau +1)={\psi }_{j,i}(\tau )+w_{j,i}(\tau )$$

(11)

where $i=1,2,\cdots ,N(h,l)$, $j=1,2\cdots ,h$.

Observation equation:

It is not difficult to draw from Figure 1:

$$\begin{array}{l}{\alpha }_j(\tau +1)={\displaystyle \sum _{i=1}^{N(h,l)}{\psi }_{j,i}}(\tau +1){\displaystyle \prod _{t=i}^l{\alpha }_t^{{\lambda }_{i,t}}}(\tau +1)\\ \hspace{1em}\hspace{1em}\hspace{1em} =H_j(\tau +1)\cdot {\left[{\psi }_{j,1}(\tau +1),{\psi }_{j,2}(\tau +1),\cdots ,{\psi }_{j,N(h,l)}(\tau +1)\right]}^{\mathrm{T}}+v_j(\tau +1)\\ \hspace{1em}\hspace{1em}\hspace{1em} =H_j(\tau +1)\cdot {\psi }_j(\tau +1)+v_j(\tau +1)\end{array}$$

Thus,

$$\begin{array}{l}\alpha (\tau +1)={\left[{\alpha }_1(\tau +1),{\alpha }_2(\tau +1),\cdots ,{\alpha }_j(\tau +1),\cdots ,{\alpha }_h(\tau +1)\right]}^T\\ \hspace{1em}\hspace{1em}\hspace{1em} =H(\tau )\cdot \psi (\tau +1)+v(\tau +1)\end{array}$$

(12)

where $H(\tau )={\left[H_1(\tau +1),H_2(\tau +1),\cdots ,H_h(\tau +1)\right]}^T$; $\psi (\tau +1)=\left[{\psi }_1(\tau +1),{\psi }_2(\tau +1),\cdots ,{\psi }_h(\tau +1)\right]$; ${\psi }_{j,i}(\tau +1)$ represents the system state at $\tau$-time, that is, the parameter status value of the kth moment; and $\beta (\tau +1)$ represents the output value of the network. It is assumed that both process noise $w(\tau )$ and $v(\tau +1)$ are Gaussian white noise during the analysis, and $Q_j=diag\left(\begin{array}{cccc}{Q}_{j,1}& Q_{j,2}& \cdots & Q_{j,N(h,l)}\end{array}\right)$ and $R_j=diag\left(\begin{array}{cccc}{R}_{j,1}& R_{j,2}& \cdots & R_{j,N(h,l)}\end{array}\right)$, which are the process noise variance and measurement noise variance, respectively. Here, we use a Kalman filter to approximate the dynamic model. As the filtering principle of Kalman filters is mentioned later in this article, please refer to Equations (20)–(24) for the detailed process.

Figure 1. Model of a multidimensional Taylor network.

Figure 1. Model of a multidimensional Taylor network.

3.3. Approximation Analysis

Given a class of nonlinear functions $\sigma (\alpha (k))$, it can be assumed that it is derivative of the rth order, but r is a relatively large number, making it difficult for us to use Taylor nets to approximate its function. The optimal approach would be to set $m,1\le m\le r$ and use the Taylor network to expand the nonlinear function to the mth order, obtain the result of Equation (16), and simultaneously ensure the higher order error term $\mathsf{\Delta }\delta \le \theta$, where $\theta$ is the acceptable error threshold. This not only makes the Taylor network fitting function process easier, but also ensures the accuracy of the fit.

4. Higher Order Extended Kalman Filter

4.1. Pseudolinearized Representation of Nonlinear Functions

For ease of description and understanding, if $l=d=2$, we can expand Equation (7) through a multidimensional Taylor network to the mth order, as follows:

$$\begin{array}{l}{\sigma }_i(\alpha (k))=\left({\omega }_{i,1,0}{\alpha }_1(\tau )+{\omega }_{i,0,1}{\alpha }_2(\tau )\right)\\ \begin{array}{l}\hspace{1em}\hspace{1em} +\left({\omega }_{i,2,0}{\alpha }_1{(\tau )}^2+{\omega }_{i,1,1}{\alpha }_1(\tau ){\alpha }_2(\tau )+{\omega }_{i,0,2}{\alpha }_2{(\tau )}^2\right)\\ \hspace{1em}\hspace{1em} +\left({\omega }_{i,3,0}{\alpha }_1{(\tau )}^3+{\omega }_{i,1,2}{\alpha }_1(\tau ){\alpha }_2^2(\tau )+{\omega }_{i,2,1}{\alpha }_1^2(\tau ){\alpha }_2(\tau )+{\omega }_{i,0,3}{\alpha }_2{(\tau )}^3\right)\\ \hspace{1em}\hspace{1em} +{\displaystyle \sum _{\begin{array}{c}{l}_1+l_2=l\\ l_1+l_2\le l\end{array}}{\omega }_{i,l_1,l_2}{\alpha }_1^{l_1}(\tau ){\alpha }_2^{l_2}(\tau )+\cdots +{\displaystyle \sum _{\begin{array}{c}{m}_1+m_2=r\\ m_1+m_2\le r\end{array}}{\omega }_{i,r_1,r_2}{\alpha }_1^{m_1}(\tau ){\alpha }_2^{m_2}(\tau )+\mathsf{\Delta }{\sigma }_i(\tau )}}\end{array}\end{array}$$

(13)

where ${\displaystyle \sum _{\begin{array}{l}{l}_1+l_2=l\\ l_1,l_2\le l\end{array}}{\alpha }_1^{l_1}(\tau ){\alpha }_2^{l_2}(\tau )}$ is the sum of all tensors of the $l$ th order and ${\omega }_{i,l_1,l_2}$ represents the weight corresponding to each order of the tensor.

Definition 1.

${\alpha }^{(l)}(\tau )=\{{\alpha }_1^{l_1}(\tau ){\alpha }_2^{l_2}(\tau )\cdots {\alpha }_2^{l_h}(\tau ), l_1+l_2+l_h=l; 0\le l_j\le l; l=0,1,\cdots ,h\}$ is a set of implicit variables of the$l$th order.

Definition 2.

${\omega }_i^{(l)}=\left[{\omega }_{i; 1}^{(l)}, {\omega }_{i; 2}^{(l)}, \cdots ,{\omega }_{i; n_l}^{(l)}\right]=\left[{\omega }_{i; l,0}, {\omega }_{i; l-1,1},\cdots , {\omega }_{i; 0,l}\right], i=1,2,\cdots l$ is the weight vector corresponding to the ith order implicit variable.

In [17], there is a detailed pseudolinearization process, so we will not repeat it in this article. In order to make the model more accurate, we treat the remainder $\mathsf{\Delta }\sigma (\tau )$ of the equation of state as latent variables. According to Definition 1 and Definition 2, the pseudolinear extended dimension form using the remainder as a hidden variable is as follows:

$${\alpha }^{(1)}(\tau +1)=W^{(1)}(\tau +1,\tau ){\alpha }^{(1)}(\tau )+{\displaystyle \sum _{l=2}^m{W}^{(l)}(\tau +1,\tau ){\alpha }^{(l)}(\tau )}+C\cdot \mathsf{\Delta }\sigma (\tau )+{\gamma }^{(1)}(\tau )$$

(14)

where ${\alpha }^{(1)}(\tau )=\left[\begin{array}{c}{\alpha }_1^{(1)}(\tau )\\ {\alpha }_2^{(1)}(\tau )\end{array}\right]$, $W^{(l)}=\left[\begin{array}{c}{\omega }_1^{(l)}\\ {\omega }_2^{(l)}\end{array}\right]$, $\gamma (\tau )={\gamma }^{(1)}(\tau )=\left[\begin{array}{c}{\gamma }_1^{(1)}(\tau )\\ {\gamma }_2^{(1)}(\tau )\end{array}\right]$, $C=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$.

Similarly, Equation (8) can be rewritten as follows:

$${\mathsf{\Gamma }}^{(1)}(\tau +1)={\chi }^{(1)}(\tau +1){\alpha }^{(1)}(\tau +1)+{\displaystyle \sum _{l=2}^m{\chi }^{(l)}(\tau +1){\alpha }^{(l)}(\tau +1)}+D\cdot \mathsf{\Delta }\sigma (\tau +1)+{\theta }^{(1)}(\tau +1)$$

(15)

where ${\mathsf{\Gamma }}^{(1)}(\tau +1)=\left[\begin{array}{c}{\mathsf{\Gamma }}_1^{(1)}(\tau )\\ {\mathsf{\Gamma }}_2^{(1)}(\tau )\end{array}\right]$, ${\chi }^{(l)}=\left[\begin{array}{c}{\chi }_1^{(l)}\\ {\chi }_2^{(l)}\end{array}\right]$, ${\theta }^{(1)}(\tau +1)=\left[\begin{array}{c}{\theta }_1^{(1)}(\tau +1)\\ {\theta }_2^{(1)}(\tau +1)\end{array}\right]$, $D=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right]$.

4.2. Linearized Representation of Nonlinear Functions

In order to transform the pseudolinear model established in Section 3.1 into a true linear form, it is necessary to establish a dynamic relationship between the lth order hidden variables and the uth order hidden variables [18]:

$${\alpha }^{(l)}(\tau +1)=W_l^{(u)}(\tau ){\alpha }^{(u)}(\tau ) l,u=2,3,\cdots ,m$$

(16)

where $W$ can be identified based on the multidimensional Taylor network in its original state; without any prior information, it can be set as follows:

$$W_l^{(u)}(\tau )=\{\begin{array}{l}I, l=u\hfill \\ 0, l\ne u\hfill \end{array}$$

(17)

Combining Definition 1, Definition 2, and Equation (19), the state model Equation (7) has the following linear matrix form:

If $A(\tau )={\left[{({\alpha }^{(1)}(\tau ))}^T, {({\alpha }^{(2)}(\tau ))}^T, \cdots , {({\alpha }^{(l)}(\tau ))}^T, \cdots , {({\alpha }^{(r)}(\tau ))}^T,\mathsf{\Delta }\sigma (\tau )\right]}^T$

$$W(\tau +1,\tau )=\left[\begin{array}{cccccccc}{W}_1^{(1)}(\tau )& W_1^{(2)}(\tau )& \cdots & W_1^{(u)}(\tau )& \cdots & W_1^{(m-1)}(\tau )& W_1^{(m)}(\tau )& C\\ W_2^{(1)}(\tau )& W_2^{(2)}(\tau )& \cdots & W_2^{(u)}(\tau )& \cdots & W_2^{(m-1)}(\tau )& W_2^{(m)}(\tau )& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮& ⋮\\ W_l^{(1)}(\tau )& W_l^{(2)}(\tau )& \cdots & W_l^{(u)}(\tau )& \cdots & W_l^{(m-1)}(\tau )& W_l^{(m)}(\tau )& 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮& ⋮& ⋮\\ W_{m-1}^{(1)}(\tau )& W_{m-1}^{(2)}(\tau )& \cdots & W_{m-1}^{(u)}(\tau )& \cdots & W_{m-1}^{(m-1)}(\tau )& W_{m-1}^{(m)}(\tau )& 0\\ W_m^{(1)}(\tau )& W_m^{(2)}(\tau )& \cdots & W_m^{(u)}(\tau )& \cdots & W_m^{(m-1)}(\tau )& W_m^{(m)}(\tau )& 0\\ 0& 0& \cdots & 0& \cdots & 0& 0& C\end{array}\right], \gamma (\tau )=\left[\begin{array}{c}{\gamma }^{(1)}(\tau )\\ {\gamma }^{(2)}(\tau )\\ ⋮\\ {\gamma }^{(l)}(\tau )\\ ⋮\\ {\gamma }^{(m-1)}(\tau )\\ {\gamma }^{(m)}(\tau )\end{array}\right]$$

then, Equation (7) has the following linearized form:

$$A(\tau +1)=W(\tau +1,\tau )A(\tau )+\gamma (\tau )$$

(18)

where $\gamma (k)$ is the modeling error.

In the same way, the linear matrix form of the measurement model can be obtained:

$$\mathsf{\Gamma }(\tau +1)=\chi (\tau +1,\tau )A(\tau +1)+\theta (\tau +1)$$

(19)

where $\chi (\tau +1,\tau )=\left[\begin{array}{ccccccccc}{\chi }_1^{(1)}(\tau +1)& {\chi }_1^{(2)}(\tau +1)& \cdots & {\chi }_1^{(u)}(\tau +1)& \cdots & {\chi }_1^{(m-1)}(\tau +1)& {\chi }_1^{(m)}(\tau +1)& 0& 0\\ {\chi }_2^{(1)}(\tau +1)& {\chi }_2^{(2)}(\tau +1)& \cdots & {\chi }_2^{(u)}(\tau +1)& \cdots & {\chi }_2^{(m-1)}(\tau +1)& {\chi }_2^{(m)}(\tau +1)& 0& 0\end{array}\right]$, $\mathsf{\Gamma }(\tau +1)=\left[\begin{array}{c}{\mathsf{\Gamma }}_1(\tau )\\ {\mathsf{\Gamma }}_2(\tau )\end{array}\right]$, and $\theta (\tau +1)=\left[\begin{array}{c}{\theta }_1(\tau +1)\\ {\theta }_2(\tau +1)\end{array}\right]$ is the modeling error.

4.3. Design of Higher Order Extended Kalman Filter

For linear models, KF-based filters are given. Given the initial value $A(0)$, when $\gamma (\tau )$ and $\theta (\tau +1)$ are Gaussian white noise with zero mean, the variances are recorded as $\vartheta$ and $\eta$, respectively.

A recursive filter can be designed as follows:

$$\widehat{A}(\tau +1|\tau )=W(\tau +1,\tau )\widehat{A}(\tau |\tau )$$

(20)

$$\lambda (\tau +1|\tau )=W(\tau +1,\tau )\lambda (\tau |\tau )W^T(\tau +1,\tau )+\vartheta (\tau )$$

(21)

$$\mathsf{{\rm K}}(\tau +1)=(\lambda (\tau +1|\tau ){\chi }^T(\tau +1)){(\chi (\tau +1)\lambda (\tau +1|\tau ){\chi }^T(\tau +1)+\eta (\tau +1))}^{-1}$$

(22)

$$\widehat{A}(\tau +1|\tau +1)=\widehat{A}(\tau +1|\tau )+\mathsf{{\rm K}}(\tau +1)(\mathsf{\Gamma }(\tau +1)-\chi (\tau +1)\widehat{A}(\tau +1|\tau ))$$

(23)

$$\lambda (\tau +1|\tau +1)=(I-\mathsf{{\rm K}}(\tau +1)\chi (\tau +1))\lambda (\tau +1|\tau )$$

(24)

5. Higher Order Extended Kalman Filter Design Based on Maximum Correlation Entropy

5.1. Non-Gaussian Modeling of State Vector Based on Multivariate Information Observation

System status $A(\tau )$ estimates $\widehat{A}(\tau |\tau )$ and estimated error covariance $\lambda (\tau |\tau )$ are obtained based on $\mathsf{{\rm K}}$ [19]. The filtering equation predicts a step prediction estimation value $\tilde{A}(\tau +1|\tau )$ and corresponding step prediction error covariance matrix $\lambda (k+1|k)$.

The step prediction estimate error of system status $A(\tau +1)$ is as follows:

$$\begin{array}{l}\tilde{A}(\tau +1|\tau )=A(\tau +1)-\widehat{A}(\tau +1|\tau )\\ \hspace{1em} =W(\tau +1,\tau )\tilde{A}(\tau |\tau )+\gamma (\tau )\end{array}$$

(25)

and it can be modified into a measurement model for system status $A(\tau +1)$ as follows:

$$\widehat{A}(\tau +1|\tau )=A(\tau +1)-\tilde{A}(\tau +1|\tau )$$

(26)

where $\widehat{A}(\tau +1|\tau )$ is a measurement of system status $A(\tau +1)$, while $\tilde{A}(\tau +1|\tau )$ is the measurement error. Finally, the combined measurement model is as follows:

$$\left[\begin{array}{c}\widehat{A}(\tau +1|\tau )\\ \mathsf{\Gamma }(\tau +1)\end{array}\right]=\left[\begin{array}{c}I\\ \chi (\tau +1)\end{array}\right]A(\tau +1)+{\varpi }^{(j)}(\tau +1)$$

(27)

where I is a unit array for the corresponding dimension, ${\varpi }^{(j)}(\tau +1)=\left[\begin{array}{c}-\tilde{A}(\tau +1|\tau ))\\ {\theta }^{(j)}(\tau +1)\end{array}\right]$, and

$$E\left[{\varpi }^{(i)}(\tau +1){({\varpi }^{(i)})}^T(\tau +1)\right]=\left[\begin{array}{cc}\tilde{\lambda }(\tau +1|\tau )& 0\\ 0& \tilde{\eta }(\tau +1)\end{array}\right]$$

(28)

According to Equation (20), a step prediction error covariance of the system state $\lambda (\tau +1)$ is received as follows:

$$\tilde{\lambda }(\tau +1|\tau )=W(\tau +1,\tau )\lambda (\tau |\tau )W^T(\tau +1,\tau )+\vartheta (\tau )$$

(29)

where $\vartheta (\tau )=diag\{{\vartheta }^{(1)}(\tau ),{\vartheta }^{(2)}(\tau ),\cdots ,{\vartheta }^{(m)}(\tau )\}$ and the ${\vartheta }^{(2)}(\tau ),\cdots ,{\vartheta }^{(m)}(\tau )$ is a covariance matrix of random error vectors ${\gamma }^{(2)}(\tau ),\cdots ,{\gamma }^{(r)}(\tau )$ when the higher order hidden variable ${\alpha }^{(2)}(\tau ),\cdots ,{\alpha }^{(m)}(\tau )$ is dynamically modeled. ${\vartheta }^{(1)}(\tau )$ is the original system status model (Equation (16)) of the non-Gaussian model error $w^{(1)}(k)$, and calculates the second-order statistic after obtaining a limited number of samples:

$${\tilde{\vartheta }}^{(1)}(\tau )=\frac{1}{N}{\displaystyle \sum _{j=1}^N\left\{[{\gamma }^{(1,j)}(\tau )-\overline{\gamma }(\tau )]{[{\gamma }^{(1,j)}(\tau )-\overline{\gamma }(\tau )]}^T\right\}}$$

(30)

In Equation (23), $\tilde{\eta }(\tau +1)$ is the calculated second-order statistic calculated after the non-Gaussian model error $\theta (\tau +1)$ in the original system measurement state model (Equation (8)), obtaining a limited sample:

$$\tilde{\eta }(\tau +1)=\frac{1}{N}{\displaystyle \sum _{j=1}^N\{[{\theta }^{(j)}(\tau +1)-\overline{\theta }(\tau +1)]{[{\theta }^{(j)}(\tau +1)-\overline{\theta }(\tau +1)]}^T\}}$$

(31)

where ${\theta }^{(j)}(\tau +1)$ is the jth realization vector of the non-Gaussian random noise vector $\theta (\tau +1)$.

5.2. The Statistical Independence Process of Each Component in the Non-Gaussian Modeling Error Vector $\varpi (\tau +1)$ in the Comprehensive Measurement Model

The vector $\varpi (\tau +1)$ in the comprehensive measurement model Equation (22) is a non-Gaussian modeling error vector, and its components are not statistically independent. In order to use the correlation entropy form of the multidimensional independent vector shown in Equation (19), the one-dimensional non-Gaussian vector $\varpi (\tau +1)$ needs to be transformed into statistical independence.

From $\lambda (\tau +1|\tau )=E\left\{\left[A(\tau +1|\tau )-A(\tau |\tau )\right]{\left[A(\tau +1|\tau )-A(\tau |\tau )\right]}^T\right\}$, $\lambda (\tau +1|\tau )$ is a positive definite matrix. Similarly, in Equation (26), $\tilde{\eta }(\tau +1)$ is also a positive definite matrix. For this reason, Equation (23) is further expressed as follows:

$$\begin{array}{l}E\left\{{\varpi }^{(i)}(\tau +1){({\varpi }^{(i)})}^T(\tau +1)\right\}=\left[\begin{array}{cc}{\mathsf{\Lambda }}_{\alpha }(\tau +1|\tau ){\mathsf{\Lambda }}_{\alpha }^T(\tau +1|\tau )& 0\\ 0& {\mathsf{\Lambda }}_{\mathsf{\Gamma }}(\tau +1){\mathsf{\Lambda }}_{\mathsf{\Gamma }}^T(\tau +1)\end{array}\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\mathsf{\Lambda }(\tau +1){\mathsf{\Lambda }}^T(\tau +1)\end{array}$$

(32)

where ${\mathsf{\Lambda }}_{\alpha }(\tau +1)$ and ${\mathsf{\Lambda }}_{\mathsf{\Gamma }}(\tau +1)$ are the Cholesky factor matrices of $\tilde{\lambda }(\tau +1|\tau )$ and $\tilde{\eta }(\tau +1)$, respectively.

Applying ${\mathsf{\Lambda }}^{-1}(\tau +1)$ to both sides of Equation (22), respectively, yields:

$${\mathsf{\Lambda }}^{-1}(\tau +1)\left[\begin{array}{c}\widehat{A}(\tau +1|\tau )\\ \mathsf{\Gamma }(\tau +1)\end{array}\right]={\mathsf{\Lambda }}^{-1}(\tau +1)\left[\begin{array}{c}I\\ \chi (\tau +1)\end{array}\right]A(\tau +1)+{\mathsf{\Lambda }}^{-1}(\tau +1){\varpi }^{(j)}(\tau +1)$$

(33)

where

$$D(\tau +1)={\mathsf{\Lambda }}^{-1}(\tau +1)\left[\begin{array}{c}\widehat{A}(\tau +1|\tau )\\ \mathrm{T}(\tau +1)\end{array}\right],S(\tau +1)={\mathsf{\Lambda }}^{-1}(\tau +1)\left[\begin{array}{c}I\\ \chi (\tau +1)\end{array}\right],e(\tau +1)={\mathsf{\Lambda }}^{-1}(\tau +1)\varpi (\tau +1)$$

The above equation can be further simplified as follows:

$$D(\tau +1)=S(\tau +1)A(\tau +1)+e(\tau +1)$$

(34)

because

$$\begin{array}{l}E\left\{e(\tau )e^T(\tau )\right\}=E\left\{\left[{\mathsf{\Lambda }}^{-1}(\tau +1)\varpi (\tau +1)\right]{\left[{\mathsf{\Lambda }}^{-1}(\tau +1)\varpi (\tau +1)\right]}^T\right\}\\ \hspace{1em} ={\mathsf{\Lambda }}^{-1}(\tau +1)E\left\{\varpi (\tau +1){\varpi }^T(\tau +1)\right\}{({\mathsf{\Lambda }}^{-1}(\tau +1))}^T\\ \hspace{1em} ={\mathsf{\Lambda }}^{-1}(\tau +1)\mathsf{\Lambda }(\tau +1){\mathsf{\Lambda }}^T(\tau +1){({\mathsf{\Lambda }}^{-1}(\tau +1))}^T\\ \hspace{1em} =I\end{array}$$

(35)

Therefore, after the non-Gaussian modeling error random variable $\varpi (\tau +1)$ undergoes the equivalent transformation of the matrix ${\mathsf{\Lambda }}^{-1}(\tau +1)$, the components of the random $e(\tau +1)$ are statistically independent.

5.3. Implementation Process of a Higher Order Extended Kalman Filter Based on Maximum Entropy

The filtering process of the extended Kalman filter (H-MCEKF) based on the maximum correlation entropy is as follows (see [20] for the specific derivation process):

• The filter initialization obtains the initial filter value $\widehat{A}(0)$ and the covariance $\lambda (0)$, choosing a suitable core bandwidth $o$ and a small positive number $\epsilon$;
• Taylor networks are used for system identification to obtain the parameters in the equations, using the expanded item and the remainder as the new hidden variables. A pseudolinearization process is performed to obtain the pseudolinear form of the system;
• Equations (20) and (21) are used to obtain $\widehat{X}(k+1|k)$ and $P(k+1|k)$, respectively, while Cholesky decomposition is used to obtain $B_p(k+1|k)$;
• $t=1$ and $\widehat{A}{(\tau +1|\tau +1)}_0=\widehat{A}(\tau +1|\tau )$ are taken, where $\widehat{A}{(\tau +1|\tau +1)}_t$ represents the estimated state of the fixed-point iteration t;
• The starting fixed-point iterative algorithm is as follows:

  $${\tilde{e}}_i\left(\tau +1\right)=d_i\left(\tau +1\right)-s_i\left(\tau +1\right)\widehat{A}{(\tau +1|\tau +1)}_{t-1}$$

  (36)

  where $e_i$ is the ith element of $e$:

  $${\tilde{C}}_{\alpha }(\tau +1)=diag\left(G_o\left({\tilde{e}}_1(\tau +1)\right),\cdots ,G_o\left({\tilde{e}}_m(\tau +1)\right)\right)$$

  (37)

  $${\tilde{C}}_{\mathsf{\Gamma }}(\tau +1)=diag\left(G_o\left({\tilde{e}}_{r+1}(\tau +1)\right),\cdots ,G_o\left({\tilde{e}}_{m+m}(\tau +1)\right)\right)$$

  (38)

  $$\tilde{H}(\tau +1)={\mathsf{\Lambda }}_{\mathsf{\Gamma }}(\tau +1){\tilde{C}}_{\mathsf{\Gamma }}^{-1}(\tau +1){\mathsf{\Lambda }}_{\mathsf{\Gamma }}^T(\tau +1)$$

  (39)

  $$\tilde{\lambda }(\tau +1\mid \tau )={\mathsf{\Lambda }}_{\alpha }(\tau +1\mid \tau ){\tilde{C}}_{\alpha }^{-1}(\tau +1){\mathsf{\Lambda }}_{\alpha }{}^T(\tau +1\mid \tau )$$

  (40)

  $$\widetilde{\mathsf{{\rm K}}}(\tau +1)=\tilde{\lambda }(\tau +1\mid \tau ){\chi }^T{\left(\chi \tilde{\lambda }(\tau +1\mid \tau ){\chi }^T+\tilde{H}(\tau +1)\right)}^{-1}$$

  (41)

  $$\widehat{A}{(\tau +1\mid \tau +1)}_t=\widehat{A}(\tau +1\mid \tau )+\widetilde{\mathsf{{\rm K}}}(\tau +1)(\mathsf{\Gamma }(\tau +1)-H\widehat{A}(\tau +1\mid \tau ))$$

  (42)

  The estimates of the current iteration step are compared with those of the previous iteration and, if satisfied,

  $$\frac{\Vert \widehat{A}{(\tau +1\mid \tau +1)}_t-\widehat{A}{(\tau +1\mid \tau +1)}_{t-1}\Vert }{\Vert \widehat{A}{(\tau +1\mid \tau +1)}_{t-1}\Vert }\le \epsilon$$

  (43)

  then $\widehat{A}(\tau +1\mid \tau +1)=\widehat{A}{(\tau +1\mid \tau +1)}_t,\lambda (\tau +1\mid \tau +1)=\tilde{\lambda }(\tau +1\mid \tau )$, and the value of the pseudovariable can be updated, or the iteration can be repeated;
• $\tau =\tau +1$, and steps (3–5) are repeated until the end of filtering.

6. Simulated Cases

This section verifies the validity of the proposed method by providing two cases: one in which the state equation is a nonlinear equation and the measurement equation is a linear equation, and one in which the state and measurement equations are both nonlinear.

6.1. Case 1

Consider a nonlinear system in which the state equation is a nonlinear model and the measurement equation is a linear model:

$$\begin{array}{l}\{\begin{array}{l}{x}_1(k+1)=\left(0.8-0.5e^{-x_1^2(k)\left(1+e^{-0.015k}\right)}\right)x_1(k)-\left(0.3+0.9e^{-x_1{}^2(k)\left(1+0.5\sin \left(\frac{\pi }{2}k\right)\right)}\right)x_2(k)+w_1(k)\hfill \\ x_2(k+1)=1.2\left(1-e^{-0.8k}\right)x_2(k)+0.11x_1(k)+\cos \left(1+x_2^2(k)\right)+e^{-0.8k}{x}_1^4(k)+w_2(k)\hfill \end{array}\hfill \\ \{\begin{array}{c}{y}_1(k+1)=x_1(k+1)+v_1(k+1)\\ y_2(k+1)=x_2(k+1)+v_2(k+1)\end{array}\hfill \end{array}$$

where the initial value $x\left(0\right)$ is a random value of $\left[0,1\right]$, the initial estimation error covariance $P\left(0|0\right)=0.1\times diag(1,1)$, and the process noise and measurement noise have the following characteristics:

$$\begin{array}{c}{w}_1(k)~0.9N(0,0.01)+0.1N(0,0.2),w_2(k)~0.9N(0,0.02)+0.1N(0,0.2)\\ v_1(k)~0.9N(0,0.01)+0.1N(0,2),v_2(k)~0.9N(0,0.02)+0.1N(0,2)\end{array}$$

Figure 2 shows a diagram of the MTN identification system, while Figure 3 shows the estimated values of state variables $x_1$ and $x_2$ under the three filtering methods. From [21], we know the influence of ε is not significant compared with the kernel bandwidth σ. The parameters are set at $\epsilon ={10}^{-6}$.

Table 1. The mean squared error using the three methods in case 1.

		$\mathbf{MSE} \mathbf{of} {\mathit{x}}_1$			$\mathbf{MSE} \mathbf{of} {\mathit{x}}_2$
$\mathit{\sigma }$	$\mathit{\epsilon }$	MCEKF	H-MCKF	H-MCKF_R	MCEKF	H-MCKF	H-MCKF_R
$\sigma =2$	$\epsilon ={10}^{-6}$	0.2073	0.1815	0.1694	0.1000	0.0774	0.0738
$\sigma =5$	$\epsilon ={10}^{-6}$	0.1974	0.1244	0.1225	0.1100	0.0964	0.0921
$\sigma =10$	$\epsilon ={10}^{-6}$	0.2282	0.1669	0.1636	0.1158	0.0925	0.0888
$\sigma =20$	$\epsilon ={10}^{-6}$	0.2244	0.1602	0.1572	0.1160	0.0916	0.0880

and

Table 2. The mean relative error using the three methods in case 1.

		$\mathbf{MRE} \mathbf{of} {\mathit{x}}_1$			$\mathbf{MRE} \mathbf{of} {\mathit{x}}_2$
$\mathit{\sigma }$	$\mathit{\epsilon }$	MCEKF	H-MCKF	H-MCKF_R	MCEKF	H-MCKF	H-MCKF_R
$\sigma =2$	$\epsilon ={10}^{-6}$	0.3372	0.2405	0.2403	0.2354	0.2202	0.2084
$\sigma =5$	$\epsilon ={10}^{-6}$	0.3462	0.2953	0.2906	0.2679	0.2485	0.2448
$\sigma =10$	$\epsilon ={10}^{-6}$	0.3658	0.3052	0.2986	0.2745	0.2469	0.2426
$\sigma =20$	$\epsilon ={10}^{-6}$	0.3634	0.3009	0.2945	0.2753	0.2466	0.2419

show the mean squared error and the mean relative error, respectively, of the estimated values under the three algorithms, which are computed as averages over 100 independent Monte Carlo runs, with each run containing 50 time steps. When $\sigma =5$, the three algorithms all obtain better filtering results. Figure 4 and Figure 5 show the probability densities of the estimation errors when estimating the states $x_1$ and $x_2$, respectively, when the parameters are $\epsilon ={10}^{-6}$ and $\sigma =5$. All of the results confirm that the proposed H-MCKF (design method for a higher order extended Kalman filter based on maximum correlation entropy and a Taylor network system) can outperform the MCEKF (maximum correntropy extended Kalman filter) significantly when the system is disturbed by non-Gaussian processes and measurement noise, and the H-MCKF_R (H-MCKF with the remainder of the state equation) further improves the filtering performance of the H-MCKF.

6.2. Case 2

Consider a nonlinear system in which the state equation and the measurement equation are both nonlinear models:

$$\begin{array}{l}\{\begin{array}{l}{x}_1(k+1)=\cos \left(0.5x_1\left(k\right)+\frac{2.5x_2\left(k\right)}{1+x_1{}^2\left(k\right)+8\cos (1.2k)}\right){+\mathrm{w}}_1(k)\hfill \\ x_2(k+1)=\sin (x_1{}^2{\left(k\right))+\mathrm{w}}_2(k)\hfill \end{array}\hfill \\ \{\begin{array}{c}{y}_1\left(k\right)=\cos (x_1\left(k\right)+\sin (x_1{}^3\left(k\right)))+v_1(k+1)\\ y_2\left(k\right)=\sin (x_2\left(k\right)-\sin (x_2{}^3\left(k\right)))+v_2(k+1)\end{array}\hfill \end{array}$$

where the initial value $x\left(0\right)$ is a random value of [0, 1], the initial estimation error covariance $P\left(0|0\right)=0.1\times diag(1,1)$, and the process noise and measurement noise have the following characteristics:

$$\begin{array}{c}{w}_1(k)~0.9N(0,0.01)+0.1N(0,0.2),w_2(k)~0.9N(0,0.02)+0.1N(0,0.2)\\ v_1(k)~0.9N(0,0.01)+0.1N(0,2),v_2(k)~0.9N(0,0.02)+0.1N(0,2)\end{array}$$

Figure 6 shows a diagram of the MTN identification system, while Figure 7 shows the estimated values of state variables $x_1$ and $x_2$ under the three filtering methods. Similar to case 1, the parameters are set at $\epsilon ={10}^{-6}$.

Table 3. The mean squared error using the three methods in case 2.

		$\mathbf{MSE} \mathbf{of} {\mathit{x}}_1$			$\mathbf{MSE} \mathbf{of} {\mathit{x}}_2$
$\mathit{\sigma }$	$\mathit{\epsilon }$	MCEKF	H-MCKF	H-MCKF_R	MCEKF	H-MCKF	H-MCKF_R
$\sigma =2$	$\epsilon ={10}^{-6}$	0.4017	0.1230	0.1219	0.2090	0.0907	0.0883
$\sigma =5$	$\epsilon ={10}^{-6}$	0.1148	0.1241	0.1233	0.2542	0.1200	0.1183
$\sigma =10$	$\epsilon ={10}^{-6}$	0.3221	0.1254	0.1248	0.2220	0.1207	0.1193
$\sigma =20$	$\epsilon ={10}^{-6}$	0.4040	0.1257	0.1251	0.2218	0.1208	0.1196

and

Table 4. The mean relative error using the three methods in case 2.

		$\mathbf{MRE} \mathbf{of} {\mathit{x}}_1$			$\mathbf{MRE} \mathbf{of} {\mathit{x}}_2$
$\mathit{\sigma }$	$\mathit{\epsilon }$	MCEKF	H-MCKF	H-MCKF_R	MCEKF	H-MCKF	H-MCKF_R
$\sigma =2$	$\epsilon ={10}^{-6}$	0.5106	0.2337	0.2306	0.3742	0.2355	0.2316
$\sigma =5$	$\epsilon ={10}^{-6}$	0.2147	0.2551	0.2530	0.3070	0.2652	0.2645
$\sigma =10$	$\epsilon ={10}^{-6}$	0.4527	0.2570	0.2553	0.3824	0.2661	0.2659
$\sigma =20$	$\epsilon ={10}^{-6}$	0.4764	0.2575	0.2558	0.3789	0.2663	0.2662

show the mean squared error and the mean relative error, respectively, of the estimated values under the three algorithms, which are computed as averages over 100 independent Monte Carlo runs, with each run containing 50 time steps. When $\sigma =5$, the three algorithms all obtain better filtering results. Figure 8 and Figure 9 show the probability densities of the estimation errors when estimating the states $x_1$ and $x_2$, respectively, when the parameters are $\epsilon ={10}^{-6}$ and $\sigma =5$. All of the results confirm that the proposed H-MCKF can outperform the MCKF significantly when the system is disturbed by non-Gaussian processes and measurement noise, and the H-MCKF_R further improves the filtering performance of the H-MCKF when the state and measurement equations are both nonlinear.

7. Conclusions

This paper considered a wide range of filter design problems for the state estimation of multivariable dynamic systems, which consist of a strong nonlinear dynamic model and a strong nonlinear observation model. Firstly, we transformed those strong nonlinear models into a higher order polynomial series using a multidimensional Taylor network. Secondly, all higher order items in the polynomial series were defined as hidden variables. Those higher order series were then rewritten as their pseudolinear equivalents. Thirdly, dynamic relationships between all hidden variables and known variables were constructed using the multidimensional Taylor network. Combining the original model of pseudolinearization with the higher order hidden variable dynamic model, linear dynamic models fitted to a standard Kalman filter were presented. Finally, considering that a finite number of samples from modeling error can be obtained, we built the higher order extended Kalman filter based on maximum correlation entropy, and acquired better filter performance than offered by the existing MCEKF [22].

Outlook: There exist several challenges worthy of further research. Firstly, the proposed higher order extended Kalman filter based on maximum correlation entropy is an online iteration process that obtains state estimation constantly, but, as such, it loses one important function possessed by the standard Kalman filter: the ability to operate in real time. Secondly, the linearized model parameters of the original nonlinear model and the hidden variable dynamic model were identified by local time period data; thus, they need to be updated with new time period data in order to fit the time dynamics of the system. Thirdly, in this paper, on the basis of defining all of the hidden variables, we established a linear form of the strong nonlinear model in an expanded state with the original variables and all hidden variables, and obtained better estimation performance than that of a standard EKF; if measurements can be expanded in the same manner as state, we believe that such a filter may offer better estimation performance than the one established by this paper.