Monte Carlo Set-Membership Filtering for Nonlinear Dynamic Systems

When underlying probability density functions of nonlinear dynamic systems are unknown, the filtering problem is known to be a challenging problem. This paper attempts to make progress on this problem by proposing a new class of filtering methods in bounded noise setting via set-membership theory and Monte Carlo (boundary) sampling technique, called Monte Carlo set-membership filter. The set-membership prediction and measurement update are derived by recent convex optimization methods based on S-procedure and Schur complement. To guarantee the on-line usage, the nonlinear dynamics are linearized about the current estimate and the remainder terms are then bounded by an optimization ellipsoid, which can be described as a semi-infinite optimization problem. In general, it is an analytically intractable problem when dynamic systems are nonlinear. However, for a typical nonlinear dynamic system in target tracking, we can analytically derive some regular properties for the remainder. Moreover, based on the remainder properties and the inverse function theorem, the semi-infinite optimization problem can be efficiently solved by Monte Carlo boundary sampling technique. Compared with the particle filter, numerical examples show that when the probability density functions of noises are unknown, the performance of the Monte Carlo set-membership filter is better than that of the particle filter.


Introduction
Filtering techniques for dynamic systems are widely used in applied fields such as target tracking, signal processing, automatic control, computer vision and economics, just to name a few. The Kalman filter [1] is well known as the recursive best linear unbiased state estimator, which is clearly established as a fundamental tool for analyzing and solving a broad class of filtering problems with linear dynamic systems. When dynamic systems are nonlinear, a few well-known generalizations are the extended Kalman filter (EKF), Gaussian sum filters and unscented Kalman filtering (UKF) (see, e.g., [2,3]). These methods are based on local linear approximations of the nonlinear system where the higher order terms are ignored.
Most recently, researchers have been attracted to a new class of filtering methods based on the sequential Monte Carlo approach for nonlinear and non-Gaussian dynamic systems. Sequential Monte Carlo methods achieve the filtering task by recursively generating weighted Monte Carlo samples of the state variables by importance sampling. The samples and their weights are then used to estimate expectation, covariance and other system characteristics. The earliest two methods is the particle filter (also called the bootstrap filter) [4] and sequential imputation for general missing data problems [5]. Subsequently, a lot of methods have been developed in different situations. A sequential importance sampling framework [6] has been proposed to unify and generalize these methods. Monte Carlo filtering techniques have caught the attention of researchers in many different fields. Many excellent results in different situations can be found in, e.g., [7], [8], [9], [10], [11], and references therein. Most of these methods are based on the assumptions that probability density functions of the state noise and measurement noise are known. When underlying probability density functions (pdf) are unknown, the filtering problem for nonlinear dynamic systems is known to be a difficult problem.
Actually, when the underlying probabilistic assumptions are not realistic (e.g., the main perturbation may be deterministic), it seems more natural to assume that the state noise and measurement noise are unknown but bounded and to characterize the set of all values of the parameter or state vector that are consistent with this hypothesis [12]. The set-membership estimation was considered first at end of 1960s and early 1970s (see [13,14]). The idea of propagating bounding ellipsoids (or boxes, polytopes, simplexes, parallelotopes, and polytopes) for systems with bounded noises has also been extensively investigated, for example, see recent papers [15,16,12,17], the book [18], and references therein. Most of these methods concentrate on the linear dynamic systems.
The set-membership filtering for nonlinear dynamic systems is known to be a challenging problem. Based on ellipsoid-bounded, fuzzy-approximated or Lipschitz-like nonlinearities, several results have been made [19,20,21,22]. These results assume that the ellipsoid bounds, the coefficients of fuzzy-approximation or Lipschitz constants are known before filtering, which limit them in real-time implementation. For example, for a typical nonlinear dynamic system in a radar, the bounds of the remainder depends on the past estimates so that they cannot be obtained before filtering. As far as we know, [23] develops a nonlinear set-membership filtering which can estimate ellipsoid bounds of nonlinearities in real-time and is capable of being on-line usage, and the filter is called the extended set-membership filter (ESMF). Specifically, the nonlinear dynamics are linearized about the current estimate and the state bounding ellipsoid is relaxed to an outer bounding box by the ellipsoid projection method, the remainder terms are then bounded using interval mathematics [24], and finally the output interval box is bounded using an outer bounding ellipsoid by minimizing the volume of the bounding ellipsoid. Moreover, the set-membership filtering algorithm is derived based on the linear set-membership filtering in the earliest work [13]. It is not difficult to see that the outer bounding ellipsoids of both the remainder and the state is conservative. The cumulative effect of the conservative bounding ellipsoid at each time step may yield disconvergence of a filtering. In fact, if the state bounding ellipsoid were not relaxed to an outer bounding box by the ellipsoid projection method and using some recent linear set-membership filtering techniques [25], it should be possible to derive the tighter outer bounding ellipsoids for both the remainder and the state of the nonlinear dynamic system. More details will be clarified in Remark 4.2 and Figure 1.
In this paper, when underlying pdfs of nonlinear dynamic systems are unknown, we attempt to make progress on the corresponding filtering problem in the bounded noise setting. We propose a new class of filtering methods via set-membership estimation theory and Monte Carlo (boundary) sampling technique, denoted by MCSMF. The set-membership prediction and measurement update of MCSMF are derived by recent convex optimization methods based on S-procedure and Schur complement. To guarantee the on-line usage, the nonlinear dynamics are linearized about the current estimate and the remainder terms are then bounded by an ellipsoid, which can be described as a semi-infinite optimization problem. In general, it is an analytically intractable problem when dynamic systems are nonlinear. However, for a typical nonlinear dynamic system in target tracking, we can analytically derive some regular properties for the remainder. Moreover, based on the remainder properties and the inverse function theorem, we prove that the boundary of the remainder set must be from the the boundary of a set {||u k || ≤ 1} when we linearize the nonlinear equations by Taylor's Theorem. Thus, when we take samples from the set {||u k || ≤ 1}, the samples on the boundary {||u k || = 1} are sufficient to derive the outer bounding ellipsoids of the remainder set. The samples in {||u k || < 1} is not necessary. Therefore, the computation complexity can be reduced much more. Compared with the particle filter and ESMF in [23], numerical examples show that when the probability density functions of noises are known, the performance of the particle filter is better than that of ESMF and MCSMF. Nevertheless, when the probability density functions of noises are unknown, the performance of MCSMF is better than that of the other two filters.
The rest of the paper is organized as follows. Preliminaries are given in Section 2. In Section 3, the prediction step and the measurement update step of the set-membership filtering for nonlinear dynamic systems are derived by solving an SDP problem based on S-procedure and Schur complement, respectively. In Section 4.1, the bounding ellipsoid of the remainder set is described as a semi-infinite optimization problem and the steps of MCSMF is summarized. In Section 4.2, for a typical nonlinear dynamic system in target tracking, some regular properties for the remainder is derived. Based on the remainder properties and the inverse function theorem, the semi-infinite optimization problem can be efficiently solved by Monte Carlo boundary sampling technique. In Section 5, numerical examples are given and discussed. In Section 6, concluding remarks are provided.

Problem formulation
We consider a nonlinear dynamic system where x k ∈ R n is the state of system at time k; y k ∈ R n 1 is the measurement. f k (x k ) and h k (x k ) are nonlinear functions of x k , w k ∈ R n is the uncertain process noise and v k ∈ R n 1 is the uncertain measurement noise. They are assumed to be confined to specified ellipsoidal sets where Q k and R k are the shape matrix of the ellipsoids W k and V k , respectively, which are known symmetric positive-definite matrices. Moreover, we assume that when the nonlinear functions are linearized, the remainder terms can be bounded by an ellipsoid. Specifically, by Taylor's Theorem, f k and h k can be linearized to where J f k = ∂f k (x k ) ∂x |x k and J h k = ∂h k (x k ) ∂x |x k are Jacobian matrices, ∆f k (u k ) and ∆h k (u k ) are high-order remainders, which can be bounded in an ellipsoid for all ||u k || ≤ 1, respectively, i.e., where e f k and e h k are the center of the ellipsoids E f k and E h k , respectively; P f k and P h k are the shape matrices of the ellipsoids E f k and E h k respectively. Note that we do not assume that the ellipsoids E f k and E h k are given before filtering. Both of them are predicated in real time.
The corresponding set-membership filtering problem can be formulated as follows. Suppose that the initial state x 0 belongs to a given bounding ellipsoid: wherex 0 is the center of ellipsoid E 0 ; P 0 is the shape matrix of the ellipsoid E 0 which is a known symmetric positive-definite matrix. At time k, given that x k belongs to a current bounding ellipsoid: wherex k is the center of ellipsoid E k ; P k is a known symmetric positive-definite matrix.
The goal of the set-membership filtering is to determine a bounding ellipsoid E k+1 based on the measurement y k+1 at time k + 1, i.e, look forx k+1 , P k+1 such that the state x k+1 belongs to whenever I) x k is in E k , II) the process and measurement noises w k , v k+1 are bounded in ellipsoids, i.e.
The key problem is how to determine the bounding ellipsoids E f k and E h k in real-time so that the filtering algorithm can be on-line usage.
Moreover, we provide a state estimation ellipsoid by minimizing its "size" which is a function of the shape matrix P and is denoted by f (P ). It is well known that tr(P ) corresponds to the sum of squares of semiaxes lengths of the ellipsoid, and logdet(P ) is related to the volume of the ellipsoid. More discussion on size of the ellipsoid can be seen in [17].

Set-membership prediction and measurement update
In this section, we derive the prediction step and the measurement step of the set-membership filtering. Both of them can be converted to solve an SDP problem based on S-procedure and Schur complement. The main results are summarized to Theorems 1-2. The proofs are given in Appendix.
Proof: See Appendix.

Measurement update step
Theorem 3.2. At time k + 1, based on measurements y k+1 , the predicted bounding ellipsoid E k+1|k and the bounding ellipsoid can be obtained by solving the optimization problem in the variables P k+1 ,x k+1 , nonnegative auxiliary variables where Proof: See Appendix.
Notice that if f (P ) = tr(P ), the optimization problem in Theorems 3.1-3.2 is an SDP problem. If f (P ) = logdet(P), it is a MAXDET problem. Both of them can also be efficiently solved in polynomialtime by interior point methods for convex programming (see, e.g., [16,26]) and related softwares [27,28].

Monte Carlo Set Membership Filtering
In this section, we discuss the key problem that how to adaptively determine a bounding ellipsoid to cover the high-order remainders. In the first subsection, for the general case, the problem can be converted to solve a SDP problem via Monte Carlo sampling. Moreover, the Monte Carlo set membership filtering is presented.
In the second subsection, for target tracking, we prove that the remainder can be bounded via Monte Carlo boundary sampling. Thus, the computation complexity Algorithm 4.3 can be reduced much more.

Ellipsoid bounding of the remainder via Monte Carlo sampling
By (3)-(4), the high-order remainders are Obviously, it is a hard problem to cover a remainder by an ellipsoid since f k and h k are generally nonlinear functions. The outer bounding ellipsoid for ∆f k (u k ) is not uniquely defined, but which can be optimized by minimizing the size f (P ) of the bounding ellipsoid. Thus, the optimization problem for the bounding ellipsoid of ∆f k (u k ) can be written as where P f k = B f k B T f k , and e f k , P f k are decision variables. It is called a semi-infinite optimization problem by [29].
For a general nonlinear dynamic system, to solve the problem (27), we may use Monte Carlo sampling by uniformly taking some samples from the boundary and interior-points of the sphere ||u k || ≤ 1 so that we can get a finite set of u 1 k , . . . , u N k , then the infinite constraint (28) can be approximated by N constraints based on u 1 k , . . . , u N k . Moreover, by Schur complement, an approximate bounding ellipsoid for ∆f k (u k ) can be derived by solving the flowing SDP optimization problem: Similarly, the outer bounding ellipsoid for h k (u k ) can be derived by solving Remark 4.1. The problem (29) is an SDP problem that can be efficiently solved using modern interiorpoint methods, which have been developed by [26] and [30]. When large number of samples are required to guarantee the bounding ellipsoid contain the remainder, the one-order optimizing algorithm [31] may be used for solving the problem (29) with a lower computation complexity. In addition, in the next subsection, we will develop boundary sampling technique for target tracking, where the samples on boundary are sufficient to derive the outer bounding ellipsoids of the remainder. Thus, computation complexity can be reduced much more. Numerical examples show that only 50 uniform samples on the boundary are enough to guarantee the bounding ellipsoid contain the remainder.
Remark 4.2. Note that the bounding ellipsoid of [23] is derived by interval mathematics. We derive the bounding ellipsoid by solving a semi-infinite optimization problem. Figure 1 illustrates the difference of two methods. It is obvious to see that the bounding ellipsoid derived by solving the SDP (29) is tighter than that obtained by interval mathematics. The cumulative effect of the conservative bounding ellipsoid at each time step may yield disconvergence of a filtering.  • Step 2: (Bounding step) Take samples u 1 k , . . . , u N k from the sphere ||u k || ≤ 1, and then determine two bounding ellipsoids to cover the remainders ∆f k and ∆h k by (29)- (30) and (31)- (32), respectively.

Monte Carlo set membership filtering for target tracking
In this subsection, for a typical nonlinear dynamic system in target tracking, we discuss that the remainder can be bounded by an ellipsoid via Monte Carlo boundary sampling for target tracking. We prove that the boundary of the remainder set {∆h k+1 (u k ) : ||u k || ≤ 1} must be from the the boundary of the sphere Let us consider the following nonlinear measurement equation [2]: where x is a four-dimensional state variable that includes position and velocity (x, y,ẋ,ẏ). Note that the h(x) only depends on the first two dimensions x(1) and x(2).
We discuss the relationship between the set {||u k || ≤ 1, u k = [u k (1) u k (2)]} and the remainder set {∆h k+1 (u k ) : ||u k || ≤ 1}.  Proof. From the definition of the function g(u), it is easy to see that g(u) is a continuously differentiable function. By simple calculations, Jacobian matrix J g of g(u) is Simplifying the determinant of J h , Thus, det(J h ) ≥ 0 and the equality holds if and only if △ 2 ▽ 1 = △ 1 ▽ 2 . Since det(J g ) = det(J h )det(E) and det(E) > 0, then det(J g ) ≥ 0 and the equality holds if and only if △ 2 ▽ 1 = △ 1 ▽ 2 , at the same time, we have g(u) = 0. Moreover, by (34)-(37), it is easy to see that , and E ij is the entry of the ith row and the jth column of the matrix E. Proof. Since S 1 ⊂ S 3 , then S 1 S 2 ⊂ S 3 S 2 . Using S 1 S 2 = S 3 S 4 , we obtain S 3 S 4 ⊂ S 3 S 2 , By definition of the set S 1 , we can divide it into two parts, i.e., and c, d are defined in Lemma 4.5. According to the expression of the set S, then, we can divide the set S into the corresponding parts, i.e., . From the definition of the set S 2 1 , we can see that det(J g ) > 0 with Lemma 4.5. Using Lemma 4.7, we can find an open set W ∈ S containing g(u), in other words, z is the interior point of S, i.e., z ∈ S i , thus, S 2 ⊂ S i . According to Lemma 4.6, we can obtain S b ⊂ S 1 .
Moreover, we prove that   (500, 1000), it is easy to check that g(u) is continuously differentiable in set S 1 = {u : u ≤ 1}. We divide S 1 into three parts, i.e.,  Remark 4.9. Note that the assumption that E is a Cholesky factorization of a positive-definite P such that {x + Eu : u ≤ 1} is not intersect with the radial x(1) <= a, x(2) = b is a weak condition. If the true target is near it, we can transform the data to a new coordinate system where the target far way the the radial, then the assumption can be satisfied.

Numerical examples in target tracking
In this section, we compare the performance between Monte Carlo set membership filter and particle filter when the underlying probability density functions of noises are known or unknown. Meanwhile, we also compare it with the extended set-membership filter (ESMF) in [23].
Considering a two-dimensional Cartesian coordinate system, we track a moving target using measured range and angle from one sensor. The system equation is as follows [2]: where The x is a four-dimensional state variable that includes position and velocity (x, y,ẋ,ẏ), T = 0.2s is the time sampling interval. The process noise and measurement noise assumed to be confined to specified ellipsoidal sets The target acceleration is σ 2 = 50. In the example, the target starts at the point (50, 30) with a velocity of From the description of the above, we can see that the condition of Algorithm 4.3 is satisfied, then, using MCSMF to calculate the error bound, which is defined as follows: where x i k andx i k are the ith true state and state estimate at time k, respectively, and m is the number of the Monte Carlo runs. When the underlying probability density functions of noises are known, we use the particle filter in [9], which is denoted by PF-T. When the underlying probability density functions of noises are unknown, we denote PF-G for the particle filter where the state noise and measurement noise are assumed the truncated Gaussian noise with zero mean. At the same time, we may assume that the noises are uniform density functions, then we still use particle filter, which is denoted by PF-U. The extended set-membership filter in [23] is denoted by ESMF. These four filters have the same initial bounding ellipsoid in this example.
The following simulation results are under Matlab R2012a with YALMIP.  Figures 4-5 show that when the probability density functions of noises are known, the performance of the particle filter is better than that of MCSMF and ESMF. The reason may be that more information of the probability density of noises is used. However, when it is unknown, the performance of the particle filter is worse than that of MCSMF. In addition, the figures also show that performance of ESMF is unstable. The reason may be that there are some uncertain parameters to be used in ESMF and the remainder is bounded by interval mathematics method, which is conservative and leads a bigger bounding ellipsoid than MCSMF.  Figures 6-7 show that the 3σ confidence bounds of particle filter is indeed tighter than that of MCSMF, but it cannot contain the true state at some time step. It is an too optimistic bound. However, the bounds of MCSMF do guarantee the containment of the true state at each time step. This is useful in some applications. For example, in a civilian air traffic control system, the confidence bounds of trajectories can be used to check the standard separation between pairs of targets for maintenance of safety conditions (collision avoidance) and regularity of traffic flow in [34].
The CPU times of MCSMF, PF-T, PF-U and PF-G are plotted as a function of number of samples and particles in Figure 8, respectively. It shows that CPU times of the three filters are increasing as the number of samples and particles is increasing. The magnitude of the CPU time of the three filters are similar.

Conclusion
We have proposed a new class of filtering methods in bounded noise setting via set-membership theory and Monte Carlo (boundary) sampling technique to determine a state estimation ellipsoid. The set-membership prediction and measurement update are derived by recent convex optimization methods based on S-procedure and Schur complement. To guarantee the on-line usage, the nonlinear dynamics are linearized about the current estimate and the remainder terms are then bounded by an ellipsoid, which can be written as a semiinfinite optimization problem. For a typical nonlinear dynamic system in target tracking, based on the remainder properties and the Inverse Function Theorem, the semi-infinite optimization problem can be efficiently solved by Monte Carlo boundary sampling technique. Numerical example shows that when the probability density functions of noises are unknown, the performance of MCSMF is better than that of the particle filter,   and which is more robust than particle filter. Future work will involve, in the setting of MCSMF, the multisensor fusion, multiple target tracking and various applications such as sensor management and placement for structures and different types of wireless networks.
where I and 0 are matrices with compatible dimensions.
By S-procedure Lemma 7.1 and Eq. (51), a sufficient condition such that the inequalities (58)-(62) imply (52) to hold is that there exist scalars τ y and nonnegative scalars τ u ≥ 0, Furthermore, (63) is written in the following compact form: where Ξ is denoted by (19).
Therefore, ifx k+1|k , P k+1|k satisfy (66) and (67), then the state x k+1 belongs to E k+1|k , whenever I) x k is in E k , II) the process and measurement noises w k , v k are bounded in ellipsoidal sets, i.e., w k ∈ W k , v k ∈ V k .
Summarizing the above results, the computation of the predicted bounding ellipsoid by minimizing a size measure f (P k+1|k ) (13) is Theorem 3.1. Theorem 3.2]: Note that we have get x k+1 ∈ E k+1|k in prediction step, which is equivalent to x k+1 =x k+1|k + E k+1|k u k+1|k , u k+1|k ≤ 1, where E k+1|k is a Cholesky factorization of P k+1|k , then,
Summarizing the above results, the computation of the measurement update bounding ellipsoid by minimizing a size measure f (P k+1 ) (64) is Theorem 3.2.