# Differential Evolution Optimized a Second-Order Divided Difference Particle Filter

Authors: Ting Cao; Huo-tao Gao; Chun-feng Sun; Guo-bao Ru
Journal: Journal of Electrical and Computer Engineering
 (2020)

---

## Abstract

In order to improve the estimation accuracy of particle filter algorithm in a nonlinear system state estimation problem, a new algorithm based on the second-order divided difference filter to generate the proposed distribution and the differential evolution algorithm for resampling is proposed. The second-order divided difference based on Strling’s interpolation formula is used to generate approximations to nonlinear dynamics, which avoids the evaluation of the Jacobian derivative matrix and is easy to implement. Cholesky factorization is used to ensure the positive definiteness of the covariance matrix. The truncated errors of the local linearization are reduced to a certain extent, and the approximation degree of the proposed distribution to the posterior probability of the system state is improved. The differential evolution algorithm is used to replace the traditional resampling algorithm, which effectively mitigates the problem of particle degradation. Monte Carlo simulation experiments show the effectiveness of the new algorithm.


---

## Body

## 1. Introduction

Particle filters (PF) are an effective solution to a nonlinear dynamic and/or measurement filtering problem, which has been widely used in many fields such as economics [1, 2], visual tracking [3, 4], navigation position [5], and radio communication [6]. The PF is a recursive Bayesian estimation method based on the Monte Carlo simulation. It exploits a set of random particles (samples) and their associated weights to represent the posterior probability density function. Despite achieving the state-of-the-art performance, existing PF approaches are restricted to two aspects: the sample degeneracy and impoverishment. The sample degeneracy appears when the weight focuses on a few particles only and the rest have negligible weight, and impoverishment is due to the loss of diversity among the particles in resampling process when high-weight particles are chosen many times. In general, there are three methods to improve the performance of the PF. The first is to increase the particle number. This is often impractical because increasing the number of particles yields very limited improved estimation accuracy, but requires significantly increased computational burden. Therefore, we rely on two other methods: (1) select an effective proposal distribution and (2) resampling.To obtain better proposal distribution, many nonlinear filtering techniques have been put forward. For instance, by linearizing nonlinear functions around the current state at first-order Taylor approximations, extended Kalman filter (EKF) is used to generate the proposal distribution, and then the extended particle filter (EPF) is proposed. However, the EKF has slow convergence, low estimated accuracy and may suffer from large errors when the system has strong nonlinearity. The unscented particle filter (UPF) [7] utilizes the unscented Kalman filter (UKF) to generate the proposal distribution. By translating the analytic integral operator into the approximate summation operator through a set of deterministic points, the UKF can accurately calculate the mean and covariance of nonlinear functions up to the second order (third in Gaussian prior) of the Taylor series expansion. Therefore, the UPF outperforms EPF in accuracy. However, as the system dimension increases, the computation burden increases rapidly, this limits the real-time performance of UPF. To address the computation load, several simplex point selection strategies have been developed. In 2003, Julier proposed spherical simplex unscented transformation (SSUT) [8]. In [9], the spherical simplex unscented particle filter resulted from using a spherical simplex unscented Kalman filter (SSUKF) for the proposal distribution generation within a particle filter framework. Zhao [10] proposed a global sampling strategy which generated proposal distribution from the mean and covariance value of the whole particles instead of from each particle. Aiming at the degeneracy problem of particle filter algorithm, many researchers proposed some improved algorithms, such as the gradient particle filtering [11] and the box particle filtering [12].The second-order divided difference (DD2) filter [13] is based on a similar linear regression to the UKF, but is able to more accurately estimate the state covariance resulting in a more precise set of weighted sample points to estimate from. Besides, the DD2 filter provides a faster processing speed than the UKF because it does not need to predict forward every positive and negative sigma point in separate stages from when the measurement prediction, measurement prediction covariance, and cross-covariance are interpreted.Resampling is a key operation in particle filtering. Some improved techniques based on intelligent optimization approaches were therefore proposed to modify the resampling scheme, which includes Genetic Particle Filter [14], Ant Colony Particle Filter [15] and Particle Swarm Particle Filter [16]. Recently, scholars are eager to the research of Artificial Neural Network (ANN) and their applications [17–22]. The application of BP or Convolution Network/PF has indicated an attractive research direction in combining particle filtering and ANN optimization.In this paper, we propose using the DD2 filter to generate the proposal distribution. Meanwhile, aiming at the problem of Particle Filter algorithm such as particles degeneracy and particles impoverishment, a differential evolution algorithm is used to make the particles of DD2 filter move toward the global optimum, which optimizes the resampling process and relieves the problem of particles degeneracy and impoverishment.The rest of the paper is organized as follows. In Section2, we briefly reviewed the general second-order divided difference filter. Section 3 gives details of the proposed particle filter. Experiments and analysis are presented in Section 4. Finally, the work is concluded in Section 5.


## 2. The Second-Order Divided Difference Filter

In this section, a generalized version of the DD2 filter will be outlined. DD2 filter exploits polynomial approximations obtained with a Stirling interpolation formula for the derivation of state estimators for nonlinear systems. Conceptually, the principle underlying the DD2 filters is much like that of the EKF and its higher order relatives. The estimator becomes more precise than estimators based on Taylor approximations; yet the implementation is significantly simpler as no derivatives are required. Additionally, by using the Cholesky factors of covariance matrices, the numerical properties of estimation algorithms are generally improved.Many nonlinear filtering problems can be described in the form of the dynamic state space (DSS) model as follows:(1)xk=fxk−1,wk−1Process  equation,(2)yk=hxk,vkMeasurement  equation,where xk and yk are the state vector and observation at time instant k. Process system noise wk−1 and the measurement noise vk are random vectors that derive the dynamic system through the nonlinear state transition function f⋅ and the nonlinear measurement function h⋅, respectively.We consider expanding the nonlinear functiony=fx, by polynomial approximations based on String’s interpolation formula. The second-order expansion around the point x=x¯ gives the following approximation:(3)fx≈fx¯+fDD′x¯x−x¯+fDD″x¯2!x−x¯2,where(4)fDD′x¯=fx¯+λ−fx¯−λ2λ,fDD″x¯=fx¯+λ+fx¯−λ−2fx¯λ2,where λ is the interval step-size for approximation. Inserting the full Taylor series in place of fx¯+λ and x¯−λ, we can obtain(5)fx¯+fDD′x¯x−x¯+fDD″x¯2!x−x¯2=fx¯+f′x¯x−x¯+f″x¯2!x−x¯2+f3x¯3!λ2+f5x¯5!λ4+⋯x−x¯+f4x¯4!λ2+f6x¯6!λ4+⋯x−x¯2.The first three terms on the right hand side of equation (5) are identical with the second order Taylor series expansion and are independent of the interval length λ. The “remainder” term corresponds to the higher order terms of the Taylor series expansion of f and is controlled by λ. Proper interval lengths can make sure that the remainder term in some sense will be close to the higher order terms of the full Taylor series. As we can see, it is the remainder term that makes the interpolation formula more accurate than Taylor approximation in some applications.The following procedure outlines the implementation of the DD2 filters.(1)
Square Cholesky factorization:(6)P^k−1=S^xS^xT,P¯k−1=S¯xS¯xT,(7)Q=SwSwT,R=SvSvT,whereS^x, S¯x, Sw, and Sv are the Cholesky factorization of the predicted state error covariance P^k−1, corrected state error P¯k−1, process noise covariance Q and measurement noise covariance R, respectively.(8)Sxx^k−11=12λfix^k−1+λS^x,j,w¯k−1−fix^k−1−λS^x,j,w¯k−1,(9)Sxwk−11=12λfix^k−1,w¯k−1+λSw,j−fix^k−1,w¯k−1−λSw,j,(10)Sxx^k−12=λ2−12λ2fix^k−1+λS^x,j,w¯k−1+fix^k−1−λS^x,j,w¯k−1−2fix^k−1,w¯k−1,(11)Sxwk−12=λ2−12λ2fix^k−1,w¯k−1+λSw,j+fix^k−1,w¯k−1−λSw,j−2fix^k−1,w¯k−1,whereS^x,j denotes the jth column of S^x and similarly for the other factors. The Cholesky factor of the predicted covariance S¯xk is obtained by the House-holder transformation of the compound matrix.(12)S¯xk=Sxx^k−11Sxwk−11Sxx^k−12Sxwk−12.(2)
The predicted statex¯k:(13)x¯k=λ2−nx−nwλ2fx^k−1,w¯k−1+12λ2∑p=1nxfx^k−1+λS^x,p,w¯k+fx^k−1−λS^x,p,w¯k−1+12λ2∑p=1nwfx^k−1,w¯k−1+λSw,p+fx^k−1,w¯k−1−λSw,p,wherenx is the size of the state dimension and nw is the size of the process noise dimension.(14)Syx¯k1=12λhix¯k+λS¯x,j,v¯k−hix¯k−λS¯x,j,v¯k,(15)Syvk1=12λhix¯k,v¯k+λSv,j−hix¯k,v¯k−λSv,j,(16)Syx¯k2=λ2−12λ2hix¯k+λS¯x,j,v¯k+hix¯k−λS¯x,j,v¯k−2hix¯k,v¯k,(17)Syvk2=λ2−12λ2hix¯k,v¯k+λSv,j+hix¯k,v¯k−λSv,j−2hix¯k,v¯k,(18)S^yk=Syx¯k1Syvk1Syx¯k2Syvk2.The predicted measurementy¯k:(19)y¯k=λ2−nx−nvλ2hx¯k,v¯k+12λ2∑p=1nxhx¯k+λS¯x,p,v¯k+hx¯k−λS¯x,p,v¯k+12λ2∑p=1nvhx¯k,v¯k+λSv,p+hx¯k,v¯k−λSv,p,wherenv is the size of the observation noise dimension.(3)
Update stage:(20)x^k=x¯k+Kkyk−y¯k,(21)S^xk=S¯xk−KkSyx¯k1KkSyvk1KkSyx¯k2KkSyvk2,(22)Kk=S^xkSyx¯k1TS^ykS^ykT−1,(23)P^k=S^xkS^xkT.


## 3. Differential Evolution Optimized a DD2 Particle Filter

Differential evolution (DE), proposed by Storn and Price [23], is an efficient and powerful population-based stochastic search technique for solving optimization problems over continuous search domain, which has been broadly applied in many scientific and engineering fields. It generates new candidate solutions by combining the parent individual and several other individuals of the same population, capable of dealing with nondifferentiable, nonlinear, and multimodal objective functions.In this paper, we adopt the DD2 filter to generate the proposal distribution. Meanwhile, the latest measurements are employed to promote the proposal distribution for generating particles more intelligently. After the DD2 importance sampling, the DE technique is utilized to optimize the resample process. The procedure of optimization views each particle as a target vector of the current population in the search space. There is a weight associated with each particle to determine if the particle is good or not. The population of the DE evolves by mutation, crossover, and selection. All of these operations can improve the valid particle number and reduce sample impoverishment. Thus, the new particles can be better approximate the true state by exploiting new areas, and the effectiveness and diversity of particles can be improved.The description of the proposed algorithm is summarized as follows:(1)
Initiation:k=0DrawN particles x0ii=1N from the prior px^0 and set weight ω0i=1/N, i=1,…,N.(2)
Importance sampling:(i)
Update particlesxkii=1N with the DD2 filter using equations (6)∼(23) to obtain x¯kii=1N and P¯kii=1N.(ii)
Sample particlesx^kii=1N from the proposal distribution qx^kix^0:k−1i,z1:k=Nx¯ki,P¯ki(iii)
Compute the particle weight and normalize:(24)ωki∝ωkipzkx^kipx^kix^k−1i−1qx^kix^0:k−1i,z1:k,(25)ω˜ki=ωki∑i=1Nωki.(3)
ResamplingIn aD-dimensional real parameter space ℛD, let the sampled particles x^kii=1N denote the state vector of current population x^i,g=x^i,1,g,…,x^i,D,g,i=1,…,N at generation G=0. Regard each particle and its weight as the target vector and fitness function, respectively. The resampling is composed of the following three steps.(i)
MutationHere, we use the “DE/rand/1” mutation strategy to generate mutant vectorVl,g.(26)Vl,g=x^r1,g+F⋅x^r2,g−x^r3,g,l≠r1≠r2≠r3,The indicesr1, r2, and r3 are mutually distinct integers randomly chosen from the range [1, N], and all are different from the base index l. F is a positive mutation scaling factor which affects the differential variation between individuals.
(ii)
Crossover
To improve the potential diversity of the population, a crossover operation comes into play after generating the donor vector through mutation.(27)Ul,g=Uj,l,g=Vj,l,g,ifrandj0,1≤Crorj=jrand,x^j,l,g,otherwise.
Wherejrand is a randomly chosen integer between 0 and 1, randj0,1 is a uniformly distributed random number in [1, D], which is generated for each j, and Crε0,1 is the crossover rate. By mutation and crossover operation, we can enlarge the search space and the particles can be shifted around real state.
(iii)
Selection
With the abovementioned crossover and mutation process, the system generates a large number of particles, and thus the new particles can better approximate the true state by exploiting new areas. Subsequently, we select particles which have large fitness (weight) value. The higher the fitness value is, the more likely the particle is selected in the DE algorithm.(28)x˙l,g+1=Ul,g,if weightUl,g≥weightx^l,g.x^l,g,otherwise.When all individuals of the new generation are calculated, the mutation step is carried out to maintain the diversity of the population. That is to say, the mutation, crossover and selection procedure repeat until the optimum is found or the maximum generation numberG is reached. For this reason, the DE adopts an iterative approach for particle generation. As the particles generally become better and better as the iterations go on, the DE can decisively improve the particle quality. Thanks to the improved particle quality, particle samples are moved towards regions where particles have larger values of posterior density function. As a result, the impoverishment problem suffered in the particle filter can be effectively handled, and the degeneracy phenomenon is alleviated obviously.(4)
Output: a set of samples that can be used to approximate the posterior distribution as follows:(29)x˜k=∑i=1Nx˙kiω˜˙ki.The basic steps of the proposed algorithm are summarized in Figure1.Figure 1
The basic steps of the proposed algorithm.


## 4. Simulation Experiments

In this section, the estimation performance and computational cost of the proposed particle filters are evaluated on a state estimation problem and compared with that of the standard sampling importance-resampling particle filter (SIR-PF), extended particle filter (EPF), unscented particle filter (UPF), and the second-order divided difference particle filter (DD2PF). A typical nonlinear process model and measurement model are as follows [24–29]:(30)xk=1+sin0.04πk+0.5xk−1+uk−1,(31)yk=0.2xk2+vk,k≤30,0.5xk+vk−2,k>30,where uk is a Gamma ςa3,2 random variable modelling the process noise. The measurement noise vk is drawn from a Gaussian distribution Nvk;0,10−5. The scale factor F, the crossover probability Cr, and the maximum number of generations G are set based on trial-and-error procedures. In this case, we set F=0.9, Cr=0.6, and G=20.To qualify a fair comparison of the estimates produced by each of the five filters, the estimates are averaged across a Monte Carlo simulation consisting of 200 runs. Each run is implemented with a different noise sample. Ten and two hundred particles are used, respectively. In order to evaluate performance of the proposed algorithm, we use the traditional measure of performance: the Root Mean Squared Error (RMSE). The RMSE of each run is defined by(32)RMSE=1M∑k=1Mx˜k−xk2.For reference, the true state and the state estimate of the different particle filter for a random Monte Carlo run are presented in Figures2 and 3. As Figures 2 and 3 show, states estimated by DD2PFDE is the most closed to the true state. Figures 4 and 5 show the variation of estimation RMSE with observation time for all particle filters. In addition, Table 1 summarizes the performance of different filters for 200 MC tests. The table demonstrates the means and variances of the mean square-error of the state estimates as well as the average processing time in seconds of each filter.Figure 2
Filtering result comparison of five filters in a random simulation run, where particle numbersN=10.Figure 3
Filtering result comparison of five filters in a random simulation run, where particle numbersN=200.Figure 4
RMSE of SIR-PF, EPF, UPF, DD2PF, and DD2PFDE for 200 Monte Carlo simulations, where particle numbersN=10.Figure 5
RMSE of SIR-PF, EPF, UPF, DD2PF, and DD2PFDE for 200 Monte Carlo simulations, where particle numbersN=200.Table 1
Estimation results and computation times averaged over 200 Monte Carlo runs.
Algorithm
Mean
Variance
Time(s)
Particles
SIR-PF
0.9431
0.0218
0.13011
10
EPF
0.4915
0.0076
0.21287
UPF
0.4911
0.0051
0.39452
DD2PF
0.4101
0.0050
0.30124
DD2PFDE
0.1598
0.0043
0.51271
SIR-PF
0.5803
0.0372
1.45103
200
EPF
0.3503
0.0107
3.55714
UPF
0.2220
0.0082
6.85312
DD2PF
0.1527
0.0041
5.56231
DD2PFDE
0.0769
0.0031
9.43107
SIR-PF
0.2190
0.0068
2.67380
1800As can be seen from Figure2, the estimates for the DD2PF, UPF, EPF, and the SIR-PF may still be far away from the true values. The reason for the failure is that only one particle may be significantly weighted and the others are much less weighted. But on the contrary, in Figure 2, the proposed particle filter DD2PFDE can track successfully throughout the whole period. Owing to the adding of mutation, crossover and selection, which can be seen in Section 3, and particles in DD2PFDE can exploit more wide areas to find suitable particles to approximate the target’s true state. And DD2PFDE can maintain the diversity of the particle set. All of these aforementioned guarantee that the proposed filter can successfully estimate the target.In Figures3 and 5, improvements on the five algorithms with 200 particles are readily seen, not only with estimation accuracy, but also in the robustness of the performance. Only with greater numbers, will more particles be considerably weighted through the likelihood function, and a posterior estimate is calculated from the normalized weighted sum of particle estimates, resulting to be much closer to the true values.As tabulated in Table1, when the number of particles used by five algorithms is the same, we can see the SIR-PF, EPF, DD2PF, UPF, and DD2PFDE runtime increases in turn. EPF and UPF take longer than SIR-PF due to the addition of EKF and UKF filtering before particle sampling. The EKF provides a faster processing speed than the UKF because it does not need to augment the state dimension, select, and calculate the sigma points in separate procedures. Comparing with the other four algorithms, the DD2PFDE is the most time consuming. This is due to the fact that both the DD2 filter and the DE iteration are introduced into DD2PFDE.In comparison with UKF, DD2 filter presents better accuracy and less computation cost and complexity. Therefore, DD2PF provides better performance than UPF. In order to evaluate the efficiency of the proposed algorithm correctly, we increase the number of PF particles to 1800. Two hundred Monte Carlo experiments show that the RMSE mean is 0.2190, and the average computation time is 2.6738 s. Note that 10 DD2PFDE particles achieve higher accuracy (the RMSE mean is 0.1598) and spend less time (the average run time is 0.51271 s) than 1800 PF particles, which indicate that the particle utilization efficiency of the proposed algorithm is higher. Because DE can mutate and select sampling with fitness function, better particles are selected and various particles are generated. By using the DE schemes, the diversity of the particles is increased. As a result, the performance of estimation is enhanced. Despite the proposed algorithm running time is longer than other algorithms; it can decrease the degeneracy phenomenon and alleviate the problem of sample impoverishment.


## 5. Conclusions

In this paper, we propose an efficient second-order divided difference particle filter based on differential evolution. The filter adopts the DD2 filter to generate the proposal distribution, and the differential evolution algorithm to optimize the resampling process. The DD2 filter linearizes the nonlinear dynamic model by using a multivariable extension of Stirling’s interpolation formula rather than the derivative-based Taylor series approximation. The filter provides excellent accuracy without the need to analytically calculate Jacobians. In addition, by using square root factorizations of covariance, the numerical properties of estimation are generally improved. Because the proposed algorithm incorporates the latest measurement information into the sampled particles, the sampled particles are optimized from the source, and the differential evolution algorithm is used to iteratively optimize the sampled particles, thus effectively alleviate the problem of particle depletion and making improved algorithm particle utilization. As a result, the accuracy of the filter would be further enhanced. We illustrated the feasibility of our proposed filtering algorithm for testing nonlinear filtering problems. The simulation results demonstrate that the proposed particle filter can perform better than the other sequential particle filters and can be easier to implement.



---

*Source: 9541624-2020-02-24.xml*