# Improved Particle Filter Using Clustering Similarity of the State Trajectory with Application to Nonlinear Estimation: Theory, Modeling, and Applications

Authors: Ziquan Jiao; Zhiqiang Feng; Na Lv; Wenjing Liu; Haijian Qin
Journal: Journal of Sensors
 (2021)

---

## Abstract

A clustering similarity particle filter based on state trajectory consistency is presented for the mathematical modeling, performance estimation, and smart sensing of nonlinear systems. Starting from an information fusion model based on the consistency principle of the spatial state trajectory, the predicted observation information of the current particle filter (original trajectory) and future multistage Gaussian particle filter (modified trajectory) are selected as the state trajectories of the sampling particles. Clustering similarity methods are used to measure the state trajectories of the sampling particles and the actual system (reference trajectory). The importance weight of a first-order Markov model is updated with the measurement results. By integrating the targeted compensation scheme of the latest measurement information into the sequential importance sampling process, the adverse effects of the particle degradation phenomenon are effectively reduced. The convergence theorems of the improved particle filter are proposed and proved. The improved filter is applied to practical cases of nonlinear process estimation, economic statistical prediction, and battery health assessment, and the simulation results show that the improved particle filter is superior to traditional filters in estimation accuracy, efficiency, and robustness.


---

## Body

## 1. Introduction

Nonlinear phenomena are common in natural engineering technology. As a popular research topic with important theoretical and practical value in solving nonlinear problems, state estimation has been applied to problems such as target tracking and navigation, fault diagnosis and detection, process feedback and control, biochemical reaction and extraction, and economic prediction and control. Nonlinear state estimation applies to a wide range of fields, especially the industrial field. Applications include longitudinal vehicle speed estimation [1], fault detection of piezoceramic actuators [2], battery health assessment [3], state detection of impending rollover [4], and state estimation of dynamic systems with hysteresis [5]. Many solutions have been proposed, such as the Luenberger observer [6], robust observer [7], Gaussian process regression, Kalman filter [8], proportional integral observer [9], unknown input observer [10], high-gain observer [11], and nonsmooth observer [12]. For example, on the premise of satisfying the Gaussian noise distribution, the Kalman filter (KF) calculates the conditional probability density of random variables using a recursive formula and iteratively updates the linear minimum variance estimate. From this come extended KF, unscented KF, invariant extended KF [13], and adaptive extended KF [14]. The above methods have many advantages and wide applications for nonlinear system state estimation, but there is still much room for optimization and improvement in terms of nonlinear complexity and environmental noise uncertainty of different practical applications.With the rapid development of computer technology, the particle filter (PF) algorithm based on Bayesian and Monte Carlo theories has shown many advantages and considerable potential to solve estimation problems involving nonlinear and non-Gaussian systems. The prototype algorithm, sequential importance sampling (SIS), was formed in the middle of the last century and is mainly applied to physics and automatic control applications. Due to inherent sample degradation and computer hardware limitations, the study of the PF algorithm slowed until 1993, when Gordon et al. [15, 16] introduced a resampling strategy to SIS and developed sequential importance resampling (SIR), which improved the method and laid the theoretical foundation for the PF. With the development of stochastic probability theory and Monte Carlo methods, the auxiliary particle filter (APF) [17] and Gaussian (sum) particle filter (GPF) [18, 19] were proposed. Introduced by Guarniero et al. [20], APF is based on the idea that the latest observation will approach the optimal proposed distribution if an auxiliary variable is imported to represent the prior probability of the current state. When the system noise is strong, the filtering accuracy is difficult to guarantee due to the lack of information. The GPF algorithm of Sun et al. [21] uses a Gaussian distribution to estimate the posterior probability density function (PDF) of a system state under the basic SIS framework, and the mean and variance are recursively obtained. The filtering effect depends heavily on the problem’s degree of nonlinearity and is limited to the dimension of the system variables [22]. With the high complexity of current natural engineering structures, the degree of nonlinearity of systems is growing. Although these algorithms somewhat improve the performance of PFs, issues remain, such as low accuracy of filter estimation and poor stability due to particle degradation and depletion, which do not meet the needs of modern engineering. Moreover, as a probabilistic method, the nonlinear estimation of PF leads to uncertainty in the result [23].For this reason, a clustering similarity particle filter (CSPF) based on the consistency principle of the spatial state trajectory is presented. The clustering similarity method is used to measure the state distance between the actual system and the sampling particles, including the observation information of the current state filtering and future multistage state prediction, to guide the generation and improvement of new distributions and update the weight calculation of the importance sampling process. This makes up for using the prior PDF instead of the importance function in the standard PF algorithm, which can prevent the occurrence of particle degradation and significantly improve the accuracy and robustness of estimation. The resampling strategy in the traditional PF algorithm is abandoned to eliminate particle depletion, which improves the quantization accuracy of uncertainty and efficiency of the algorithm. The above methods adopt the idea of simply modifying the proposed distribution. The designed method uses clustering theory to measure the similarity of observation information corresponding to multistage (fromk to k+L+l) state trajectories so as to guide the generation and updating of the latest proposal distribution, which significantly improves the computational complexity of the designed method. To ensure the efficiency of the improved method for different nonlinear state estimation applications, the following two aspects can be improved. (a) The order of the state trajectory should be selected reasonably. Theoretically, the higher the order of trajectory selection, the more accurate the corresponding observation information can express the actual state, and the more accurate the estimation result, but the computational efficiency is greatly decreased. (b) The number of sampling particles N should be reduced appropriately. With the increase of the number of particles, the sampling probability density function will gradually approach the probability distribution of the actual state. While improving the accuracy of state estimation, the computational effort will increase. For these reasons, it is necessary to coordinate the contradiction between estimation accuracy and computational efficiency. The nonlinear state estimation results are largely affected by the signal information, which involves the quality and scale of the research object dataset, appropriate parameter identification, and state tracking training methods. In different research and application objects, increasing the quality and scale of experimental datasets containing more physical model information can improve the state estimation performance of data mining; an appropriate parameter training method can ensure that the model can obtain as much useful information as possible from the dataset and help to establish an appropriate state space model. Based on the above measures, compared to traditional estimation methods under the premise of consistent preconditions, the designed method can greatly improve the accuracy of state results, and it also improves computational efficiency.The remainder of this paper is structured as follows. Section2 discusses nonlinear system theory. In Section 3, an improved CSPF algorithm is proposed based on the analysis of the defects of the PF. Section 4 provides a theoretical explanation of the improved algorithm and proves the relevant theorems. Section 5 compares the simulation results of the proposed algorithm and the traditional improved PF algorithm. We discuss our conclusions in Section 6.


## 2. Theory Statement

We summarize the basic definitions and properties of the state space and optimal Bayesian recursion theory of nonlinear systems.


### 2.1. State Space Model

We setΩ,P,F as a random probability space and define two actual vector stochastic processes: X=Xk,k∈N and Y=Yk,k∈N\0, where sample space Ω is the set of all possible outcomes, event space F is a set of outcomes in the sample space, probability function P assigns a probability to each event in the event space, X is the state process, and Y is the observation information. Let nx and ny be the dimensions of the state X and observation information Y, respectively, corresponding to the state space, and define BRn as the Borelσ‐algebra set on n-dimensional Euclidean space Rn. Most nonlinear systems can take the form of a dynamic state space (DSS) [24]:
(1)Xk=FXk−1,ωkstatemodel,Yk=HXk,υkobservationequation,where k is the discrete-time (stage) index, Xk∈Rn is the set of system states at time k, and Yk∈Rn is the observation information corresponding to state Xk. F∗ and H∗ are known state transition and observation functions, respectively, corresponding to the state transition kernel PDF KXk∣Xk−1 and observation likelihood PDF gYk∣Xk in the statistical description. The system shift noise ωk∈Rn and measurement noise υk∈Rm are independent and identically distributed (i.i.d.) sequences that obey any PDF form.The state spaceX follows the first-order Markov process; i.e., the state Xk of the current moment is only related to the state Xk−1 of the previous moment. Assuming an initial distribution X0~μdX0, the probability density functions of the state transition kernel PDF KXk∣Xk−1 and observation likelihood PDF gYk∣Xk are Lebesgue measures:
(2)KXk∣Xk−1=PωXk−FXk−1,gYk∣Xk=PυYk−HXk,where Pω and Pυ are the probability functions under the influence of shift noise ωk and measurement noise υk, respectively.


### 2.2. Optimal Bayesian Recursion Theory

In Bayesian theory [15], the state Xk of a nonlinear system at time k is updated based on the observation information Y1:k to obtain a minimum mean squared error (MSE) estimate. The optimal state estimation of the system is the conditional expected value EXk∣Y1:k of the posterior PDF PXk∣Y1:k. Based on the premise that the state variable and observation function follow a first-order Markov process, the posterior PDF PXk∣Y1:k is obtained by two steps of recursive iterations: prediction and updating. Xk:l≜Xk,Xk+1,⋯,Xl and Yk:l≜Yk,Yk+1,⋯,Yl are defined as the space path information of the state process and observation likelihood from time k to l, respectively.


#### 2.2.1. Prediction

Combined with the transition kernel PDFKXk∣Xk−1 calculated by the state space equation, the prior PDF PXk∣Y1:k−1 is predicted using the Chapman–Kolmogorov equation and the posterior PDF PXk−1∣Y1:k−1 at time k−1:
(3)PXkY1:k−1=∫PXk,Xk−1∣Y1:k−1dXk−1=∫PXk∣Xk−1,Y1:k−1PXk−1∣Y1:k−1dXk−1=∫KXk∣Xk−1PXk−1∣Y1:k−1dXk−1,where the state PDF PXk,Xk−1∣Y1:k−1 specifies the conditional probability of Xk∩Xk−1 given Y1:k−1, and PXk∣Xk−1,Y1:k−1 specifies the conditional probability of Xk given Xk−1∩Y1:k−1.


#### 2.2.2. Updating

The prior PDFPXk∣Y1:k−1 is updated using the observation likelihood PDF gYk∣Xk at time k, and the posterior PDF PXk∣Y1:k of state Xk is obtained as
(4)PXkY1:k=PYk∣Xk,Y1:k−1PXk∣Y1:k−1PYk∣Y1:k−1=gYk∣XkPXk∣Y1:k−1∫gYk∣XkPXk∣Y1:k−1dXk,where the observation PDF PYk∣Xk,Y1:k−1 specifies the conditional probability of Yk given Xk∩Y1:k−1, and PYk∣Y1:k−1 specifies the conditional probability of Yk given Y1:k−1.The system state PDF measureπk:l∣mdXk:l is defined as
(5)πk:l∣mdXk:l≜PXk:l∣Y1:m.The joint posterior PDF measureπ0:k∣kdX0:k can be obtained by Bayes’ theorem (Equations (3) and (4)):
(6)π0:k∣k−1dX0:k=∫RnxKdXk∣Xk−1π0:k−1∣k−1dX0:k−1prediction,(7)π0:k∣kdX0:k=gYk∣Xkπ0:k∣k−1dX0:k∫RnxgYk∣Xkπ0:k∣k−1dX0:kupdating.The marginal posterior PDF measureπk∣k−1dXk is obtained by
(8)πk∣k−1dXk=∫RnxKdXk∣Xk−1πk−1∣k−1dXk−1prediction,(9)πk∣kdXk=gYk∣Xkπk∣k−1dXk∫RnxgYk∣Xkπk∣k−1dXkupdating.Similarly,xk:l≜xk,xk+1,⋯,xl and yk:l≜yk,yk+1,⋯,yl are, respectively, defined as the sampling particle path information of the state process and observation likelihood from time k to l.Definition 1.
Suppose thatν is a probability measure and φ represents an arbitrary function, Θ and Ψ are arbitrary function variables, K is the PDF of the transfer kernel satisfying the Markov process, and the following calculation method is defined:
(10)ν,φ≜∫φν,νKΘ≜∫νdxKΘ∣x,Kφx≜∫KdΨ∣xφΨ.
According to the above symbols, for any functionφ, Bayesian theory (prediction and updating processes) can be redefined, using Equation (8), as
(11)πk∣k−1,φ=∫Rnxpxk∣y1:k−1φxkdxk=∫Rnx∫RnxKxk∣xk−1pxk−1∣y1:k−1dxk−1φxkdxk=∫Rnxpxk−1∣y1:k−1∫RnxKxk∣xk−1φxkdxkdxk−1=πk−1∣k−1,Kφprediction.
From Equation (9), it is concluded that
(12)πk∣k,φ=∫Rnxpxk∣y1:kφxkdxk=∫Rnxgyk∣xkpxk∣y1:k−1φxkdxk∫Rnxgyk∣xkpxk∣y1:k−1dxk=πk∣k−1,gφπk∣k−1,gupdating.
Except for a small number of dynamic models, it is difficult to obtain an analytic solution in Bayesian theory (Equations (6)–(9), (11), and (12)) and the exact solution of the posterior probability for general nonlinear and non-Gaussian systems.


## 2.1. State Space Model

We setΩ,P,F as a random probability space and define two actual vector stochastic processes: X=Xk,k∈N and Y=Yk,k∈N\0, where sample space Ω is the set of all possible outcomes, event space F is a set of outcomes in the sample space, probability function P assigns a probability to each event in the event space, X is the state process, and Y is the observation information. Let nx and ny be the dimensions of the state X and observation information Y, respectively, corresponding to the state space, and define BRn as the Borelσ‐algebra set on n-dimensional Euclidean space Rn. Most nonlinear systems can take the form of a dynamic state space (DSS) [24]:
(1)Xk=FXk−1,ωkstatemodel,Yk=HXk,υkobservationequation,where k is the discrete-time (stage) index, Xk∈Rn is the set of system states at time k, and Yk∈Rn is the observation information corresponding to state Xk. F∗ and H∗ are known state transition and observation functions, respectively, corresponding to the state transition kernel PDF KXk∣Xk−1 and observation likelihood PDF gYk∣Xk in the statistical description. The system shift noise ωk∈Rn and measurement noise υk∈Rm are independent and identically distributed (i.i.d.) sequences that obey any PDF form.The state spaceX follows the first-order Markov process; i.e., the state Xk of the current moment is only related to the state Xk−1 of the previous moment. Assuming an initial distribution X0~μdX0, the probability density functions of the state transition kernel PDF KXk∣Xk−1 and observation likelihood PDF gYk∣Xk are Lebesgue measures:
(2)KXk∣Xk−1=PωXk−FXk−1,gYk∣Xk=PυYk−HXk,where Pω and Pυ are the probability functions under the influence of shift noise ωk and measurement noise υk, respectively.


## 2.2. Optimal Bayesian Recursion Theory

In Bayesian theory [15], the state Xk of a nonlinear system at time k is updated based on the observation information Y1:k to obtain a minimum mean squared error (MSE) estimate. The optimal state estimation of the system is the conditional expected value EXk∣Y1:k of the posterior PDF PXk∣Y1:k. Based on the premise that the state variable and observation function follow a first-order Markov process, the posterior PDF PXk∣Y1:k is obtained by two steps of recursive iterations: prediction and updating. Xk:l≜Xk,Xk+1,⋯,Xl and Yk:l≜Yk,Yk+1,⋯,Yl are defined as the space path information of the state process and observation likelihood from time k to l, respectively.


### 2.2.1. Prediction

Combined with the transition kernel PDFKXk∣Xk−1 calculated by the state space equation, the prior PDF PXk∣Y1:k−1 is predicted using the Chapman–Kolmogorov equation and the posterior PDF PXk−1∣Y1:k−1 at time k−1:
(3)PXkY1:k−1=∫PXk,Xk−1∣Y1:k−1dXk−1=∫PXk∣Xk−1,Y1:k−1PXk−1∣Y1:k−1dXk−1=∫KXk∣Xk−1PXk−1∣Y1:k−1dXk−1,where the state PDF PXk,Xk−1∣Y1:k−1 specifies the conditional probability of Xk∩Xk−1 given Y1:k−1, and PXk∣Xk−1,Y1:k−1 specifies the conditional probability of Xk given Xk−1∩Y1:k−1.


### 2.2.2. Updating

The prior PDFPXk∣Y1:k−1 is updated using the observation likelihood PDF gYk∣Xk at time k, and the posterior PDF PXk∣Y1:k of state Xk is obtained as
(4)PXkY1:k=PYk∣Xk,Y1:k−1PXk∣Y1:k−1PYk∣Y1:k−1=gYk∣XkPXk∣Y1:k−1∫gYk∣XkPXk∣Y1:k−1dXk,where the observation PDF PYk∣Xk,Y1:k−1 specifies the conditional probability of Yk given Xk∩Y1:k−1, and PYk∣Y1:k−1 specifies the conditional probability of Yk given Y1:k−1.The system state PDF measureπk:l∣mdXk:l is defined as
(5)πk:l∣mdXk:l≜PXk:l∣Y1:m.The joint posterior PDF measureπ0:k∣kdX0:k can be obtained by Bayes’ theorem (Equations (3) and (4)):
(6)π0:k∣k−1dX0:k=∫RnxKdXk∣Xk−1π0:k−1∣k−1dX0:k−1prediction,(7)π0:k∣kdX0:k=gYk∣Xkπ0:k∣k−1dX0:k∫RnxgYk∣Xkπ0:k∣k−1dX0:kupdating.The marginal posterior PDF measureπk∣k−1dXk is obtained by
(8)πk∣k−1dXk=∫RnxKdXk∣Xk−1πk−1∣k−1dXk−1prediction,(9)πk∣kdXk=gYk∣Xkπk∣k−1dXk∫RnxgYk∣Xkπk∣k−1dXkupdating.Similarly,xk:l≜xk,xk+1,⋯,xl and yk:l≜yk,yk+1,⋯,yl are, respectively, defined as the sampling particle path information of the state process and observation likelihood from time k to l.Definition 1.
Suppose thatν is a probability measure and φ represents an arbitrary function, Θ and Ψ are arbitrary function variables, K is the PDF of the transfer kernel satisfying the Markov process, and the following calculation method is defined:
(10)ν,φ≜∫φν,νKΘ≜∫νdxKΘ∣x,Kφx≜∫KdΨ∣xφΨ.
According to the above symbols, for any functionφ, Bayesian theory (prediction and updating processes) can be redefined, using Equation (8), as
(11)πk∣k−1,φ=∫Rnxpxk∣y1:k−1φxkdxk=∫Rnx∫RnxKxk∣xk−1pxk−1∣y1:k−1dxk−1φxkdxk=∫Rnxpxk−1∣y1:k−1∫RnxKxk∣xk−1φxkdxkdxk−1=πk−1∣k−1,Kφprediction.
From Equation (9), it is concluded that
(12)πk∣k,φ=∫Rnxpxk∣y1:kφxkdxk=∫Rnxgyk∣xkpxk∣y1:k−1φxkdxk∫Rnxgyk∣xkpxk∣y1:k−1dxk=πk∣k−1,gφπk∣k−1,gupdating.
Except for a small number of dynamic models, it is difficult to obtain an analytic solution in Bayesian theory (Equations (6)–(9), (11), and (12)) and the exact solution of the posterior probability for general nonlinear and non-Gaussian systems.


## 2.2.1. Prediction

Combined with the transition kernel PDFKXk∣Xk−1 calculated by the state space equation, the prior PDF PXk∣Y1:k−1 is predicted using the Chapman–Kolmogorov equation and the posterior PDF PXk−1∣Y1:k−1 at time k−1:
(3)PXkY1:k−1=∫PXk,Xk−1∣Y1:k−1dXk−1=∫PXk∣Xk−1,Y1:k−1PXk−1∣Y1:k−1dXk−1=∫KXk∣Xk−1PXk−1∣Y1:k−1dXk−1,where the state PDF PXk,Xk−1∣Y1:k−1 specifies the conditional probability of Xk∩Xk−1 given Y1:k−1, and PXk∣Xk−1,Y1:k−1 specifies the conditional probability of Xk given Xk−1∩Y1:k−1.


## 2.2.2. Updating

The prior PDFPXk∣Y1:k−1 is updated using the observation likelihood PDF gYk∣Xk at time k, and the posterior PDF PXk∣Y1:k of state Xk is obtained as
(4)PXkY1:k=PYk∣Xk,Y1:k−1PXk∣Y1:k−1PYk∣Y1:k−1=gYk∣XkPXk∣Y1:k−1∫gYk∣XkPXk∣Y1:k−1dXk,where the observation PDF PYk∣Xk,Y1:k−1 specifies the conditional probability of Yk given Xk∩Y1:k−1, and PYk∣Y1:k−1 specifies the conditional probability of Yk given Y1:k−1.The system state PDF measureπk:l∣mdXk:l is defined as
(5)πk:l∣mdXk:l≜PXk:l∣Y1:m.The joint posterior PDF measureπ0:k∣kdX0:k can be obtained by Bayes’ theorem (Equations (3) and (4)):
(6)π0:k∣k−1dX0:k=∫RnxKdXk∣Xk−1π0:k−1∣k−1dX0:k−1prediction,(7)π0:k∣kdX0:k=gYk∣Xkπ0:k∣k−1dX0:k∫RnxgYk∣Xkπ0:k∣k−1dX0:kupdating.The marginal posterior PDF measureπk∣k−1dXk is obtained by
(8)πk∣k−1dXk=∫RnxKdXk∣Xk−1πk−1∣k−1dXk−1prediction,(9)πk∣kdXk=gYk∣Xkπk∣k−1dXk∫RnxgYk∣Xkπk∣k−1dXkupdating.Similarly,xk:l≜xk,xk+1,⋯,xl and yk:l≜yk,yk+1,⋯,yl are, respectively, defined as the sampling particle path information of the state process and observation likelihood from time k to l.Definition 1.
Suppose thatν is a probability measure and φ represents an arbitrary function, Θ and Ψ are arbitrary function variables, K is the PDF of the transfer kernel satisfying the Markov process, and the following calculation method is defined:
(10)ν,φ≜∫φν,νKΘ≜∫νdxKΘ∣x,Kφx≜∫KdΨ∣xφΨ.
According to the above symbols, for any functionφ, Bayesian theory (prediction and updating processes) can be redefined, using Equation (8), as
(11)πk∣k−1,φ=∫Rnxpxk∣y1:k−1φxkdxk=∫Rnx∫RnxKxk∣xk−1pxk−1∣y1:k−1dxk−1φxkdxk=∫Rnxpxk−1∣y1:k−1∫RnxKxk∣xk−1φxkdxkdxk−1=πk−1∣k−1,Kφprediction.
From Equation (9), it is concluded that
(12)πk∣k,φ=∫Rnxpxk∣y1:kφxkdxk=∫Rnxgyk∣xkpxk∣y1:k−1φxkdxk∫Rnxgyk∣xkpxk∣y1:k−1dxk=πk∣k−1,gφπk∣k−1,gupdating.
Except for a small number of dynamic models, it is difficult to obtain an analytic solution in Bayesian theory (Equations (6)–(9), (11), and (12)) and the exact solution of the posterior probability for general nonlinear and non-Gaussian systems.


## 3. Particle Filter

To solve the complex problem in the above optimal Bayesian filtering algorithm, Monte Carlo sampling is used instead of an integral operation [25]. The idea is to use a discrete distribution with a series of random samples and their corresponding weights to approximate the posterior PDF measure πk∣k and calculate the expected value of the samples to estimate the actual system state Xk. The importance PDF qxk∣y1:k is generally used to represent the discrete distribution to obtain the sampling particle set x0:ki:i=1,2,⋯,N to calculate the posterior empirical measure distribution πk∣kN:
(13)πkkNdxk≜1N∑i=1Nδxkidxk,where δxki∗ is the Dirac delta function. With the sampling number N→∞, the empirical measure πk∣kN is infinitely close to the actual posterior PDF measure πk∣k.


### 3.1. Sequential Importance Resampling (SIR) Filter

Since the posterior PDF distributionpxk∣y1:k is unknown, it is necessary to construct the importance PDF qxk∣xk−1,y1:k to satisfy the requirements of the Monte Carlo sampling method and make up for the shortcoming that sampling cannot be carried out in the target distribution, and
(14)qxkxk−1,y1:k=Kxkxk−1,is typically selected during the SIR process.Assuming that the posterior PDF measureπk−1∣k−1Ndxk−1 at time k−1 is known and the particle set is xkii=1N∼πk−1k−1Ndxk−1Kdxkxk−1 at time k, the prediction measure πk∣k−1N of the prediction stage can be obtained as
(15)πkk−1Ndxk≜1N∑i=1Nδxk−xki.When the number of sample particlesN is large enough, the prediction measure πk∣k−1N is infinitely close to the actual state πk∣k−1. The Monte Carlo approximate posterior measure πk∣kN is obtained by substituting the prediction measure πk∣k−1N into Equation (9):
(16)πkkNdxk≜gykxkπkk−1Ndxk∫Rnxgykxkπkk−1Ndxk=∑i=1Ngykxkiδxkidxk∑i=1Ngykxki.The above formula is equivalent to(17)πkkNdxk=∑i=1Nωkiδxkidxk,∑i=1Nωki=1,where ωki∝gykxki is the weight of the importance PDF qxk∣y1:k after normalization of all sampling particles x0:kii=1N, and the posterior measure πk∣kN is the weighted sum of the Dirac delta function. The above process is called SIS filtering.After several updating iterations, the weights of some particles in the SIS process may be small enough to ignore, which cannot be avoided due to the shortcomings of the algorithm. To overcome this, resampling is usually used to solve the degradation problem of the standard PF algorithm. By duplicating particles with higher weights and discarding those with smaller weights, the particle setx~0:ki:i=1,2,⋯,N is gathered in the high-probability posterior region to obtain the approximate value of the unweighted empirical distribution measure π~kkN:
(18)π~kkN=1N∑i=1Nδx~kidx~k.It can be inferred that the essence of resampling is realized usingN sampling iterations in the empirical distribution measure πk∣kN, and the new particle set obtained by this method approximates the actual posterior measure πk∣k. Common resampling methods are random, system, polynomial, and residual resampling. The process of the standard PF algorithm is shown as Algorithm 1.Algorithm 1: Standard particle filter.
Step 1. Particle initializationAt timek=0, for i=1,2,⋯,NSamplex0i∼π00dx0At timek≥1, for i=1,2,⋯,NStep 2. Sequential importance samplingPredictπkk−1N, sample xki∼πk−1k−1NKdxkEvaluate the normalized importance weightsωkiωki∝gykxki, and ∑i=1Nωki=1Step 3. Resampling strategyCalculate effective samplesNeff=∑i=1Nωki2−1Compare resampling thresholdNthifNeff<Nthx~ki,ω~ki=RESAMPLExki,ωki, and ω~ki=1/Nelsex~ki,ω~ki=xki,ωkiStep 4. State estimationxk=∑i=1Nω~kix~ki


### 3.2. Clustering Similarity Particle Filter (CSPF)

The standard PF algorithm is simple in structure and easy to execute. Under the optimal estimation, the approximate estimated value of the algorithm converges to the actual state value. However, there are some issues in practical engineering applications.


#### 3.2.1. Particle Degradation Phenomenon

The standard PF introduces the importance PDF distribution in the SIS process, which causes the variance of the particle weight to accumulate with each iteration. The importance weights corresponding to most particles tend to zero, resulting in a particle degradation phenomenon [26]. The above effects lead to a significant waste of computing resources, with the result that the approximate estimation cannot accurately describe the posterior distribution of the actual state. This degradation phenomenon cannot be avoided due to defects of the algorithm.


#### 3.2.2. Particle Depletion Problem

A resampling strategy is an effective and important method to improve particle degradation. By resampling the discrete approximate posterior PDF distribution obtained by the importance sampling process, samples with larger weights are duplicated many times under the guidance of the particle motion and the distribution of the state at the previous moment so that the number of effective particles increases and degradation is suppressed. However, resampling is likely to cause the abandonment or loss of some low-weight particles, which causes the resampled particles to prematurely move away from the actual state posterior region. This results in sample dilution [27] and eventually in the increase of state estimation variance, which greatly diminishes filtering performance.In view of the above problems, our improved PF algorithm relies on the consistency principle of the spatial state trajectory [28]; i.e., the closer the state trajectory of a particle is to the actual state trajectory, the more likely the particle state represents the actual state. By using clustering similarity theory to measure the degree of trajectory consistency, the higher the degree of consistency similarity, the closer it is to the actual state, and the particle weights of the SIS process are updated to improve particle degradation. The improved algorithm abandons the resampling strategy, which can fundamentally eliminate the particle depletion problem.The particle setXk•ii=1N=xj∘i:j=k,⋯,k+L;xj′i:j=k+L+1,⋯,k+L+li=1N from time k to time k+L+l is selected as the state trajectory at time k, where L and l are predefined constants. The original trajectory set xj∘i follows the SIS filtering process [15], and the modified trajectory set xj′i complies with the GPF prediction algorithm [18]. Because the actual state is unknown, the observation likelihood information is used to represent the consistency parameter of the state trajectory. Depending on the particle state trajectory, the corresponding observed likelihood trajectory set Yk•ii=1N=yj∘i:j=k,⋯,k+L;yj′i:j=k+L+1,⋯,k+L+li=1N is determined as
(19)Yk•i=HXk•i,υk,where the measurement noise is υk=0; i.e., the observation equation H∗ is a known function determined by the specific research objects without noise interference. The observed likelihood trajectory corresponding to the actual state (reference trajectory) is Yk=yj:j=k,⋯,k+L+l. In this work, a clustering method using distance-based similarity is selected to analyze the trajectory consistency, and the distance similarity measurement dk•i [29] of the observed likelihood trajectories of the actual state and sampling particles is calculated as
(20)dk•i=disYk,Yk•ii=1N=∑j=kk+Lyj−yj∘is+∑j=k+L+1k+L+lyj−yj′iss,where dis∗ is the distance similarity function, dis∗≥0, and s is the measurement type parameter. To increase the reliability of the algorithm, the distance similarity function is transformed to an exponential similarity function:
(21)dki=eλ×dk•i,i=1,⋯,N,where λ is the gradient factor. The importance weights wk•i and wk+L+l•i corresponding to times k and k+L+l can be calculated as
(22)wk•i=dk−1i×pυYk−HXk•,υk,i=1,⋯,N,(23)wk+L+l•i=dk−1i×pυYk+L+l−HXk+L+l•,υk+L+l,i=1,⋯,N,where pυ∗ is the PDF of the observed likelihood noise. Using the above algorithm, the original trajectory set xk∘ii=1N of the SIS process and the corresponding importance weights wk•ii=1N at time k represent the posterior PDF πk∣kN of the system state, modified trajectory set xk+L+l′ii=1N of the GPF-predicted distribution, and corresponding importance weights wk+L+l•ii=1N at time k+L+l, which can approximately represent the predicted PDF πk+L+lk. Therefore, the state estimate value xk can be obtained by the filtering operation, and the state estimate xk+L+l can be calculated by the prediction step. The implementation of the improved PF algorithm is as follows.(1) Estimation. This step is consistent with the SISGPF estimation process used to extract the particle distribution set Xk•ii=1N.(2) Updating. The weights wk•ii=1N and wk+L+l•ii=1N are determined and normalized to wkii=1N and wk+L+lii=1N, respectively, to estimate and predict system states xk and xk+L+l:
(24)xk≈x~k=∑i=1Nwkixk∘ifiltering,xk+L+l≈x~k+L+l=∑i=1Nwk+L+lixk+L+l′iprediction.The above steps constitute an iterative process of the improved algorithm. Unlike the standard PF, this method uses an estimate-update-filter (prediction) process without resampling. The improved PF algorithm (CSPF) is shown as Algorithm2.Algorithm 2: Clustering similarity particle filter.
Step 1. Particle initializationAt timek=0, for i=1,2,⋯,NSamplex0i∼π00dx0At timek≥1, for i=1,2,⋯,NStep 2. Importance samplingPredictπkk−1N, sample xki∼Kxkxk−1iStep 3. Similarity measurementParticle state trajectoriesXk•iDrawxj∘i∼SIS,j=k,k+1,⋯,k+Lxj′i∼GPF,j=k+L+1,⋯,k+L+lCalculatedk•i=disYk,Yk•ii=1NExponential transformationdki=eλ×dk•iStep 4. Recursive importance weightsRevised proposal distributionqxkxk−1,yk=dkiKxkxk−1iUpdate weightswk•i∝dk−1igykxkiNormalizewki=wk•i/∑i=1Nwk•iStep 5. State estimationxk=∑i=1Nwkixki


## 3.1. Sequential Importance Resampling (SIR) Filter

Since the posterior PDF distributionpxk∣y1:k is unknown, it is necessary to construct the importance PDF qxk∣xk−1,y1:k to satisfy the requirements of the Monte Carlo sampling method and make up for the shortcoming that sampling cannot be carried out in the target distribution, and
(14)qxkxk−1,y1:k=Kxkxk−1,is typically selected during the SIR process.Assuming that the posterior PDF measureπk−1∣k−1Ndxk−1 at time k−1 is known and the particle set is xkii=1N∼πk−1k−1Ndxk−1Kdxkxk−1 at time k, the prediction measure πk∣k−1N of the prediction stage can be obtained as
(15)πkk−1Ndxk≜1N∑i=1Nδxk−xki.When the number of sample particlesN is large enough, the prediction measure πk∣k−1N is infinitely close to the actual state πk∣k−1. The Monte Carlo approximate posterior measure πk∣kN is obtained by substituting the prediction measure πk∣k−1N into Equation (9):
(16)πkkNdxk≜gykxkπkk−1Ndxk∫Rnxgykxkπkk−1Ndxk=∑i=1Ngykxkiδxkidxk∑i=1Ngykxki.The above formula is equivalent to(17)πkkNdxk=∑i=1Nωkiδxkidxk,∑i=1Nωki=1,where ωki∝gykxki is the weight of the importance PDF qxk∣y1:k after normalization of all sampling particles x0:kii=1N, and the posterior measure πk∣kN is the weighted sum of the Dirac delta function. The above process is called SIS filtering.After several updating iterations, the weights of some particles in the SIS process may be small enough to ignore, which cannot be avoided due to the shortcomings of the algorithm. To overcome this, resampling is usually used to solve the degradation problem of the standard PF algorithm. By duplicating particles with higher weights and discarding those with smaller weights, the particle setx~0:ki:i=1,2,⋯,N is gathered in the high-probability posterior region to obtain the approximate value of the unweighted empirical distribution measure π~kkN:
(18)π~kkN=1N∑i=1Nδx~kidx~k.It can be inferred that the essence of resampling is realized usingN sampling iterations in the empirical distribution measure πk∣kN, and the new particle set obtained by this method approximates the actual posterior measure πk∣k. Common resampling methods are random, system, polynomial, and residual resampling. The process of the standard PF algorithm is shown as Algorithm 1.Algorithm 1: Standard particle filter.
Step 1. Particle initializationAt timek=0, for i=1,2,⋯,NSamplex0i∼π00dx0At timek≥1, for i=1,2,⋯,NStep 2. Sequential importance samplingPredictπkk−1N, sample xki∼πk−1k−1NKdxkEvaluate the normalized importance weightsωkiωki∝gykxki, and ∑i=1Nωki=1Step 3. Resampling strategyCalculate effective samplesNeff=∑i=1Nωki2−1Compare resampling thresholdNthifNeff<Nthx~ki,ω~ki=RESAMPLExki,ωki, and ω~ki=1/Nelsex~ki,ω~ki=xki,ωkiStep 4. State estimationxk=∑i=1Nω~kix~ki


## 3.2. Clustering Similarity Particle Filter (CSPF)

The standard PF algorithm is simple in structure and easy to execute. Under the optimal estimation, the approximate estimated value of the algorithm converges to the actual state value. However, there are some issues in practical engineering applications.


### 3.2.1. Particle Degradation Phenomenon

The standard PF introduces the importance PDF distribution in the SIS process, which causes the variance of the particle weight to accumulate with each iteration. The importance weights corresponding to most particles tend to zero, resulting in a particle degradation phenomenon [26]. The above effects lead to a significant waste of computing resources, with the result that the approximate estimation cannot accurately describe the posterior distribution of the actual state. This degradation phenomenon cannot be avoided due to defects of the algorithm.


### 3.2.2. Particle Depletion Problem

A resampling strategy is an effective and important method to improve particle degradation. By resampling the discrete approximate posterior PDF distribution obtained by the importance sampling process, samples with larger weights are duplicated many times under the guidance of the particle motion and the distribution of the state at the previous moment so that the number of effective particles increases and degradation is suppressed. However, resampling is likely to cause the abandonment or loss of some low-weight particles, which causes the resampled particles to prematurely move away from the actual state posterior region. This results in sample dilution [27] and eventually in the increase of state estimation variance, which greatly diminishes filtering performance.In view of the above problems, our improved PF algorithm relies on the consistency principle of the spatial state trajectory [28]; i.e., the closer the state trajectory of a particle is to the actual state trajectory, the more likely the particle state represents the actual state. By using clustering similarity theory to measure the degree of trajectory consistency, the higher the degree of consistency similarity, the closer it is to the actual state, and the particle weights of the SIS process are updated to improve particle degradation. The improved algorithm abandons the resampling strategy, which can fundamentally eliminate the particle depletion problem.The particle setXk•ii=1N=xj∘i:j=k,⋯,k+L;xj′i:j=k+L+1,⋯,k+L+li=1N from time k to time k+L+l is selected as the state trajectory at time k, where L and l are predefined constants. The original trajectory set xj∘i follows the SIS filtering process [15], and the modified trajectory set xj′i complies with the GPF prediction algorithm [18]. Because the actual state is unknown, the observation likelihood information is used to represent the consistency parameter of the state trajectory. Depending on the particle state trajectory, the corresponding observed likelihood trajectory set Yk•ii=1N=yj∘i:j=k,⋯,k+L;yj′i:j=k+L+1,⋯,k+L+li=1N is determined as
(19)Yk•i=HXk•i,υk,where the measurement noise is υk=0; i.e., the observation equation H∗ is a known function determined by the specific research objects without noise interference. The observed likelihood trajectory corresponding to the actual state (reference trajectory) is Yk=yj:j=k,⋯,k+L+l. In this work, a clustering method using distance-based similarity is selected to analyze the trajectory consistency, and the distance similarity measurement dk•i [29] of the observed likelihood trajectories of the actual state and sampling particles is calculated as
(20)dk•i=disYk,Yk•ii=1N=∑j=kk+Lyj−yj∘is+∑j=k+L+1k+L+lyj−yj′iss,where dis∗ is the distance similarity function, dis∗≥0, and s is the measurement type parameter. To increase the reliability of the algorithm, the distance similarity function is transformed to an exponential similarity function:
(21)dki=eλ×dk•i,i=1,⋯,N,where λ is the gradient factor. The importance weights wk•i and wk+L+l•i corresponding to times k and k+L+l can be calculated as
(22)wk•i=dk−1i×pυYk−HXk•,υk,i=1,⋯,N,(23)wk+L+l•i=dk−1i×pυYk+L+l−HXk+L+l•,υk+L+l,i=1,⋯,N,where pυ∗ is the PDF of the observed likelihood noise. Using the above algorithm, the original trajectory set xk∘ii=1N of the SIS process and the corresponding importance weights wk•ii=1N at time k represent the posterior PDF πk∣kN of the system state, modified trajectory set xk+L+l′ii=1N of the GPF-predicted distribution, and corresponding importance weights wk+L+l•ii=1N at time k+L+l, which can approximately represent the predicted PDF πk+L+lk. Therefore, the state estimate value xk can be obtained by the filtering operation, and the state estimate xk+L+l can be calculated by the prediction step. The implementation of the improved PF algorithm is as follows.(1) Estimation. This step is consistent with the SISGPF estimation process used to extract the particle distribution set Xk•ii=1N.(2) Updating. The weights wk•ii=1N and wk+L+l•ii=1N are determined and normalized to wkii=1N and wk+L+lii=1N, respectively, to estimate and predict system states xk and xk+L+l:
(24)xk≈x~k=∑i=1Nwkixk∘ifiltering,xk+L+l≈x~k+L+l=∑i=1Nwk+L+lixk+L+l′iprediction.The above steps constitute an iterative process of the improved algorithm. Unlike the standard PF, this method uses an estimate-update-filter (prediction) process without resampling. The improved PF algorithm (CSPF) is shown as Algorithm2.Algorithm 2: Clustering similarity particle filter.
Step 1. Particle initializationAt timek=0, for i=1,2,⋯,NSamplex0i∼π00dx0At timek≥1, for i=1,2,⋯,NStep 2. Importance samplingPredictπkk−1N, sample xki∼Kxkxk−1iStep 3. Similarity measurementParticle state trajectoriesXk•iDrawxj∘i∼SIS,j=k,k+1,⋯,k+Lxj′i∼GPF,j=k+L+1,⋯,k+L+lCalculatedk•i=disYk,Yk•ii=1NExponential transformationdki=eλ×dk•iStep 4. Recursive importance weightsRevised proposal distributionqxkxk−1,yk=dkiKxkxk−1iUpdate weightswk•i∝dk−1igykxkiNormalizewki=wk•i/∑i=1Nwk•iStep 5. State estimationxk=∑i=1Nwkixki


## 3.2.1. Particle Degradation Phenomenon

The standard PF introduces the importance PDF distribution in the SIS process, which causes the variance of the particle weight to accumulate with each iteration. The importance weights corresponding to most particles tend to zero, resulting in a particle degradation phenomenon [26]. The above effects lead to a significant waste of computing resources, with the result that the approximate estimation cannot accurately describe the posterior distribution of the actual state. This degradation phenomenon cannot be avoided due to defects of the algorithm.


## 3.2.2. Particle Depletion Problem

A resampling strategy is an effective and important method to improve particle degradation. By resampling the discrete approximate posterior PDF distribution obtained by the importance sampling process, samples with larger weights are duplicated many times under the guidance of the particle motion and the distribution of the state at the previous moment so that the number of effective particles increases and degradation is suppressed. However, resampling is likely to cause the abandonment or loss of some low-weight particles, which causes the resampled particles to prematurely move away from the actual state posterior region. This results in sample dilution [27] and eventually in the increase of state estimation variance, which greatly diminishes filtering performance.In view of the above problems, our improved PF algorithm relies on the consistency principle of the spatial state trajectory [28]; i.e., the closer the state trajectory of a particle is to the actual state trajectory, the more likely the particle state represents the actual state. By using clustering similarity theory to measure the degree of trajectory consistency, the higher the degree of consistency similarity, the closer it is to the actual state, and the particle weights of the SIS process are updated to improve particle degradation. The improved algorithm abandons the resampling strategy, which can fundamentally eliminate the particle depletion problem.The particle setXk•ii=1N=xj∘i:j=k,⋯,k+L;xj′i:j=k+L+1,⋯,k+L+li=1N from time k to time k+L+l is selected as the state trajectory at time k, where L and l are predefined constants. The original trajectory set xj∘i follows the SIS filtering process [15], and the modified trajectory set xj′i complies with the GPF prediction algorithm [18]. Because the actual state is unknown, the observation likelihood information is used to represent the consistency parameter of the state trajectory. Depending on the particle state trajectory, the corresponding observed likelihood trajectory set Yk•ii=1N=yj∘i:j=k,⋯,k+L;yj′i:j=k+L+1,⋯,k+L+li=1N is determined as
(19)Yk•i=HXk•i,υk,where the measurement noise is υk=0; i.e., the observation equation H∗ is a known function determined by the specific research objects without noise interference. The observed likelihood trajectory corresponding to the actual state (reference trajectory) is Yk=yj:j=k,⋯,k+L+l. In this work, a clustering method using distance-based similarity is selected to analyze the trajectory consistency, and the distance similarity measurement dk•i [29] of the observed likelihood trajectories of the actual state and sampling particles is calculated as
(20)dk•i=disYk,Yk•ii=1N=∑j=kk+Lyj−yj∘is+∑j=k+L+1k+L+lyj−yj′iss,where dis∗ is the distance similarity function, dis∗≥0, and s is the measurement type parameter. To increase the reliability of the algorithm, the distance similarity function is transformed to an exponential similarity function:
(21)dki=eλ×dk•i,i=1,⋯,N,where λ is the gradient factor. The importance weights wk•i and wk+L+l•i corresponding to times k and k+L+l can be calculated as
(22)wk•i=dk−1i×pυYk−HXk•,υk,i=1,⋯,N,(23)wk+L+l•i=dk−1i×pυYk+L+l−HXk+L+l•,υk+L+l,i=1,⋯,N,where pυ∗ is the PDF of the observed likelihood noise. Using the above algorithm, the original trajectory set xk∘ii=1N of the SIS process and the corresponding importance weights wk•ii=1N at time k represent the posterior PDF πk∣kN of the system state, modified trajectory set xk+L+l′ii=1N of the GPF-predicted distribution, and corresponding importance weights wk+L+l•ii=1N at time k+L+l, which can approximately represent the predicted PDF πk+L+lk. Therefore, the state estimate value xk can be obtained by the filtering operation, and the state estimate xk+L+l can be calculated by the prediction step. The implementation of the improved PF algorithm is as follows.(1) Estimation. This step is consistent with the SISGPF estimation process used to extract the particle distribution set Xk•ii=1N.(2) Updating. The weights wk•ii=1N and wk+L+l•ii=1N are determined and normalized to wkii=1N and wk+L+lii=1N, respectively, to estimate and predict system states xk and xk+L+l:
(24)xk≈x~k=∑i=1Nwkixk∘ifiltering,xk+L+l≈x~k+L+l=∑i=1Nwk+L+lixk+L+l′iprediction.The above steps constitute an iterative process of the improved algorithm. Unlike the standard PF, this method uses an estimate-update-filter (prediction) process without resampling. The improved PF algorithm (CSPF) is shown as Algorithm2.Algorithm 2: Clustering similarity particle filter.
Step 1. Particle initializationAt timek=0, for i=1,2,⋯,NSamplex0i∼π00dx0At timek≥1, for i=1,2,⋯,NStep 2. Importance samplingPredictπkk−1N, sample xki∼Kxkxk−1iStep 3. Similarity measurementParticle state trajectoriesXk•iDrawxj∘i∼SIS,j=k,k+1,⋯,k+Lxj′i∼GPF,j=k+L+1,⋯,k+L+lCalculatedk•i=disYk,Yk•ii=1NExponential transformationdki=eλ×dk•iStep 4. Recursive importance weightsRevised proposal distributionqxkxk−1,yk=dkiKxkxk−1iUpdate weightswk•i∝dk−1igykxkiNormalizewki=wk•i/∑i=1Nwk•iStep 5. State estimationxk=∑i=1Nwkixki


## 4. Convergence Proof

The proposed algorithm is based on bootstrap filtering theory, and Bayesian state estimation can be realized by a weighted bootstrap method [15, 30]. It is assumed that the sampling particle set Xk•ii=1N is derived from a continuous PDF, Rx. The posterior PDF πkk and RxSx are proportional, and Sx is a known function. If the sample number N→∞, then the discrete distribution of particles composed of Xk•ii=1N and its corresponding weights SXk•i/∑i=1NSXk•ii=1N can be regarded as approaching the actual posterior PDF πkk. Referring to Equation (16), the posterior PDF πkk of system state xk is proportional to the product of the observation likelihood function and prior PDF πkk−1, which can be equivalent to Rx, and the weight wk•i in Equation (22) can be regarded as the observation likelihood function equivalent to Sx, which follows bootstrap filtering theory and is reasonable and effective.


### 4.1. Convergence of the Improved Algorithm

Suppose that the probability density measure space on setE=ΓRnx is the probability measure set on the largest-dimensional Euclidean space Rnx with convergence topology and set E,d is the measure space. akk=1∞ and bkk=1∞ are two continuous function sequences: E→E. In the stochastic filtering setup, space E will be all probability measure spaces on n-dimensional Euclidean space Rnx.Definition 2.
ak and bk, respectively, represent the mapping relationships of measure πk−1k−1→πkk−1 and of measure πkk−1→πkk. We define ak as the mapping relation (prediction) satisfied on measure set ΓRnx→ΓRnx:
(25)akνdxk≜νKdxk=∫RnxνKdxkxk−1dxk−1.
This holds for any measureν. Therefore, substituting the continuous function φ∈ℂbRnx in the prediction Equation (11), we obtain
(26)akν,φ=∬RnxνKdxkxk−1φxkdxk−1=ν,Kφ.
The prediction measure expression can be obtained as(27)πkk−1=akπk−1k−1.Definition 3.
Referring to Equation (12), we define bk as the mapping relation (updating) satisfied on measure set ΓRnx→ΓRnx:
(28)bkν,φ=ν,gφν,g.
The Bayesian filtering process can be expressed as(29)πkk=bkπkk−1≜bk⊙akπk−1k−1,where the operator “⊙” represents the composite mapping function.Definition 4.
Settinghk and h1:k as the conversion functions of measure πk−1k−1→πkk and of measure π00→πkk, respectively, the Bayesian filtering process can be expressed as
(30)hk≜bk⊙ak,h1:k≜hk⊙hk−1⊙⋯⊙h1.
In an abstract environment, the PF algorithm uses the Monte Carlo method to solve a problem for which it is difficult to obtain the exact analytical integral solution in Bayesian theory. The principle is to generate a series of samples from the target distribution to approximately estimate the partial characteristics of the actual state, and the estimation result is only the expectation of a “good performance” function, which can be approximated as the average value:(31)Eπφx≈φx1+φx2+⋯+φxNN=E¯φ,N.
WhenN→∞, the estimated value E¯φ,N converges to the expected value Eπφx. It can be assumed that
(32)Eπφx=μexpectation,Varφx=σ2variance,where μ is the most basic digital feature to measure the centralized position or average level of a random variable x, and σ2 is a numeric characteristic of the dispersion of the random variable x.
Based on the law of large numbers and the central limit theorem [31], it can be concluded that
(33)E¯φ,N→μ,(34)PrlimN→∞μ−E¯φ,N=0=1,where Pr∗ is the probability function.
Therefore, for the analytical solution of the integral operation, the disturbance caused by the Monte Carlo sampling method is inevitable, mainly because the estimated value is based on a random and limited sample set. However, under the guarantee of the law of large numbers and the central limit theorem, when the number of particles tends to infinity, the disturbance is minimal and satisfies the following Gaussian distributionN∗∗:
(35)E¯φ,N∼NE¯φ,Nμ,σ2/N.
WhenN→∞, the state estimate E¯φ,N converges to the real expected value μ, and the estimation variance decreases with the increase of the number of sample particles.
From the above analysis, it can be concluded that the particle filter is based on Bayesian filtering and can be combined with the Monte Carlo sampling method to generate a sampling disturbance functioncN [32]. The perturbation Equations (33) and (34) can be expressed as
(36)limN→∞cNE¯φ,N=μ.
The process formulas (29) and (30) of the particle filter algorithm can be expressed as
(37)πkkN=bk⊙cN⊙akπk−1k−1N=hkNπk−1k−1N=h1:kNκ,where κ is the initial value π00. Our improved algorithm uses clustering analysis to measure the similarity of multistage measurement information [33] as the proposed distribution qxkxk−1,yk to replace the prior PDF pxk∣xk−1 in the SIS process:
(38)qxkxk−1,yk=dkiKdxkxk−1.
The importance weight calculation is updated and modified as follows:(39)wki∝gykxkdk−1i=g⌢ykxk.
Substituting this in Equation (28), the updating formula of the improved algorithm becomes
(40)b⌢kνN,φ=νN,g⌢φνN,g⌢,where b⌢k represents the mapping relationship of the improved algorithm measure πkk−1N→πkkN and Monte Carlo measure limN→∞νN=ν. Referring to Equations (37) and (30), the improved PF can be expressed as
(41)π⌢kkN=b⌢k⊙cN⊙akπk−1k−1N=h⌢kNπk−1k−1N=h⌢1:kNκ.Theorem 5.
It is assumed that the state transition kernel functionK satisfies the first-order Markov process, and the observation likelihood function g is continuous, bounded, and strictly positive in xk∈Rnx. Under the condition of the Monte Carlo sampling disturbance cN, the improved PF algorithm measure π⌢kkN converges to the theoretical value (actual state value) πk∣k of Bayesian optimal filtering:
(42)limN→∞π⌢kkN=πkk.Proof.
In the PF algorithm, the Monte Carlo sampling disturbance is random and uncertain.cN,ϑ is set as a random disturbance, the sample number N>0, and the independent variable ϑ∈Ω. For all measures ν∈ΓRnx,
(43)cN,ϑν=1N∑j=1NδVjϑ,where Vj:Ω→Rd is an i.i.d. random variable with measure ν. According to the algorithm and the simplification of Equations (29) and (41), we can obtain
(44)πkk=bkπkk−1=bk⊙akπk−1k−1,(45)π⌢kkN=b⌢kπ⌢kk−1N=b⌢k⊙cN⊙akπ⌢k−1k−1N=b⌢k⊙cNak.
At timek, the measure of the Bayesian prediction stage is ak, and the sampling disturbance measure of the PF prediction stage is cNak. Using the i.i.d. random variables Vj, ∀φ∈Λ, and Λ=φii>0∈CbRnx, we can obtain
(46)EcNak,φi−ak,φi4=1N4E∑j=1NφiVj−ak,φi4=1N4∑j=1NEφiVj−ak,φi4=4N4∑j1,j2=1,j1≠j2NEφiVj1−ak,φi2φiVj2−ak,φi2≤N2φi4+3NN−12φi4N4≤48φi4N2,where E∗ is the solution function of set expectation, ∗ represents the supremum norm in the domain CbRnx, and φ≜supx∈Rnxφx. The summed expectation of the number of sampled particles from 1 to ∞ is
(47)E∑N=1∞cNak,φi−ak,φi4≤48φi4∑N=1∞1N2<∞.
Hence,(48)∑N=1∞cNak,φi−ak,φi4<∞.
This implies that for the prediction stage, the measure at a certain time can be expressed as(49)limN→∞cNak,φi−ak,φi=0.
Referring to Equation (40), for ν∈ΓRnx and any function φ, the updating stage measure can be obtained as
(50)limN→∞b⌢kνN,φ=limN→∞νN,g⌢φlimN→∞νN,g⌢=limN→∞dk−1iνN,gφlimN→∞dk−1iνN,g=ν,gφν,g=bkν,φ.
This result is compared with Equations (28) and (50), and it is concluded that the improved algorithm has the same measures as the Bayesian filter in the updating stage, i.e.,
(51)limN→∞b⌢kνN=bkν.
Combining Equations (44), (45), (49), and (51), we can obtain
(52)limN→∞π⌢kkN=πkkN.
Therefore, the improved PF algorithm (CSPF) based on the clustering similarity of the state trajectories still converges to the actual state under the interference of Monte Carlo sampling disturbances. Theorem5 is proved.☐


### 4.2. Convergence of the Mean Squared Error of Results

Combined with the conditions and conclusions of Section4.1, we analyze the convergence of the results by calculating the boundary of the MSE of the improved algorithm [34]. We demonstrate that the convergence of the reasoning process is related to the number of sample particles at each stage of the algorithm. Suppose that in the neighborhood of E=ΓRnx, νNwN=1∞ represents a measure sequence of random probability and satisfies limN→∞νN=ν. For any function φ∈CbRnx, it can be obtained from Theorem 5 that
(53)limN→∞EνN,φ−ν,φ2=0.Theorem 6.
It is assumed that the state transition kernel functionK satisfies the first-order Markov process, and the observation likelihood function g is continuous, bounded, and strictly positive in xk∈Rnx. For any function φ, there must be a real constant Ck∣k satisfying
(54)Eπ⌢kkN,φ−πkk,φ2≤Ckkφ2N,where φ≜supx∈Rnxφx.Proof.
According to the improved PF algorithm, the proof is divided into prediction and updating parts.☐Lemma 7.
Refer to the prediction stage in Algorithm2 (steps 1 and 2) and assume that the conditions set by Theorem 6 are met. When k≥0, there must be a real constant Ckk−1, and the prediction measure πkk−1N satisfies
(55)Eπkk−1N,φ−πkk−1,φ2≤Ckk−1φ2N.
We use induction to complete the proof. When timek=0, for any function φ, there must be a real constant C00 satisfying
(56)Eπ00N,φ−π00,φ2≤C00φ2N.
From step 1 in Algorithm2, when k=0, N i.i.d. particles are sampled from the prior PDF measure π00, x0i∼π00dx0. Using the Marcinkiewicz–Zygmund inequality [35], we obtain
(57)Eπ00N,φ−π00,φ2=E1N∑i=1Nφxi−π00,φ2=1N2E∑i=1Nφxi−Nπ00,φ2≤C00φ2N.
Thus, whenk=0, Equation (55) is proved, and π00N converges to π00.
At timek−1, for any function φ, there must be a real constant Ck−1∣k−1 satisfying
(58)Eπk−1k−1N,φ−πk−1k−1,φ2≤Ck−1k−1φ2N.
At timek, step 2 in Algorithm 2 can be derived:
(59)Eπkk−1N,φ−πkk−1,φ2≤Ckk−1φ2N.
By substituting the prediction stage formula (11) in the above formula, we obtain
(60)πkk−1N,φ−πkk−1,φ=πkk−1N,φ−πk−1k−1N,Kφ+πk−1k−1N,Kφ−πk−1k−1,Kφ≤πkk−1N,φ−πk−1k−1N,Kφ+πk−1k−1N,Kφ−πk−1k−1,Kφ.
Settingζk−1 as the σ‐field generated by particle set xkii=1N and combining this with the Monte Carlo method, we obtain
(61)Eπkk−1N,φζk−1=πk−1k−1N,Kφ.
Substitutingkφ≤φ in the above formula,
(62)Eπkk−1N,φ−πk−1k−1N,Kφ2=Eπkk−1N,φ−Eπkk−1N,φζk−12ζk−1=Eπkk−1N,φ2ζk−1−2E2πkk−1N,φζk−1+E2πkk−1N,φζk−1=Eπkk−1N,φ2ζk−1−E2πkk−1N,φζk−1≤φ2N.
Referring to Equation (58), there must be a real constant:
(63)Eπk−1k−1N,Kφ−πk−1k−1,Kφ2≤Ck−1k−1Kφ2N≤Ck−1k−1φ2N.
Using the Minkowski inequality, we obtain(64)Eπkk−1N,φ−πkk−1,φ21/2≤Eπkk−1N,φ−πk−1k−1N,Kφ21/2+Eπk−1k−1N,Kφ−πk−1k−1,Kφ21/2≤Ckk−1φN,where Ckk−1=1+Ck−1k−12. Lemma 7 is proved.Lemma 8.
Refer to the updating stage in Algorithm2 (steps 3–5) and assume that the conditions set by Theorem 6 and Lemma 7 are met. When k≥0, for any function φ, there must be a real constant Ck∣k, and the prediction measure π⌢k∣kN satisfies
(65)Eπkk−1N,φζk−1=πk−1k−1N,Kφ,where φ≜supx∈Rnxφx, setting ζk−1 as the σ‐field generated by particle set xkii=1N.
Combined with the updating stage, we can obtain the following using Equation (12):
(66)π⌢kkN,φ−πkk,φ=πkk−1N,g⌢φπkk−1N,g⌢−πkk−1,gφπkk−1,g=πkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g+πkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g.
Substituting Equation (39) in this result yields
(67)πkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g=πkk−1N,g⌢φπkk−1,g−πkk−1N,gφπkk−1N,g⌢πkk−1N,g⌢πkk−1,g=πkk−1N,dk−1gφπkk−1,g−πkk−1N,gφπkk−1N,dk−1gπkk−1N,dk−1gπkk−1,g=πkk−1N,gφπkk−1N,gπkk−1,gπkk−1,g−πkk−1N,g≤φπkk−1,gπkk−1,g−πkk−1N,g.
Similarly, using the Minkowski inequality and Equation (64), we obtain
(68)Eπ⌢kkN,φ−πkk,φ21/2≤Eπkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g21/2+Eπkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g21/2≤Eπkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g21/2+Eπkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g21/2≤φπkk−1,gEπkk−1,g−πkk−1N,g21/2+1πkk−1,gEπkk−1N,gφ−πkk−1,gφ21/2≤2gπkk−1,gCkk−1φN≤CkkφN,where 2g/πkk−1,gCkk−1=Ckk.Lemma 8 is proved, and the MSE of the improved PF algorithm is convergent.


## 4.1. Convergence of the Improved Algorithm

Suppose that the probability density measure space on setE=ΓRnx is the probability measure set on the largest-dimensional Euclidean space Rnx with convergence topology and set E,d is the measure space. akk=1∞ and bkk=1∞ are two continuous function sequences: E→E. In the stochastic filtering setup, space E will be all probability measure spaces on n-dimensional Euclidean space Rnx.Definition 2.
ak and bk, respectively, represent the mapping relationships of measure πk−1k−1→πkk−1 and of measure πkk−1→πkk. We define ak as the mapping relation (prediction) satisfied on measure set ΓRnx→ΓRnx:
(25)akνdxk≜νKdxk=∫RnxνKdxkxk−1dxk−1.
This holds for any measureν. Therefore, substituting the continuous function φ∈ℂbRnx in the prediction Equation (11), we obtain
(26)akν,φ=∬RnxνKdxkxk−1φxkdxk−1=ν,Kφ.
The prediction measure expression can be obtained as(27)πkk−1=akπk−1k−1.Definition 3.
Referring to Equation (12), we define bk as the mapping relation (updating) satisfied on measure set ΓRnx→ΓRnx:
(28)bkν,φ=ν,gφν,g.
The Bayesian filtering process can be expressed as(29)πkk=bkπkk−1≜bk⊙akπk−1k−1,where the operator “⊙” represents the composite mapping function.Definition 4.
Settinghk and h1:k as the conversion functions of measure πk−1k−1→πkk and of measure π00→πkk, respectively, the Bayesian filtering process can be expressed as
(30)hk≜bk⊙ak,h1:k≜hk⊙hk−1⊙⋯⊙h1.
In an abstract environment, the PF algorithm uses the Monte Carlo method to solve a problem for which it is difficult to obtain the exact analytical integral solution in Bayesian theory. The principle is to generate a series of samples from the target distribution to approximately estimate the partial characteristics of the actual state, and the estimation result is only the expectation of a “good performance” function, which can be approximated as the average value:(31)Eπφx≈φx1+φx2+⋯+φxNN=E¯φ,N.
WhenN→∞, the estimated value E¯φ,N converges to the expected value Eπφx. It can be assumed that
(32)Eπφx=μexpectation,Varφx=σ2variance,where μ is the most basic digital feature to measure the centralized position or average level of a random variable x, and σ2 is a numeric characteristic of the dispersion of the random variable x.
Based on the law of large numbers and the central limit theorem [31], it can be concluded that
(33)E¯φ,N→μ,(34)PrlimN→∞μ−E¯φ,N=0=1,where Pr∗ is the probability function.
Therefore, for the analytical solution of the integral operation, the disturbance caused by the Monte Carlo sampling method is inevitable, mainly because the estimated value is based on a random and limited sample set. However, under the guarantee of the law of large numbers and the central limit theorem, when the number of particles tends to infinity, the disturbance is minimal and satisfies the following Gaussian distributionN∗∗:
(35)E¯φ,N∼NE¯φ,Nμ,σ2/N.
WhenN→∞, the state estimate E¯φ,N converges to the real expected value μ, and the estimation variance decreases with the increase of the number of sample particles.
From the above analysis, it can be concluded that the particle filter is based on Bayesian filtering and can be combined with the Monte Carlo sampling method to generate a sampling disturbance functioncN [32]. The perturbation Equations (33) and (34) can be expressed as
(36)limN→∞cNE¯φ,N=μ.
The process formulas (29) and (30) of the particle filter algorithm can be expressed as
(37)πkkN=bk⊙cN⊙akπk−1k−1N=hkNπk−1k−1N=h1:kNκ,where κ is the initial value π00. Our improved algorithm uses clustering analysis to measure the similarity of multistage measurement information [33] as the proposed distribution qxkxk−1,yk to replace the prior PDF pxk∣xk−1 in the SIS process:
(38)qxkxk−1,yk=dkiKdxkxk−1.
The importance weight calculation is updated and modified as follows:(39)wki∝gykxkdk−1i=g⌢ykxk.
Substituting this in Equation (28), the updating formula of the improved algorithm becomes
(40)b⌢kνN,φ=νN,g⌢φνN,g⌢,where b⌢k represents the mapping relationship of the improved algorithm measure πkk−1N→πkkN and Monte Carlo measure limN→∞νN=ν. Referring to Equations (37) and (30), the improved PF can be expressed as
(41)π⌢kkN=b⌢k⊙cN⊙akπk−1k−1N=h⌢kNπk−1k−1N=h⌢1:kNκ.Theorem 5.
It is assumed that the state transition kernel functionK satisfies the first-order Markov process, and the observation likelihood function g is continuous, bounded, and strictly positive in xk∈Rnx. Under the condition of the Monte Carlo sampling disturbance cN, the improved PF algorithm measure π⌢kkN converges to the theoretical value (actual state value) πk∣k of Bayesian optimal filtering:
(42)limN→∞π⌢kkN=πkk.Proof.
In the PF algorithm, the Monte Carlo sampling disturbance is random and uncertain.cN,ϑ is set as a random disturbance, the sample number N>0, and the independent variable ϑ∈Ω. For all measures ν∈ΓRnx,
(43)cN,ϑν=1N∑j=1NδVjϑ,where Vj:Ω→Rd is an i.i.d. random variable with measure ν. According to the algorithm and the simplification of Equations (29) and (41), we can obtain
(44)πkk=bkπkk−1=bk⊙akπk−1k−1,(45)π⌢kkN=b⌢kπ⌢kk−1N=b⌢k⊙cN⊙akπ⌢k−1k−1N=b⌢k⊙cNak.
At timek, the measure of the Bayesian prediction stage is ak, and the sampling disturbance measure of the PF prediction stage is cNak. Using the i.i.d. random variables Vj, ∀φ∈Λ, and Λ=φii>0∈CbRnx, we can obtain
(46)EcNak,φi−ak,φi4=1N4E∑j=1NφiVj−ak,φi4=1N4∑j=1NEφiVj−ak,φi4=4N4∑j1,j2=1,j1≠j2NEφiVj1−ak,φi2φiVj2−ak,φi2≤N2φi4+3NN−12φi4N4≤48φi4N2,where E∗ is the solution function of set expectation, ∗ represents the supremum norm in the domain CbRnx, and φ≜supx∈Rnxφx. The summed expectation of the number of sampled particles from 1 to ∞ is
(47)E∑N=1∞cNak,φi−ak,φi4≤48φi4∑N=1∞1N2<∞.
Hence,(48)∑N=1∞cNak,φi−ak,φi4<∞.
This implies that for the prediction stage, the measure at a certain time can be expressed as(49)limN→∞cNak,φi−ak,φi=0.
Referring to Equation (40), for ν∈ΓRnx and any function φ, the updating stage measure can be obtained as
(50)limN→∞b⌢kνN,φ=limN→∞νN,g⌢φlimN→∞νN,g⌢=limN→∞dk−1iνN,gφlimN→∞dk−1iνN,g=ν,gφν,g=bkν,φ.
This result is compared with Equations (28) and (50), and it is concluded that the improved algorithm has the same measures as the Bayesian filter in the updating stage, i.e.,
(51)limN→∞b⌢kνN=bkν.
Combining Equations (44), (45), (49), and (51), we can obtain
(52)limN→∞π⌢kkN=πkkN.
Therefore, the improved PF algorithm (CSPF) based on the clustering similarity of the state trajectories still converges to the actual state under the interference of Monte Carlo sampling disturbances. Theorem5 is proved.☐


## 4.2. Convergence of the Mean Squared Error of Results

Combined with the conditions and conclusions of Section4.1, we analyze the convergence of the results by calculating the boundary of the MSE of the improved algorithm [34]. We demonstrate that the convergence of the reasoning process is related to the number of sample particles at each stage of the algorithm. Suppose that in the neighborhood of E=ΓRnx, νNwN=1∞ represents a measure sequence of random probability and satisfies limN→∞νN=ν. For any function φ∈CbRnx, it can be obtained from Theorem 5 that
(53)limN→∞EνN,φ−ν,φ2=0.Theorem 6.
It is assumed that the state transition kernel functionK satisfies the first-order Markov process, and the observation likelihood function g is continuous, bounded, and strictly positive in xk∈Rnx. For any function φ, there must be a real constant Ck∣k satisfying
(54)Eπ⌢kkN,φ−πkk,φ2≤Ckkφ2N,where φ≜supx∈Rnxφx.Proof.
According to the improved PF algorithm, the proof is divided into prediction and updating parts.☐Lemma 7.
Refer to the prediction stage in Algorithm2 (steps 1 and 2) and assume that the conditions set by Theorem 6 are met. When k≥0, there must be a real constant Ckk−1, and the prediction measure πkk−1N satisfies
(55)Eπkk−1N,φ−πkk−1,φ2≤Ckk−1φ2N.
We use induction to complete the proof. When timek=0, for any function φ, there must be a real constant C00 satisfying
(56)Eπ00N,φ−π00,φ2≤C00φ2N.
From step 1 in Algorithm2, when k=0, N i.i.d. particles are sampled from the prior PDF measure π00, x0i∼π00dx0. Using the Marcinkiewicz–Zygmund inequality [35], we obtain
(57)Eπ00N,φ−π00,φ2=E1N∑i=1Nφxi−π00,φ2=1N2E∑i=1Nφxi−Nπ00,φ2≤C00φ2N.
Thus, whenk=0, Equation (55) is proved, and π00N converges to π00.
At timek−1, for any function φ, there must be a real constant Ck−1∣k−1 satisfying
(58)Eπk−1k−1N,φ−πk−1k−1,φ2≤Ck−1k−1φ2N.
At timek, step 2 in Algorithm 2 can be derived:
(59)Eπkk−1N,φ−πkk−1,φ2≤Ckk−1φ2N.
By substituting the prediction stage formula (11) in the above formula, we obtain
(60)πkk−1N,φ−πkk−1,φ=πkk−1N,φ−πk−1k−1N,Kφ+πk−1k−1N,Kφ−πk−1k−1,Kφ≤πkk−1N,φ−πk−1k−1N,Kφ+πk−1k−1N,Kφ−πk−1k−1,Kφ.
Settingζk−1 as the σ‐field generated by particle set xkii=1N and combining this with the Monte Carlo method, we obtain
(61)Eπkk−1N,φζk−1=πk−1k−1N,Kφ.
Substitutingkφ≤φ in the above formula,
(62)Eπkk−1N,φ−πk−1k−1N,Kφ2=Eπkk−1N,φ−Eπkk−1N,φζk−12ζk−1=Eπkk−1N,φ2ζk−1−2E2πkk−1N,φζk−1+E2πkk−1N,φζk−1=Eπkk−1N,φ2ζk−1−E2πkk−1N,φζk−1≤φ2N.
Referring to Equation (58), there must be a real constant:
(63)Eπk−1k−1N,Kφ−πk−1k−1,Kφ2≤Ck−1k−1Kφ2N≤Ck−1k−1φ2N.
Using the Minkowski inequality, we obtain(64)Eπkk−1N,φ−πkk−1,φ21/2≤Eπkk−1N,φ−πk−1k−1N,Kφ21/2+Eπk−1k−1N,Kφ−πk−1k−1,Kφ21/2≤Ckk−1φN,where Ckk−1=1+Ck−1k−12. Lemma 7 is proved.Lemma 8.
Refer to the updating stage in Algorithm2 (steps 3–5) and assume that the conditions set by Theorem 6 and Lemma 7 are met. When k≥0, for any function φ, there must be a real constant Ck∣k, and the prediction measure π⌢k∣kN satisfies
(65)Eπkk−1N,φζk−1=πk−1k−1N,Kφ,where φ≜supx∈Rnxφx, setting ζk−1 as the σ‐field generated by particle set xkii=1N.
Combined with the updating stage, we can obtain the following using Equation (12):
(66)π⌢kkN,φ−πkk,φ=πkk−1N,g⌢φπkk−1N,g⌢−πkk−1,gφπkk−1,g=πkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g+πkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g.
Substituting Equation (39) in this result yields
(67)πkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g=πkk−1N,g⌢φπkk−1,g−πkk−1N,gφπkk−1N,g⌢πkk−1N,g⌢πkk−1,g=πkk−1N,dk−1gφπkk−1,g−πkk−1N,gφπkk−1N,dk−1gπkk−1N,dk−1gπkk−1,g=πkk−1N,gφπkk−1N,gπkk−1,gπkk−1,g−πkk−1N,g≤φπkk−1,gπkk−1,g−πkk−1N,g.
Similarly, using the Minkowski inequality and Equation (64), we obtain
(68)Eπ⌢kkN,φ−πkk,φ21/2≤Eπkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g21/2+Eπkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g21/2≤Eπkk−1N,g⌢φπkk−1N,g⌢−πkk−1N,gφπkk−1,g21/2+Eπkk−1N,gφπkk−1,g−πkk−1,gφπkk−1,g21/2≤φπkk−1,gEπkk−1,g−πkk−1N,g21/2+1πkk−1,gEπkk−1N,gφ−πkk−1,gφ21/2≤2gπkk−1,gCkk−1φN≤CkkφN,where 2g/πkk−1,gCkk−1=Ckk.Lemma 8 is proved, and the MSE of the improved PF algorithm is convergent.


## 5. Practical Applications

The factors that affect the performance of the state estimation method based on the PF are mainly determined by the construction of the spatial state model and the quality of the algorithm. The accurate architecture of the spatial state model is mainly realized by the reasonable selection of prior knowledge of the physical mechanism, as reflected in empirical physical equations, scientific parameter identification and training, multifeature search, and optimization of noise distribution. The method’s advantages and disadvantages are limited by the particle degradation caused by an unreasonable proposed distribution and particle depletion caused by the resampling strategy. Therefore, to improve and obtain the best performance of state estimation methods for a fair comparison depends on the following. (a) The architecture of the spatial state model is consistent. The same empirical physical equation and noise distribution are used in the application cases. The state tracking ability of each method is fully used, and the model parameter identification is trained. (b) The advantages of different methods are fully exploited. Commonly used and effective strategies are random, polynomial, system, and residual resampling. The applicability and accuracy of resampling strategies vary for different complex nonlinear systems. Considering the advantages of the above state estimation methods, according to the application case studied in this paper, through the analysis of simulation test results, it is concluded that the system resampling strategy performs best.Therefore, based on the performance improvements of the method, the improved method can be flexibly applied to practical application fields to solve engineering problems on the premise of application scenario guarantees, such as the accurate architecture of the state space model for the research object, scientific identification and training of physical parameters, multifeature search, and optimization of noise distribution. We compared the clustering similarity particle filter (CSPF) to several traditional filtering algorithms on two examples. All methods were implemented in MATLAB, and the root mean squared error (RMSE) was used to evaluate accuracy.Example 1.
The advantages of the improved algorithm can be seen in both the Gaussian and non-Gaussian nonlinear systems. A typical univariate nonstationary growth model (UNGM) [36] with highly nonlinear and bimodal characteristics of the state distribution is selected, which is widely used in social and economic fields such as the evaluation of urban development and short-term predictions of insurance stocks and bank interest. The following model was selected to verify the effectiveness of the CSPF algorithm:
(69)xk=αxk−1+βxk−11+xk−12+γcosδ∗k−1+ωk,yk=xk2ε+υk,where the measurement noise υk∼N0,Rk, and Rk=1. The initial state was x0=0.1, α=0.5, β=25, γ=8, δ=1.2, and ε=20. The number of particles was N=100. The simulation steps were k=1,2,⋯,T, where T=5000. The correlation coefficients were L=2 and l=1, and the gradient factor was λ=1.3. SIR, APF, GPF, and the proposed algorithm were implemented in MATLAB, and 100 trials were conducted. The computer processor speed was 2.7 GHz, and the memory capacity was 8 GB. The RMSE,
(70)RMSE=1T∑k=1Txktur−xkest,was used to evaluate the performance of the algorithms, where xktur and xkest are, respectively, the real and state estimation values of the system.The basic framework of the particle filter has state estimation and update parts. State estimation realizes state prediction of the current time based on the state prediction value of the previous time and the system state shift equation (including shift noise). Because the state variables are inevitably disturbed by various noises in the working environment, it is necessary to modify the estimated values based on actual observations to be as close to the real state values as possible so as to realize the update process of state estimation according to the observation values and measurement noise. Furthermore, the purpose of constructing the probability density function of the unknown system state using prior knowledge and actual observation data is realized. It can be seen that the shift and measurement noise will directly and greatly affect the nonlinear estimation method. This paper takes this case as an example, setting different kinds of shift noise to analyze the performance level of the nonlinear estimation methods. The two types of shift noise were Gaussian noiseωk_G∼Nωk;0,1 and non-Gaussian noise ωk_NG∼ξNωk;0,σω12+1−ξNωk;0,σω22, where ξ=0.8, σω12=1, and σω22=10. The non-Gaussian noise ωk_NG was a heavy-tailed distribution, and the acceptance-rejection sampling procedure used in Monte Carlo simulations [25] was adopted to achieve noise sampling with a confidence level of 97.5%. The simulation results are shown in Figures 1–3 and Tables 1 and 2.Figure 1
Comparison of state estimation based on different algorithms in the Gaussian noise environment. Note: SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; CSPF = clustering similarity particle filter.Figure 2
Comparison of state estimation based on different algorithms in the non-Gaussian noise environment. Note: SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; CSPF = clustering similarity particle filter.Figure 3
RMSE of SIR, APF, GPF, ECSPF, and CCSPF for 15 random simulation runs. Note: SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter.Table 1
Mean and variance of RMSE for five nonlinear algorithms running 100 times under both the Gaussian and non-Gaussian noise models.
Filter estimationSIRAPFGPFECSPFCCSPFRMSE mean (Gaussian)3.47713.43683.56771.10031.2027RMSE variance (Gaussian)0.01710.01350.02210.00250.0045RMSE mean (non-Gaussian)4.33324.39394.47261.92732.0941RMSE variance (non-Gaussian)0.02220.01550.02350.01080.0133Note: SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter; RMSE = root mean squared error.Table 2
Accuracy improvement percentage of mean and variance of RMSE for five nonlinear algorithms running 100 times under both the Gaussian and non-Gaussian noise models.
FilterRMSEGaussian noiseNon-Gaussian noiseSIRAPFGPFSIRAPFGPFECSPFMean68.4%68.0%69.2%55.5%56.1%56.9%Variance85.4%81.5%88.7%51.4%30.3%54.0%CCSPFMean65.4%65.0%66.3%51.7%52.3%53.2%Variance73.7%66.7%79.6%40.1%14.2%43.4%Note: SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter; RMSE = root mean squared error.The consistency of the state trajectory of the improved PF algorithm (clustering similarity particle filter (CSPF)) was measured by the Euclidean spatial distance (ECSPF),s=2, and Chebyshev spatial distance (CCSPF), s→∞. Figures 1 and 2 compare various PF algorithms for a single operation of the nonlinear system in Gaussian and non-Gaussian noise environments (the improved algorithm is represented by ECSPF), which can directly evaluate the effect of state estimation. Figure 3 compares the RMSE of the simulation results after 15 simulation runs. Algorithm labels ending in the letter G are denoted by solid curves and correspond to Gaussian system noise ωk_G, and those ending in the letters NG are denoted by dotted curves and correspond to non-Gaussian system noise ωk_NG. Tables 1 and 2 show the RMSE mean and variance of state estimation results after 100 simulation runs to compare the accuracy of the improved algorithm and other algorithms.Figures1 and 2 show that the state estimation accuracy of CSPF was significantly better than that of the other algorithms in the two noise environments, and the fluctuation trends of the CSPF estimated curve at any time continued moving and approaching the fluctuation trends of the actual state curve. By contrast, the other algorithms were not consistent at each time point. As shown in Figure 3 and Tables 1 and 2, the RMSE means and variances of the two improved algorithms were lower than those of SIR, APF, and GPF. For example, the prediction accuracy (variance) of CSPF was improved by 65%–69% (67%–89%) in the Gaussian noise environment and by 52%–57% (14%–54%) in the non-Gaussian noise environment. Therefore, the improved algorithm had greater accuracy and more stability, and the performance of ECSPF was particularly significant.Under the same running time, the accuracy of the state estimation of the different algorithms was determined by adjusting the number of sampling particlesN and the order of the state trajectory. From this, the efficiency of the improved algorithm was verified by comparing the corresponding operational cost (see Table 3 for details). In Gaussian and non-Gaussian noise environments, the improved algorithm with Euclidean and Chebyshev distance similarity measures was compared with SIR, APF, and GPF. For the same operational cost, i.e., for the same computing and simulation times, the improved algorithm not only used the fewest sampling particles but also had significantly better accuracy than the other three algorithms. The RMSE was reduced by factors of about 3 and 2 in Gaussian and non-Gaussian noise environments, respectively. This verified that the improved algorithm had higher operational efficiency and substantial advantages in computation time.Example 2.
We assessed the health status of a lithium-ion battery. Through the state tracking and capacity training of historical samples, the physical model parameters of the empirical degradation and distribution information of the noise were identified and optimized [3]. Different nonlinear methods were used to establish state tracking and remaining useful life (RUL) prediction models for the battery. The performances of various algorithms were evaluated based on the state tracking effect and prediction accuracy.Table 3
Computation time analysis for five nonlinear algorithms running 100 times under Gaussian and non-Gaussian noise models.
FilterGaussian noise (T=5000)FilterNon-Gaussian noise (T=1000)NRMSE meanAVG time (s)NRMSE meanAVG time (s)SIR5003.25561.5SIR3204.09042.3APF2503.32111.5APF1604.26972.3GPF8003.30931.5GPF3404.15362.3ECSPF1001.10031.5CCSPF1002.12152.3Note: AVG = average; RMSE = root mean squared error; SIR = sequential importance resampling; APF = auxiliary particle filter; GPF = Gaussian particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter.In order to solve the problem regarding battery health assessment, it is particularly important to ensure the accurate structure of the decay physical model (state space model) and scientific parameter identification; this is done under the premise of the reasonable selection of the state tracking training set (noise distribution) and initial parameters. In this paper, parameter identification is carried out based on the attenuation information of state tracking historical samples to construct the RUL prediction model [37]. Since the accuracy of state tracking and RUL prediction largely depends on the physical model of the battery capacity degradation, the development of the model requires physical knowledge of the system [38]. This is usually represented by information collected by sensors, including battery parameters (e.g., charge-discharge voltage and current, power, electrochemical impedance spectroscopy, frequency, and temperature), to build an equivalent circuit model [39] to characterize the degradation trend of battery capacity, as shown in Figure 4, which includes the open-circuit voltage (OCV), electrolyte resistance RE, polarization current IR, double-layer voltage VC, capacitance CDL, charge transfer resistance RCT, Warburg impedance RW, load current IL, and terminal voltage VL.Figure 4
Battery equivalent circuit model.Relying on the attenuation mechanism of the electrochemical characteristics, the relationship between capacity degradation and internal impedance of the battery was determined using statistical regression. The simulated attenuation characteristic of the impedance increased with the number of charge-discharge cycles to obtain a double-exponential empirical degradation model [40, 41],
(71)Cap=a∗expb∗k+c∗expd∗k,where Cap is the battery capacity A·h and k is the number of charge-discharge cycles. The unknown model parameters a and c are related to the battery impedance, b and d are related to the rate of capacity degradation, and exp∗ is the exponential function.A tracking and training model was constructed using historical samples of the battery capacity to estimate the model parameters of the empirical degradation in real time and optimize the multicharacteristic noise [42]. The identified physical model effectively converged to the gradual trend of actual battery degradation, which provided reasonable identification parameters and effective initial values to establish an RUL prediction model. This enabled the updating of the state space model composed of the state transition equation (physical model of the battery system), which can represent the recursion law of the system, and the observation equation (data characteristic measurement relationship), which can transform the implicit information of the system to visible output. The state estimation of the current time is obtained by the state transfer process and the prediction results of the previous time. Using the error between the actual measurement (noise interference) and the estimated observation information at the current moment, a weighted correction term is generated to realize the updating process. This allows the state estimation and prediction system model to be obtained with high reliability, and the process of battery state tracking and RUL predictive evaluation can be realized. A flowchart is shown in Figure 5.Figure 5
Battery state tracking and remaining useful life (RUL) prediction method for the Li-ion battery based on different algorithms.Using the state tracking process to identify the model parameters and battery RUL prediction, the SIR, APF, and CSPF algorithms were compared (since the application of the algorithm is limited by the dimension of the system variables, the GPF is not suitable for this example). State tracking parameters were identified before the charge-discharge cycleS; i.e., we estimated the optimal model parameters that can minimize the error between the predicted value of the algorithm and the actual value of the experiment. After the charge-discharge cycle S, the battery capacity was predicted to determine whether the failure threshold was exceeded. The original data of the capacity degradation of the lithium-ion battery were from the open-source experimental data of the Center for Advanced Life Cycle Engineering (CALCE) of the University of Maryland. An Arbin BT2000 battery system was the experimental platform, and the test data were stored in Excel format. A data sample with a normalized capacity of 0.90 Ah, M8, was selected as the test set for comparison of the state assessment methods. The battery was charged at the standard constant current of 0.5 C until it reached 4.2 V and was switched to constant voltage charging until the current decayed to 0.05 A, and the discharge was considered complete when the voltage dropped to 2.7 V during the charge-discharge cycle [43], as shown in Table 4. While collecting experimental data, due to the error from the accuracy of measuring equipment and human misoperation, a small amount of abnormal data was extracted from the battery dataset, as it would have affected the quality of the history dataset. The size of the dataset would also affect the operational cost, and simplifying the data (e.g., taking one point from each n point) would reduce this cost. Therefore, the original data were preprocessed, filtered, and simplified.Table 4
Experimental test conditions of CALCE dataset M8.
Battery IDExperimental platformTemperatureChargingEnd conditionDischargeEnd conditionM8Arbin BT200020-25°CConstant current/voltage4.2 V/0.05 AConstant current2.7 VGenerally, 80% of the rated capacity was taken as the end-of-life thresholdT, and the actual failure threshold was 133 cycles according to open-source data, where the state tracking set was S=65 and the number of particles was N=100. The simulation results are shown in Figure 6 and Tables 5 and 6. Figure 6 can be used to compare the effect of state tracking and the accuracy of capacity prediction of different algorithms. Tables 5 and 6 compare the state tracking and prediction results obtained by multiple simulation calculations. The evaluation indicators of the state tracking effect are the sum of squared error (SSE), MSE, RMSE, and coefficient of determination (R2). The closer SSE, MSE, and RMSE were to 0, the closer R2 was to 1, and the better was the effect of state tracking. The evaluation indicators of the prediction accuracy were the mean, variance, relative error, and RMSE of the prediction failure threshold.Figure 6
State tracking and RUL prediction results for the battery based on SIR, APF, ECSPF, and CCSPF. Note: SIR = sequential importance resampling; APF = auxiliary particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter.
(a)
SIR(b)
APF(c)
ECSPF(d)
CCSPFTable 5
Comparison of the average state tracking effect for four nonlinear algorithms running 1000 times.
FilterSSEMSERMSER2SIR1.3526×10−52.0810×10−71.4328×10−57.8318×10−1APF1.8923×10−52.9112×10−71.6408×10−56.4347×10−1ECSPF8.7368×10−71.3441×10−83.6547×10−69.8419×10−1CCSPF1.0425×10−61.6038×10−83.9975×10−69.8281×10−1Note: SIR = sequential importance resampling; APF = auxiliary particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter; SSE = sum of squared error; MSE = mean squared error; RMSE = root mean squared error;R2 = coefficient of determination.Table 6
Comparison of RUL prediction results based on four nonlinear algorithms running 100 times.
Comparison itemsSIRAPFCSPFECSPFCCSPFMean113115123122Variance156.3476.65.96.2Error15.0%13.5%7.5%8.3%RMSE23.5%26.8%10.6%11.1%Note: SIR = sequential importance resampling; APF = auxiliary particle filter; ECSPF = Euclidean distance clustering similarity particle filter; CCSPF = Chebyshev distance clustering similarity particle filter; RMSE = root mean squared error.Figure6 and Table 5 show that compared with the SIR and APF algorithms, the state tracking indicators and effect of the proposed CSPF algorithm were superior, and the Euclidean distance measure yielded the best results. The excellent state tracking effect of the CSPF method also guarantees that the attenuation mechanism of state tracking in the early period is the same as in later-stage predictions over the whole life cycle of the battery. As shown in Figure 6 and Table 6, after 100 simulations using SIR, APF, ECSPF, and CCSPF, the average failure thresholds of RUL prediction were 113, 115, 123, and 122 cycles, respectively, and the relative error was within 15.0%. Taking ECSPF as an example, the variance, relative error, and RMSE of the RUL prediction failure threshold were significantly lower than those of SIR and APF, and the prediction trend was relatively closer to the actual capacity degradation curve. The relative error of ECSPF was better than those of SIR and APF by about 50% and 45%, respectively, and the threshold variances were lower by factors of about 26 and 80, respectively. RMSE was about 55% and 60% greater, respectively. The results show that the prediction accuracy of CSPF was relatively higher than those of the other algorithms, the discrete degree of the state particle set of prediction results was the minimum, the degree of uncertainty expression was the lowest, and the algorithm was more robust.


## 6. Conclusions

Facing the difficulties of the standard PF algorithm in nonlinear system state estimation, such as low precision, instability, and low computational efficiency, we proposed an improved PF algorithm (CSPF) based on the consistency principle of the spatial state trajectory. Relying on the model construction of spatial trajectories between sampled particles and actual states, current and future multistage measurement information was predicted by SIS and GPF to form trajectories combining the original and modified trajectories. The similarity of the combined trajectories was measured by clustering analysis to guide the generation of new distributions and update the particle importance weights to mitigate the particle degradation phenomenon. Because resampling is not adopted in the improved algorithm, the problem of particle depletion is fundamentally eliminated. The convergence theorem of the improved algorithm of the CSPF and MSE of the results was proved. The effectiveness of CSPF was verified by comparison with current methods in cases of socioeconomic prediction and battery health assessment. The application experiments showed that, compared to SIR, APF, and GPF, CSPF had higher accuracy and better robustness for the state estimation of nonlinear systems under the influence of Gaussian or non-Gaussian noise. Primary conclusions are summarized as follows:(1)
In Gaussian and non-Gaussian noise environments, the prediction accuracy of CSPF was improved significantly, more than 50% in typical cases(2)
The latest measurement information was used to update the new proposal distribution, and the computational cost of CSPF was reduced by adjusting the number of sampling particles and the order of the state trajectory. In a typical case, theRMSE of CSPF was reduced by a factor of more than 2 in the same computing time(3)
The first-order Markov model was modified by the clustering similarity of state trajectories so that the indicators of state tracking effects were relatively higher, which provided an accuracy guarantee for the prediction model. In the battery application, the prediction accuracy of CSPF was the best, and the relative error of the RUL failure threshold was maintained within 8%



---

*Source: 9916339-2021-05-26.xml*