Abstract
For the state estimation problem, Bayesian approach provides the most general formulation. However, most existing Bayesian estimators for dynamic systems do not take constraints into account, or rely on specific approximations. Such approximations and ignorance of constraints may reduce the accuracy of estimation. In this paper, a new methodology for the states estimation of constrained systems with nonlinear model and non-Gaussian uncertainty which are commonly encountered in practice is proposed in the framework of particles filter. The main feature of this method is that constrained problems are handled well by a sample size test and two particles handling strategies. Simulation results show that the proposed method can outperform particles filter and other two existing algorithms in terms of accuracy and computational time.
1. Introduction
With the development of advanced control theory and its application, more and more information can be obtained from most industry processes. Therefore the data-driven based approach has been used to improve the stochastic system performance and robust fault detection [1]. Bayesian methods have been studied and applied extensively in many different stochastic systems, because of their distinct ability to express uncertainties and performance optimization [2]. In practical dynamic estimation problems, nonlinear and non-Gaussian processes with constraints are commonly encountered. However, most existing Bayesian methods for such systems either ignore the constraints or rely on the assumptions of linearity and Gaussianity, such as Extended Kalman Filter [3] Moving Horizon Estimator [4]. Apparently, such ignorance of constraints and approximations of model may cause the loss of the precision of the estimates.
Sequential Monte Carlo (MC) filter that is known as particle filter (PF) [5] relaxes the assumptions of the nonlinear or the simplification of specific distribution. For the state estimation problem of nonlinear and non-Gaussian system, it provides a suboptimum solution conveniently. In many practical applications, due to the physical laws or model restrictions, the constraints of states or noise in the forms of algebraic equality and inequality are normal. There is also no question that using such constraints may lead to the improvement of estimation performance. According to such problems, an accept/reject scheme was proposed in [6]. This scheme enables the algorithm to accept the particles only on the constraint surface and all the violated particles will be discarded, ensuring all the used particles are generated from the true probability densities. However, some disadvantages are brought by this scheme. Rejecting all the particles violating constraints may reduce the number of particles and yield poor estimation. In addition, for the systems with poor prior information, all the particles may locate outside of the constraint region and cause the method failure. On this basis, paper [7] came up with a constrained PF algorithm based on hybrid use of acceptance/rejection and optimization strategies. In this algorithm, the particles outside the constrained region will be optimized only when the estimation performance based on the particles inside the region fails a chi-square test. Although the estimation performance can be improved due to this scheme, but the measurement error must be Gaussian distribution in the chi-square test which is also a similar harsh hypothesis. Meanwhile, a quantitative criterion of the particles number which is a key issue is not considered. It is believed that in some dynamic systems with poor prior information choosing the numbers of particles at pleasure will lead to a consequence that most particles will fall outside the constrained region and bringing these particles into the region requires much computational resource.
In this paper, an adaptive particle filter algorithm which chooses and optimizes particles heuristically is proposed. The novelty of this proposed method is that a sample size test is used to replace the chi-square test raised in [7]. In the sample size test, the particles outside the region will be chosen to optimize only when the number of particles in the region is less than the sampling size which can be determined on line through the Kullback-Leiler distance (KLD) [8] and other sampling size adaptive methods, such as Boers’s sampling [9] and likelihood-based sampling [10]. The constraints and particles number problem are both treated in the sample size test in order to obtain a relatively better performance with least computational cost.
The content of this paper is organized as follows. Section 2 presents the constrained state estimation problem. Section 3, introduces the principles of the Bayesian estimator and generic particle filter. The main result of this paper is given in Section 4. Section 5 provides two simulation examples. Finally, Section 6 concludes this paper.
2. Problem Statement
The dynamic of the state in discrete time is modeled by where and are nonlinear functions. , , , and are state, measurements, process noise, and measurement noise respectively. , and initial state may all follow non-Gaussian distribution. All vectors are assumed to be of appropriate dimensions. Figure 1 shows a typical example of the constraints imposed on states which occurred regularly in practice. Where the circle means the constraint in a two-dimensional space, + remarks the particles outside constrained region, with * expressing the inside particles. The main contents discussed in this paper are how to take full use of all particles to get more accurate estimation with less computation.
3. Bayesian Estimator
3.1. The Theorem of Bayesian Estimator
The objective of Bayesian estimators is to reconstruct the posteriori probability density function (PDF) with the likelihood function and priori PDF. For convenience, let and denote the full state and complete observations sets for the period of interest. represents the posterior joint PDF of the state conditioned on the measurements . Expanding it with Bayes’ rule where denote the prior PDF of the dynamic system state which is based on the prior knowledge of the model function and parameters. If the states follow the Markov property, we can rewrite as is known as the likelihood function. If the observations are conditional independent given the states , the likelihood function implies denotes the marginal PDF of the observations which can also be expressed as By (3), (4), and (5), at time , (2) becomes It is implied that the prior PDF of the state is updated by combining the likelihood function to yield the posterior PDF which are considered as a synthesis of different sources of information.
Obviously, the Bayesian filer for stochastic system involves high dimensional integration. When the system is unconstrained linear models with independent Gaussian errors, Kalman filter has been a common practice, but for constrained systems with nonlinear model and non-Gaussian uncertainty, it is impossible to get an optimal solution through Kalman filter as in linear Gaussian model, and some other suboptimal methods rely on crude approximations of the model or probability distributions are resorted to.
3.2. Generic Particle Filter
In recent years, Bayesian methods, especially Bayesian estimation of dynamic system through particle filter, have attracted attention of many researchers [11–14]. Comparing with other Bayesian estimators, particle filter is able to approximate posterior PDF by a finite number of samples and associated weights without relying on a fixed functional form of the posterior. At time , if we are able to sample random samples , from , MC estimate of this distribution would be given by where is the Dirac delta function. Unfortunately, it is unavailable to sample efficiently from the true posterior PDF. The main alternative solution to this problem is defining an importance sampling distribution [15] and generating particles from ; then the can be depicted as where the importance weight is equal to Normalized the importance weights, (9) is given by The importance weights are computed in the following form: In practice, must be selected as close as possible to the target distribution . For dynamic system, article [16] has proposed and proved that the optimal importance distribution is Hence, the important weight can be updated recursively as
It has been proved that in the entire filtering process there may exist many particles with small weights which contribute little to the estimation performance but waste most computational resource and cause the particles degeneracy problem [17]. To such problem the common solution is increasing a resampling step after the importance sampling. Thus, the estimation of the state can be expressed as where . Owing to the above analysis, the steps of the generic particle filter algorithm are summarized as follows.
Step 1. From time , sampling particles from a priori distribution .
Step 2. Sampling from importance sampling distribution , where the priori particles are , compute and normalize the weights by (13) and (11).
Step 3. In accordance with the rules of resampling, obtain the posterior particle.
Step 4. The estimation of the state can be obtained through (14) based on the posterior particles.
4. Constrained Adaptive Particle Filter
4.1. Determine the Number of Particles
It is apparent that increasing the particles number will improve the estimation performance. However, especially in real time application, blindly increasing particles does not necessarily provide greater profit to the estimated accuracy. These uncertainties are expressed in the following aspects. To increase the particles is bound to bringing greater computational burden, paying more time cost. At the same time, it may introduce more random errors which bring a negative impact to affect the estimation accuracy. Therefore, automatically determining the number of sampling particles is not only essential for reducing unnecessary computation and increasing the computational efficiency, but also useful in improving the estimation performance. In this section, the principles of determining the sampling size through Kullback-Leiler distance will be introduced. The K-L distance is defined to measure the difference between two probability distributions and : This variable is always positive and equal to zero only when the distributions are identical. Suppose the true posterior distribution is the discrete multinomial distribution. can be used to denote the number of samples drawn from different subspaces with the true probability . The maximum likelihood estimation of can be specified as , and the logarithm of the likelihood ratio is given by where then we have Let means the probability that the KLD between the true distribution and the maximum likelihood estimation base on samples is less or equal to . One can easily get the following derivation: The quantile of a chi-square distribution can be defined as: If we choose send equal to which can be written as: Combining (18) and (19), we get Approached (20) by Wilson-Hilferty transformation [18], which yields where is the upper quantile at () of the standard normal distribution. To summarize, (22) gives the number of particles that guarantees with probability that KLD is less than .
4.2. The Discard and Optimization Scheme
4.2.1. The Discard Scheme
All the particles violating the constraint region are rejected in this scheme. Let denote the constraints which can be imposed onto state or noise. The notation defined as denotes the indicator function of the region . The function is equal to unity when the state or noise at time is in the region and zero otherwise. Hence, the formula of the weights based on (13) in this scheme becomes The weights of outside particles are equal to zero when the number of particles passes the sample size test which will be introduced in the following section. With using this mechanism, all unnecessary particles will be effectively abandoned and enable us to reduce the computational resource remarkably.
4.2.2. The Optimization Scheme
By the basis of Shao et al. [7], an equation as is used to move the deviating particles to the most likely place inside the region, where , , respectively denotes the optimized particles and prior particles. and are estimated error distribution and measurement error distribution. In most application cases, truncated Gaussian or Gaussian mixtures are generally chosen to depict and [19, 20].
4.3. A Sample Size Test
In the above section, a method used to determine the particle number through the KLD has been briefly descripted. However, it must be recognized that the KLD is not the only way to choose the number of particles on line. For different requirements, different methods can be applied; one can refer to [21] for more detail.
In this section, a sample size test that combines with the generic particle filter will be introduced to handle the constraints as well as particles number problem. The particles number determined through KLD is denoted as , and then the sample size is chosen as , where is the reserve capacity determined by the user. In most situations, 1.1~1.2. The test can be expressed as follows: first, decide whether the sample particles fall into the restricted area, the statistical number of particles within the restricted area is recorded as ; if ; means that the estimation accuracy based on the particles within the region meets the requirements, particles outside the region can be discarded without hesitation. Otherwise, particles outside the region with relatively large weights are selected to optimize. The novelty of this test is that the number of particles is selected as the threshold, which determines which proposed scheme will be effective. Due to this threshold, the test not only solves the problem of determining the particles number, but also avoids the optimization of large numbers of particles which may bring with enormous computing load.
4.4. Adaptive Constrained Particle Filter
So far we have explained how to handle the constraints existing in the dynamic system by a sample size test. The adaptive constrained particle filter for dynamic system can now be summarized as follows.
Step 1. Set and generate the initial particles , , where is the initial particles number.
Step 2. Consider . To meet the estimation accuracy defined by user, (22) and other sample size adaptive algorithms can be used to determine the number of particles within the regions. We choose the true sampling size as , and then generate prior particles from the important sampling distribution.
Step 3. Evaluate the importance weights for every particles using (24) and then normalize the importance weights.
Step 4. If the number of the particles inside the region passes the sample size test, the discard scheme is executed to reject all the outside particles, then jump to Step 5. Otherwise, jump to Step 5 after the optimization scheme.
Step 5. In order to avoid the particles degeneracy problem, some resampling strategies such as multiply/suppress samples with high/low normalized important weights are employed. Based on these strategies, generate the posterior particles.
Step 6. The estimation of the state based on particles is given by (14). Repeat Step 2~Step 6 until reaching the end of the filter process.
It is notable that sometimes in order to ensure the accuracy of the filtering, a minimum particles number is employed. When the number of particles calculated by the algorithms described in Section 4.1 is less than , then will be replaced by .
5. Simulation Result
To test the efficacy of our method, two examples are given in this section. The simulations were carried out in Matlab 2009a on a Pentium IV PC with 2.5 GHz CPU and 1 G memory. The root mean square error (RMSE) is defined to assess the performance of the algorithm; the smaller the value, the higher the estimation accuracy. The true states and estimated states in time and th are, respectively, corresponding to and . Thus, the RMSE based on times Monte-Carlo simulations is characterized as The adaptive constrained particle filter (ACPF), generic particle filter (GPF), Lang’s particle filter (LPF) [6], and constrained particle filter (CPF) proposed by [7] were used to estimate the states and simulated for 30 times with 50 time steps.
Example 1. The first scenario is a typical nonlinear system originally taken from [22], where the following system model is considered: where As in [23], impose a nonnegative constraint on the process noise. is zero mean Gaussian noise with covariance and the measurement noise . For the methods of fixed sampling size, 300 particles are used. The particles number of ACPF is updated by KLD method with , confidence interval , and initial particles number and .
It is clear to see from Figure 2 that compared to GPF and CPF, the average quantity of samples in ACPF can be significantly reduced.
From Table 1 and Figure 3, CPF and ACPF yield improvements in estimation performance significantly over GPF due to considering the constraints in the estimation process which is consistent with the analysis. At the same time, relative to CPF, the CPU time of ACPF can be notably reduced for almost the same level of performance which is attributed to the sample size test.
Example 2. In order to compare conveniently, consider a bath reactor system as where . The states are constrained according to Define the initial condition ; the measurement noise and the process noise are zero mean Gaussian noise with covariance matrix The initial guess is far from the true value which is expressed as means poor prior information of this system. LPF and CPF are simulated with the sample size .
From the simulation results we found, LPF is inappropriate for this system due to the stringency of the restricted area which, leading to the number of particles, is almost equal to zero. Figures 4, 5, 6, and 7 compare the estimations of three different algorithms. It is revealed that since without considering the constraints, the performance of GPF is worse than ACPF and CPF. From Table 2, under the same number of particles, the RMSE of CPF is smaller than GPF, but the incremental accuracy is at the expense of a large amount of computing time, as also clearly reflected in Table 2. It shows that the ACPF overcomes the drawbacks of CPF with using a sample size test to get better estimation performance with less time cost.
As summarized, results of the two examples are confirmed that compared with GPF, the existence of two kinds of particles processing mechanisms results in the improvement of the estimation accuracy and robustness of the poor prior information. The using of sample size test in ACPF makes the algorithm avoid numerous computational loads brought by the optimization mechanism in the poor priori circumstance. According to the two simulations, the proposed algorithm reaches the comparative estimation accuracy of CPF with only 1/4 of its computation time, making the new method more suitable for real-time application.
6. Conclusion
In this paper, a practical approach to particle filter based on Bayesian estimation to the estimation process of dynamic system is described. The key idea of this method is that some adaptive sample size algorithms and two particles handling strategies are simultaneously used to solve the constrained problem and improve the performance of estimation. The simulation results show that the proposed method overcomes the drawback of computational complexity by a sample size test. Using the discard and optimization scheme not only ensures all the samplings representing the true distribution which contributes to estimation performance, but also increases the robustness of this algorithm to the poor prior information. Another significant advantage of this algorithm is that different adaptive sample size algorithms can be chosen for different requirements of estimation accuracy and computing consumption. For the constrained dynamic systems, the practical value of the raised method is enormous.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank the financial support by the National Natural Science Foundation of China (nos. 21206053, 21276111, and 61273131), and this study was partly supported by the 111 Project (B12018), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).