(238f) Sequential Monte Carlo Moving Horizon Estimation

State estimation is a fundamental element of engineering operations because it provides the knowledge for implementing process monitoring, control, optimization and safety related tasks. State estimation must be solved in the presence of often limited and noisy measurements as well as incomplete and uncertain process models. Many engineering systems tend to be nonlinear in nature and the uncertainties may be characterized by non-Gaussian probability density functions (pdf). The general solution to the state estimation problem is given by the marginal conditional pdf of the state from which point estimates are drawn as representative of the unknown state. Most operational tasks such as control and optimization require a point estimate for decision-making. It may not be generally claimed that the conditional mean estimate or the conditional mode estimate is superior to the other. In cases such as the existence of multiple steady states in state space, the mean estimate may lie in region of low probability between multiple peaks. In problems where a binary loss function is chosen for estimation, the mode is the only available choice. The nature of the system and the utility of the estimate will determine the choice of the point estimator.

Recently there is a surge of interest in solving the nonlinear non-Gaussian conditional mean estimation using the sequential Monte Carlo approach (SMC). This implementation of the recursive Bayesian filter is also referred to as particle filters. There has been a flurry of research activity in process systems engineering community regarding nonlinear and constrained estimation using the Monte Carlo approach with various simplifications such as Gaussianity and statistical linearization. Much of the literature in particle filtering relates to the approximation of the marginal conditional pdf and the estimation of the conditional mean. The conditional mode is difficult to estimate with small particle sets, while the mean is more robust.

There is a considerable body of literature on maximum a posteriori (MAP) state estimation, based on the maximization of the joint pdf of the state trajectory conditioned on all available measurements. The estimation of the most probable state trajectory becomes a computationally intractable problem as more measurements are accumulated and the number of decision variables increases. In practice, moving horizon estimation (MHE) limits the MAP optimization problem to a moving horizon of measurements to solve for finite sequences of the modal trajectory. For each horizon, the joint conditional pdf includes the a priori state pdf, which summarizes the uncertainty about the state at the beginning of the horizon. Traditionally, the a priori pdf is assumed to be a Gaussian pdf, which is recursively updated using linearized system dynamics such as the extended Kalman filter. In recent years, more accurate sampling based filters are suggested for determining the a priori pdf, for example, the SMC particle filters and the unscented Kalman filter. This area of research is active in what is known as the arrival cost computation in MHE. Several researchers identified potential for improving the MHE implementation using the SMC particle filter and yet others saw the synergy as the MHE providing improvements to the SMC particle filter.

In this paper the concept of approximating unknown state pdfs with Monte Carlo samples and the concept of solving a finite moving horizon MAP problem are treated in the same framework of sequential Monte Carlo. The novelty of the proposed sequential Monte Carlo Moving horizon estimation (SMCMHE) lies in (1) implementing a dynamic programming approach using the sequentially evolving particles for solving the MAP problem, specifically the Viterbi algorithm combined with the iterated dynamic programming method, which has desirable convergence properties, (2) using kernel density estimation (KDE) with fast bandwidth selection for estimating the a priori state pdf, thus eliminating the need for a separate filter running in the background for updating the arrival cost and (3) imposing multivariate linear and nonlinear constraints into the iterated dynamic programming approach.

The evolving set of particles in sequential Monte Carlo methods form a family of potential discrete state trajectories. This family of trajectories naturally generates a non-uniform and randomized dynamic grid that is adaptive according to the evolution of the state conditional pdf, which conveniently identifies the feasible regions of the state space. Dynamic programming can be effectively used to maximize the joint conditional density on this feasible grid. The dynamic randomized particle grid is computationally more favorable for dynamic programming compared to fixed grids that have to face the crippling effects of “ the menace of the expanding grids.” In this paper, the MAP trajectory identified by the Viterbi algorithm is used as an initial condition for the iterated dynamic programming method to converge to the globally optimum modal trajectory in continuous state space. The computational complexity is limited to appending computed values of the objective function with the latest measurement and dropping the earliest values of the horizon as it moves forward and the sweep searches of the Viterbi algorithm. Kernel density estimation is a nonparametric method to estimate the underlying density of a set of random numbers of particles. KDE is a convolution of the kernel function with the data, which is efficiently implemented by the fast Fourier transform (FFT). The width of the kernel, or bandwidth, must be optimally chosen for which a novel iterative method is proposed (subject of a companion paper submitted to the conference, Ungarala (2012)). The arrival cost of each horizon is computed via KDE for non-Gaussian a priori pdfs. Constraints are handled by generating importance densities over constraint surfaces in state space from which the particles are drawn. This random sampling approach is in contrast with a number of computationally demanding methods proposed in recent literature, where each particle is individually subject to constraints by further optimization methods.

The proposed approach is demonstrated with simulation examples taken from literature including benchmark nonlinear systems such as a pendulum subject to total energy conservation and a chemical reactor system subject to multivariate equality constraints on species compositions. The computational demands are compared with solutions obtained from using standard nonlinear programming code for the MAP optimization problem.

References:

Ungarala, S. (2012) Kernel density estimation using orthogonal polynomials and Goertzel DFT. Submitted to 10C07 Process Monitoring, Fault Detection and Diagnosis I, AIChE Annual Meeting, Pittsburgh, PA.