Performance Assessment of Nonlinear State Filters

8th IFAC Symposium on Advanced Control of Chemical Processes The International Federation of Automatic Control Singapore, July 10-13, 2012

Performance Assessment of Nonlinear State Filters A. Tulsyan ∗ , B. Huang ∗ , R.B Gopaluni ∗∗ , J.F Forbes ∗ ∗ Computer Process Control Group, Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Canada, T6G-2G6 (e-mail: {tulsyan; biao.huang; fraser.forbes}@ ualberta.ca) ∗∗ Process Modeling and Control Lab, Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, Canada, V6T-1Z3 (e-mail: [email protected])

Abstract: Nonlinear state filters of different approximations and capabilities have been developed in the last decade. The quality of different nonlinear filters, in terms of the mean squared error (MSE) of the estimates, depends on the approximations used in the filtering algorithm; however, there are no known methods for effectively evaluating the relative performance of these filters. A new method which measures the performance of different state filters against the theoretical posterior Cram´er-Rao lower bound (PCRLB) is proposed. The complex high-dimensional integrals in PCRLB are approximated using sequential Monte-Carlo (SMC) methods. Efficacy of the proposed method is illustrated through a simulation example. 1. INTRODUCTION In process industries, nonlinear filtering methods find wide applications in model based control and monitoring. Under the Bayesian framework, a filtering problem aims at constructing a posterior filter density at each filtering time (Doucet et al. [2001]). In general, nonlinear filters require infinite number of parameters to analytically represent the posterior filter density (Ristic et al. [2004]). In the last few decades, several tractable filtering solutions based on analytical and statistical approximation techniques (e.g., extended Kalman filters (EKFs) and particle filters (PFs), respectively) have been developed for state estimation in nonlinear systems (Arulampalam et al. [2002]). Popular filters such as EKFs and PFs are efficient in estimating in stochastic nonlinear systems; however, their performance is often limited or affected by various numerical and statistical approximations. Performance assessment of these sub-optimal, but tractable filters is, therefore, important for: (i) assessing the quality of the approximations involved; and (ii) comparing the filtering performance against that of an optimal filter. Despite the great practical interest in evaluating performance of various sub-optimal filters, it still remains one of the most complex problems ˇ in estimation theory (Simandl et al. [2001]). Selecting an appropriate benchmark is crucial for assessment of different filtering algorithms. The conventional Cram´er-Rao lower bound (CRLB), defined as an inverse of the Fisher information matrix (FIM) provides a lower bound on the second order error (MSE) obtained with any maximum-likelihood (ML) based unbiased state or parameter estimator. Analogous extension of CRLB to the class of Bayesian derived estimators, (of which, filtering is a sub-class) was derived by Trees [1968], which is commonly referred to as PCRLB. Note that a full statistical characterization of any non-Gaussian density 978-3-902823-05-2/12/$20.00 © 2012 IFAC

371

require higher-order moments, in addition to the mean and covariance. As a result, PCRLB does not fully characterize the accuracy of nonlinear filtering algorithms. Nevertheless, PCRLB is an important tool, as it only depends on the fundamental properties of the process, such as: the process dynamics; prior knowledge about the initial condition of the states; and system noise characteristics (Bergman [2001]). Furthermore, PCRLB is independent of any particular realization of the input-output data (see Section 2.2), thereby motivating its use as a benchmark. Computation of the PCRLB for nonlinear filters, as derived by Trees [1968], is based on batch data, which often renders it impractical for systems with a large number of states. To overcome the issue of high computational requirements, several recursive forms of the PCRLB were developed. Recently, Tichavsk´ y et al. [1998] proposed a recursive approach to compute the PCRLB in general nonlinear state-space models with non-Gaussian noise. Although the formulation in Tichavsk´ y et al. [1998] allows one to recursively compute the lower bound; however, obtaining a tractable solution for the same is non-trivial. This is due to the complex multi-dimensional integrals involved that are not amenable to any easy analytical solution. Several attempts have been made in the past to approximate the PCRLB. In Bergman [2001], an SMC method is used to approximate the PCRLB for systems described by linear state dynamics with Gaussian noise and nonˇ linear measurement model. In Simandl et al. [2001], the PCRLB approximation for linear and nonlinear models with additive Gaussian process and measurement noise is reported. To the author’s bet knowledge, there are no tractable solution to approximate the PCRLB for general nonlinear state-space models with non-Gaussian noise. This paper deals with performance assessment of various nonlinear state filters for a class of processes described by general nonlinear state-space equations with non-Gaussian 10.3182/20120710-4-SG-2026.00046

8th IFAC Symposium on Advanced Control of Chemical Processes Singapore, July 10-13, 2012

noise. The proposed approach compares the MSE obtained with different nonlinear filtering methods against the theoretical PCRLB, for which an SMC based recursive approximation algorithm is developed. 2. PROBLEM FORMULATION 2.1 Stochastic Nonlinear system Consider the following discrete-time stochastic nonlinear state-space model: xt+1 = ft (xt , ut , θ, t, vt ), (1a) yt = gt (xt , ut , θ, t, wt ), (1b) where: xt ⊆ X ∈ Rn×1 and yt ⊆ Y ∈ Rm×1 are the state and measurement variables, respectively; ut ⊆ U ∈ Rp×1 are the time-varying control variables; θ ⊆ Θ ∈ Rr×1 are the model parameters; ft (·) is a n-dimensional state mapping function; and gt (·) is a m-dimensional output mapping function, each being nonlinear in xt and θ, such that ft := X × U × Θ × [1, T ] × Rn×1 → Rn×1 and gt := X × U × Θ × [1, T ] × Rm×1 → Rm×1 , where T is the final filtering time. The state and measurement noise are represented as vt ∈ Rn×1 and wt ∈ Rm×1 , respectively. The model structure in (1) represents a class of general stochastic nonlinear state-space models. For notational simplicity, explicit dependency on ut and θ are not shown in the rest of this article. Assumption 1. vt and wt are mutually independent and i.i.d sequences from the probability density functions (pdf) px (·) and py (·), respectively. px (·) and py (·) are known (e.g., Gaussian; uniform) a priori and are parametrized in finite number of moments (e.g., mean; variance). 2.2 PCRLB: Preliminaries A lower bound on the MSE of the estimated states at time t in (1) is given by (Trees [1968], Tichavsk´ y et al. [1998])   T bt|t )(xt − x bt|t ) ≥ Jt−1 , (2) Pt|t , Ep(x0:t ,y1:t ) (xt − x where: Pt|t is a n×n matrix, x bt|t is an unbiased estimate of xt for a sub-optimal filter, given a measurement sequence y1:t = {y1 , . . . , yt }; Jt is a n × n posterior FIM for the optimal filter; p(x0:t , y1:t ) is a joint probability distribution for the states and measurements constructed based on information available until time t; the superscript (·)T is the transpose operator; and Ep(·) [·] is the expectation operator with respect to p(·). Alternatively, the lower bound in (2) can also be written in the scalar form as:   Ep(x0:t ,y1:t ) (xt − x bt|t )T (xt − x bt|t ) ≥ tr[Jt−1 ], (3) where tr[·] is the trace operator. For n = 1, (2) and (3) are the same. The PCRLB representation in (2) implies that the difference Pt|t − Jt−1 is a positive semi-definite matrix for all t ∈ [1, T ]. If x bt|t in (2) or (3) is selected as the minimum MSE (MMSE) estimates, such that x bt|t , Ep(xt |y1:t ) [xt ] (Trees [1968]) then   Pt|t =Ep(x0:t ,y1:t ) (xt − x bt|t )(xt − x bt|t )T , (4a)   T =Ep(y1:t ) Ep(x0:t |y1:t ) (xt − x bt|t )(xt − x bt|t ) , (4b)   =Ep(y1:t ) Ep(xt |y1:t ) (xt − x bt|t )(xt − x bt|t )T , (4c) =Ep(y1:t ) Vp(xt |y1:t ) [xt ] , (4d) 372

where V[·] is the covariance operator. From (4a) and (4b), p(x0:t , y1:t ) = p(x0:t |y1:t )p(y1:t ) is used; and from (4b) and (4c) the probability marginalization property is used. Finally, Pt|t in (4d) can be understood as the expected conditional covariance of xt ∼ p(xt |y1:t ) distributed according to the posterior distribution p(xt |y1:t ), over all possible realizations of the measurements y1:t ∼ p(y1:t ). ˇ In Tichavsk´ y et al. [1998] and Simandl et al. [2001], a recursive method to compute Jt in (2) or (3) is given by Jt+1 = Dt22 − Dt21 (Jt + Dt11 )−1 Dt12 where: J0 = − Ep(x0 ) [∆xx00 log p(x0 )]; Dt11 Dt21 Dt12 Dt22

=− =− =− =− −

∀t ∈ [1, T ],

(5) (6a)

Ep(x0:t+1 ,y1:t+1 ) [∆xxtt log p(xt+1 |xt )]; Ep(x0:t+1 ,y1:t+1 ) [∆xxtt+1 log p(xt+1 |xt )]; Ep(x0:t+1 ,y1:t+1 ) [∆xxtt+1 log p(xt+1 |xt )]; t+1 Ep(x0:t+1 ,y1:t+1 ) [∆xxt+1 log p(xt+1 |xt )] xt+1 Ep(x0:t+1 ,y1:t+1 ) [∆xt+1 log p(yt+1 |xt+1 )];

(6b) (6c) (6d) (6e)

and: ∆ is a second order partial derivative operator such ∂ that ∆yx , ∇x ∇Ty with ∇x , ∂x ; p(xt+1 |xt ) and p(yt |xt ) model the process and the measurement noise pdf, respectively; and p(x0 ) is the pdf of the states xt at t = 0. The expectation operators in (6b) through (6e) are defined with respect to p(x0:t+1 , y1:t+1 ), which make the matrices {Dtij ; i = {1, 2}; j = {1, 2}} independent of any particular realization of the states or the measurements. Also, note that in the above recursion, Dt12 = [Dt21 ]T . See Tichavsk´ y et al. [1998] for detailed derivation of the recursion in (5). From (6), it is evident that (5) depends only on the system dynamics, the process and measurement noise distributions and the choice of p(x0 ). In fact, the posterior FIM in (5) is independent of the choice of any filtering algorithm. Thus (5) provides a true lower bound on the MSE or the expected conditional covariance of the posterior distribution. This motivates its use as a benchmark for performance assessment of different sub-optimal filters. 3. APPROXIMATING PCRLB As discussed in Section 1, computing an analytical solution to the posterior FIM in (5) for the model and noise form considered in (1) and Assumption 1, respectively, is non-trivial. Obtaining a tractable solution to Jt in (5) requires approximating the following expected values of the derivatives of the pdfs: It11 ,Ep(x0:t+1 |y1:t+1 ) [∆xxtt log p(xt+1 |xt )];

It21 It22,a It22,b

,Ep(x0:t+1 |y1:t+1 ) [∆xxtt+1 t+1 ,Ep(x0:t+1 |y1:t+1 ) [∆xxt+1 t+1 ,Ep(x0:t+1 ,y1:t+1 ) [∆xxt+1

(7a)

log p(xt+1 |xt )];

(7b)

log p(xt+1 |xt )];

(7c)

log p(yt+1 |xt+1 )].

(7d) Dt11 = [Dt21 ]T ;

Using the same argument as in (4a) and (4b): −Ep(y1:t+1 ) [It11 ]; Dt21 = −Ep(y1:t+1 ) [It21 ]; Dt12 = and Dt22 = −Ep(y1:t+1 ) [It22,a ] − It22,b . As discussed in Section 1, SMC method will be used to approximate the integrals in (7). The SMC approach developed in Robert and Casella [1999] is based on N random sampling of the state trajectories {xi0:t+1 ; i = 1, . . . , N } distributed

8th IFAC Symposium on Advanced Control of Chemical Processes Singapore, July 10-13, 2012

according to the pdf p(x0:t+1 |y1:t+1 ). For example, SMC approximation of It11 based on this approach is given by It11 = Ep(x0:t+1 |y1:t+1 ) [∆xxtt log p(xt+1 |xt )] ≈

N 1 X xt [∆ log p(xit+1 |xit )]. (8) N i=1 xt

As discussed in Andrieu et al. [2004], the approximation in (8) is not efficient as the likelihood of occurrence of the sampled state trajectories xi0:t+1 are not considered. An approach suggested by Andrieu et al. [2004] approximates the integral It11 as Ep(x0:t+1 |y1:t+1 ) [∆xxtt log p(xt+1 |xt )] ≈ N 1 X i xt w [∆xt log p(xit+1 |xit )], N i=1

(9)

where wi are the weights corresponding to the sampled trajectories xi0:t+1 . As pointed in Andrieu et al. [2004] and Gopaluni [2008], the approximation in (9) requires sampling from a high dimensional pdf p(x0:t+1 |y1:t ), which can be inefficient as t increases. An alternate approach presented in Gopaluni [2008] requires sampling from a low dimensional pdf. In the following section, similar lower dimensional approximation approach will be developed for the integrals in (7).

weights proportional to occurrence of xit|t and xit+1|t+1 , respectively; and δ(·) is the Dirac delta function. Now, substituting (12a) and (12b) into (11), an empirical estimate of p(xt:t+1 |y1:t+1 ) is obtained as PN i p(xt+1 |xt ) i=1 wt|t δ(xt − xit|t ) pb(xt:t+1 |y1:t+1 ) = R PN m δ(x − xm )dx p(xt+1 |xt ) m=1 wt|t t t t|t xt ×

X

X

× dxt:t+1 , Z = [∆xxtt log p(xt+1 |xt )]p(xt:t+1 |y1:t+1 )dxt:t+1 . X

(10)

To approximate It11 in (10), random samples from the joint distribution p(xt:t+1 |y1:t+1 ) is required. Unfortunately, direct sampling from p(xt:t+1 |y1:t+1 ) is non-trivial. To overcome this, p(xt:t+1 |y1:t+1 ) in (10) is simplified using Bayes’ theorem such that (see Appendix A for the derivation) p(xt+1 |xt )p(xt |y1:t )p(xt+1 |y1:t+1 ) R p(xt:t+1 |y1:t+1 ) = . p(xt+1 |xt )p(xt |y1:t )dxt xt (11) Using SMC estimates of p(xt |y1:t ) and p(xt+1 |y1:t+1 ) based on the procedure outlined in Sch¨ on et al. [2011] yields, p(xt |y1:t ) ≈ pb(xt |y1:t ) =

N X i=1

p(xt+1 |y1:t+1 ) ≈ pb(xt+1 |y1:t+1 ) =

N X

i wt|t δ(xt − xit|t ), (12a)

j wt+1|t+1 δ(xt+1 − xjt+1|t+1 ),

(12b)

j=1

where: xit|t and xit+1|t+1 are the samples of the states drawn from the posterior distributions p(xt |y1:t ) and i i p(xt+1 |y1:t+1 ), respectively; wt|t and wt+1|t+1 are the 373

j wt+1|k+1 δ(xt+1 − xjt+1|k+1 ).

j=1

(13)

Several algebraic manipulations of (13) yield pb(xt:t+1 |y1:t+1 ) = PN PN j,i j j i i j=1 i=1 wt+1|t+1 wt|t p(xt+1|t+1 |xt|t )δ(zt+1,t − zt+1,t ) , PN m m m=1 wt|t p(xt+1 |xt|t ) =

N X N X

j,i j,i ), δ(zt+1,t − zt+1,t wt+1,t

(14)

j=1 i=1

where:

j,i zt+1,t , [xjt+1 ; xit ];

zt+1,t , [xt+1 ; xt ]; j,i wt+1,t ,

3.1 Approximating matrix: It11 Using the definition of expectation and from the Markov property of the model in (1), the dimensionality of the integration operator in (7a) can be reduced as follows Z It 11 = [∆xxtt log p(xt+1 |xt )]p(x0:t+1 |y1:t+1 )dx0:t+1 , X Z  Z xt = [∆xt log p(xt+1 |xt )] p(x0:t+1 |y1:t+1 )dx0:t−1

N X

j i wt+1|k+1 wt|t p(xjt+1|k+1 |xit|k ) . PN j m m m=1 wt|t p(xt+1 |xt|t )

(14) provides an approximate SMC estimate for the required pdf p(xt:t+1 |y1:t+1 ), which can be made arbitrarily accurate by increasing N . The computational complexity j,i of the weights wt+1,t in (14) is of the order O(N 2 ). As suggested in Gopaluni [2008], complexity of order O(N ) is obtained by replacing (14) with the following pb(xt:t+1 |y1:t+1 ) =

where:

i,i wt+1,t

N X

(15)

i=1

i,i ζt+1,t

= PN

i,i ζt+1,t =

i,i i,i wt,t+1 δ(zt+1,t − zt+1,t ),

;

j,j j=1 ζt+1,t i i wt+1|k+1 wt|t p(xit+1|t+1 |xit|t ) . PN m i m m=1 wt|t p(xt+1 |xt|t )

Now, substituting (15) into (10), we obtain Z  xt  Ibt11 = ∆xt log p(xt+1 |xt ) pb(xt:t+1 |y1:t+1 )dxt:t+1 , X

=

N X i=1

 xt  i,i ∆xt log p(xt+1 |xt ) {zt+1,t =zi,i wt+1,t

t+1,t

}

. (16)

b 11 = −Ep(y The SMC estimate of Dt11 is given as D [Ibt11 ]. t 1:t+1 ) The expectation over measurements y1:t+1 can again be approximated using the SMC method as given below pb(y1:t+1 ) =

M 1 X q δ(y1:t+1 − y1:t+1 ), M q=1

(17)

q where M is the total number of output trajectories y1:t+1 simulated randomly using (1). Substituting (17) into the b t11 , yields expression for the matrix estimate D

8th IFAC Symposium on Advanced Control of Chemical Processes Singapore, July 10-13, 2012

M X N X b t11 = − 1 D wi,i,q × M q=1 i=1 t+1,t

[∆xxtt log p(xt+1 |xt )]{zt+1,t =zi,i,q } , t+1,t

(18)

i,i,q i,i,q are the particles and correspondand wt+1,t where: zt+1,t ing weights given in (15) based on the output trajectory q . Finally, (18) provides an SMC approximation to the y1:t+1 matrix Dt11 used in (6b). Similar SMC approximation for other matrices in (7b) through (7d) are given in Appendix B. Finally, substituting the matrix approximations in (18); (B.2a); (B.2b); and (B.5) into (5), yields an estimate of the posterior FIM such that b t22 − D b t21 (Jbt + D b 11 )−1 D b 12 . Jbt+1 = D (19) t

t

It is further emphasized that the computation of the approximate PCRLB solution as required in (2) or (3) is straightforward. Applying the matrix inversion Lemma (see Horn and Johnson [1985]) on (19) yields, b t22 ]−1 − [D b t22 ]−1 D b t21 × [Jbt+1 ]−1 = [D h i−1 b t12 [D b t22 ]−1 D b t21 − (Jbt + D b t11 ) b t12 [D b t22 ]−1 . (20) D D The expression in (20) can thus be used recursively approximate the theoretical PCRLB. 4. FINAL ALGORITHM In this section, a formal procedure for approximating the recursive PCRLB in general state-space models is outlined. Algorithm 1: PCRLB Approximation

q (1) Generate a set of random trajectories of y1:T ∀q ∈ [1, M ] by simulating (1). (2) Compute N particles with corresponding weights i,q {xi,q t|t ; wt|t ; i = 1, . . . , N }∀t ∈ [1, T ] based on each q output trajectory y1:T ∀q ∈ [1, M ] using any standard particle filtering algorithm. See Arulampalam et al. [2002] for available options. i,q (3) Resample N particles {xi,q t|t ; wt|t ; i = 1, . . . , N }∀t ∈ [1, T ]; ∀q ∈ [1, M ] using any standard resampling algorithm. See Ristic et al. [2004] for available options. (4) Compute J0 based on the pdf p(x0 ) selected a priori. −1 If p(x0 ) is Gaussian then J0 = P0|0 . Set t ← 1. (5) while t 6= T do Set q ← 1. (6) while q 6= M do b t11 ;D b t21 ;D b t12 and; (7) Compute the MC estimate of D 22 b using the resampled particles. Set q ← q + 1. D t (8) end while (9) Compute the MC estimate of the FIM (Jbt ) and the PCRLB (Jbt−1 ) using (19) and (20), respectively. Set t ← t + 1. (10) end while

It is important to emphasize that the proposed approach makes multiple approximations in deriving a tractable solution and therefore, the quality of the solution depends on the accuracy of the approximations involved. Also, due to multiple approximations, the positive semi-definite condition on Pbt|t − Jbt−1 is not guaranteed to hold for all t ∈ [1, T ]. It is beyond the scope of this work to validate the quality of the bound approximation; however, it should 374

be noted that all the particle approximations used in this work can be made arbitrarily accurate simply by increasing the number of particles (N ) and the MC simulations (M ). The choice of N and M are user defined, which can be selected based on the required bound accuracy and available computing speed. It is also emphasized that an appropriate choice of the resampling method (see Step 3, Algorithm 1) can also significantly improve the quality of the approximate solution in Algorithm 1. 5. SIMULATION EXAMPLE In this section, a simulation example is used for assessing the performance of two popular nonlinear filters: sequential-importance-resampling (SIR) filter (Doucet et al. [2001]) and EKF filter (Ristic et al. [2004]). Consider the following nonlinear model (Sch¨on et al. [2011]) xt 25xt xt+1 = + + 8 cos(1.2t) + vt , 2 1 + x2t yt = 0.05x2t + wt , (21) where: vt and wt are i.i.d processes with Gaussian distributions such that vt ∼ N (vt |0, Q) and wt ∼ N (wt |0, R). The algorithm parameters are selected as T = 20 seconds; N = 1000 particles; and M = 200. In this simulation study, a nonlinear system with Gaussian noise is considered for illustrative purposes; however, it should be noted that the proposed algorithm is far more general and can handle non-Gaussian noise. Performance assessment of SIR and EKF filters are considered for the following two cases (1) Case I: The variance of the process and measurement noise are selected as Q = 0.01 and R = 0.01. (2) Case II: The variance of the process and measurement noise are selected as Q = 0.1 and R = 0.1. In both Case I and II, measure of the expected conditional covariance of the posterior distribution (Pt|t ∀t ∈ [1, T ], see Section 2.2) computed using SIR filter and EKF are compared against the PCRLB solution. It is highlighted that obtaining an analytical solution to Pt|t ∀t ∈ [1, T ] is non-trivial for the model form considered in (21). Note that, EKF filter recursively computes an estimate of the conditional covariance of the posterior distribution at each filtering time which when averaged q over all output trajectories y1:T ∀q ∈ [1, M ] yields an estimate of Pt|t ∀t ∈ [1, T ]. See Ristic et al. [2004] for further details on EKF and its implementation. In SIR filter, an estimate of Pt|t ∀t ∈ [1, T ] can be obtained using SMC approach such that,   Ep(x ,y ) (xt − x bt|t )T (xt − x bt|t ) ≈ Pbt|t 0:t

1:t

=

M N 1 X X i,q i,q w (x − x bqt|t )T (xi,q bqt|t ), (22) t|t − x M q=1 i=1 t|t t|t

i,q where: x bqt|t is the MMSE estimate; and {xi,q t|t ; wt|t } are the set of particles and weights computed using SIR filter q based on the output trajectory y1:t . See Doucet et al. [2001] for further details on computing (22).

For simulation, prior distribution for both the filters are defined as x0 ∼ N (x0 |0, 5). Larger prior variance ensures that the true initial state is defined within the chosen prior.

8th IFAC Symposium on Advanced Control of Chemical Processes Singapore, July 10-13, 2012

EKF, the SIR filter quickly converges to the true state. The case is best illustrated in Fig.2(a) for t ≤ 7.

30 True SIR

20

EKF

In Fig.1(b) and Fig.2(b), trajectories of Pbt|t are compared against PCRLB solution computed using Algorithm 1. As expected, the PCRLB is lower than the trajectories of Pbt|t computed based on SIR and EKF filters. Furthermore, the computed PCRLB is positive for all t ∈ [1, T ]. Even though the validation of the quality of the PCRLB solution is not included in the scope of this work, the simulation results reaffirms the practical reliability of Algorithm 1.

x bt|t

10

0

−10 −20

2

4

6

8

10 12 (a) Time

0.06

14

16

18

20

PCRLB SIR EKF

Pbt|t

0.04

0.02

0

4

6

8

10

12 (b) Time

14

16

18

20

Fig. 1. (a) State trajectory estimation based on 200 MC simulations.(b) Comparing the PCRLB against the trajectory estimates of Pt|t (t = 1, 2 are not shown due to scaling issues). The trajectories are computed based on SIR and EKF filters under Case I.

Interestingly, in Fig.2 (at t = 7, 8, for instance), compared to the SIR filter, the EKF algorithm delivers a relatively lower Pbt|t but poorer point estimate for the state. The result reflects bias issues with the EKF under poor choice of the prior distribution and strong system nonlinearities.

20 True SIR

10 x bt|t

EKF

0

−10 −20

2

4

6

8

10 12 (a) Time

0.4

14

16

18

Finally, in Fig.3, the PCRLB computed for Case I and II are compared against each other. As one would expect, the theoretical lower bound for Case II is relatively higher than that for Case I. The results in Fig.3 highlights the estimation difficulties associated with higher noise levels, irrespective of the choice of a nonlinear filter.

20 PCRLB SIR

0.3 Pbt|t

EKF

0.2 0.1 0

4

6

8

10

12 (b) Time

14

16

18

In Fig.1(b), the trajectories of Pbt|t computed using the SIR and the EKF filter almost perfectly match the PCRLB solution except at certain filtering time. This is due to the low noise variance considered in Case I. For Case II, in Fig.2(b), the trajectories of Pbt|t are larger for both the filters compared to Case I but are still in the neighbourhood of the PCRLB solution. Note that the performance of the SIR filter for both Case I and II can be improved further by employing an efficient resampling algorithm and by increasing the number of particles (N ).

20

Fig. 2. (a) State trajectory estimation based on 200 MC simulations.(b) Comparing the PCRLB against the trajectory estimates of Pt|t (t = 1, 2 are not shown due to scaling issues). The trajectories are computed based on SIR and EKF filters under Case II.

The results in this simulation study appear promising in our attempt to develop a tool for performance assessment of different nonlinear state filters; however, it should be noted that the recursive PCRLB is a function of time and system inputs. In such situations, the performance assessment of filtering algorithms based on point-by-point comparison with the PCRLB at each filtering time might not be reliable. To circumvent this, in future, a time and input independent performance index will be considered for quantifying the overall performance of different nonlinear state filtering algorithms.

0.25 PCRLB−I PCRLB−II

6. CONCLUSIONS

0.2

PCRLB

0.15

0.1

0.05

0

2

4

6

8

10

12

14

16

18

20

Time

Fig. 3. Comparing the PCRLB solution computed using Algorithm 1 for Case I and II. In Fig.1(a) (Case I) and Fig.2(a) (Case II), the trajectory of the true state is compared against the state estimates computed using EKF and SIR filters. Compared to the 375

A performance assessment method for different nonlinear state filters based on PCRLB is proposed for a general stochastic nonlinear state-space models with non-Gaussian noise. The proposed approach provides a tractable solution for computing the theoretical PCRLB. Simulation results suggest that the computed PCRLB solution is reasonably accurate considering the multiple approximations involved. Estimation difficulties associated with any filtering algorithm under higher process and measurement noise is also illustrated through a simulation example. ACKNOWLEDGEMENT This work was supported by the Natural Sciences and Engineering Research Council (NSERC), Canada.

8th IFAC Symposium on Advanced Control of Chemical Processes Singapore, July 10-13, 2012

Appendix A. JOINT DISTRIBUTION

Now, to obtain an MC estimate of It22,b , substituting (12b) and (17) into (B.3) yields

Using the definition of conditional probability and by appealing to the Bayes’ theorem several times, we obtain: p(xt:t+1 |y1:t+1 ) = p(xt |xt+1 , y1:t , yt+1 )p(xt+1 |y1:t+1 ), p(yt+1 |xt , xt+1 , y1:t )p(xt |xt+1 , y1:t )p(xt+1 |y1:t+1 ) = , p(yt+1 |xt+1 , y1:t ) p(yt+1 |xt+1 , y1:t )p(xt |xt+1 , y1:t )p(xt+1 |y1:t+1 ) = , p(yt+1 |xt+1 , y1:t ) =p(xt |xt+1 , y1:t )p(xt+1 |y1:t+1 ), p(xt+1 |xt , y1:t )p(xt |y1:t )p(xt+1 |y1:t+1 ) = , p(xt+1 |y1:t ) p(xt+1 |xt )p(xt |y1:t )p(xt+1 |y1:t+1 ) R . (A.1) = p(xt+1 |xt )p(xt |y1:t )dxt xt

With A.1, the joint pdf p(xt:t+1 |y1:t+1 ) can be expressed as a product of filtering densities and state transition density. Appendix B. OTHER APPROXIMATIONS

Ibt21 = Ibt22,a =

 xt+1  i,i wt+1,t ∆xt log p(xt+1 |xt ) {z

N X

h i i,i wt+1,t ∆xxt+1 log p(xt+1 |xt ) t+1

i=1

i=1

i,i t+1,t =zt+1,t }

.

(B.1a)

i,i {zt+1,t =zt+1,t }

.

(B.1b) Dt21

Using (B.1a), SMC estimate for the matrices and Dt12 in (6c) and (6d), respectively, are computed as given below N M X X b t21 = − 1 D wi,i,q × M q=1 i=1 t+1,t

[∆xxt+1 log p(xt+1 |xt )]{zt+1,t =zi,i,q t

b 12 =[D b 21 ]T D t t

t+1,t

}

(B.2a) (B.2b)

Approximating the matrix It22,b in (7d) requires a slightly different approach than one considered in Section 3.1. x This is because the integrand ∆xt+1 t+1 log p(yt+1 |xt+1 )] is no longer independent of p(y1:t+1 ). Therefore, using the Markov property of (1), (7d) can be simplified as follows It22,b = − Ep(x0:t+1 ,y1:t+1 ) [∆xxt+1 log p(yt+1 |xt+1 )], t+1 Z =− [∆xxt+1 log p(yt+1 |xt+1 )]p(x0:t+1 |y1:t+1 ) t+1 X ,Y

× p(y1:t+1 )dx0:t+1 dy1:t+1 , Z =− [∆xxt+1 log p(yt+1 |xt+1 )]p(y1:t+1 )dz1:t+1 t+1 X ,Y Z  × p(x0:t+1 |y1:t+1 )dx0:t dxt+1 , Z X =− [∆xxt+1 log p(yt+1 |xt+1 )]p(xt+1 |y1:t+1 ) t+1 X ,Y

× p(y1:t+1 )dxt+1 dy1:t+1 .

N

q {yt+1 =yt+1 ; xt+1 =xi,q } t+1

(B.3) 376

, (B.4)

i,q where: xi,q t+1|t+1 and wt+1|t+1 are the particles and corresponding weights, respectively, from the posterior distribuq q tion p(xt+1 |y1:t+1 ) based on the output trajectory y1:t+1 . 22,a b t22 = −Ep(y b An estimate of Dt22 is given as D [ I ] t 1:t+1 ) 22,b b −It , which can be obtained by substituting (B.1b) and b t22 such that, (B.4) into the expression for D b 22 = D t

−

M N 1 X X h i,i,q [∆xt+1 log p(xt+1 |xt )]{zt+1,t =zi,i,q } + w t+1,t M q=1 i=1 t+1,t xt+1 i i,q t+1 [∆xxt+1 log p(yt+1 |xt+1 )]{yt+1 =yq ; xt+1 =xi,q } . wt+1|t+1 t+1

Analogous to the approximation method adopted in Section 3.1, It21 and It22,a in (7b) and (7c), respectively, can be approximated as follows N X

M

1 X X i,q Ibt 22,b = w × M q=1 i=1 t+1|t+1 i h log p(y |x ) ∆xxt+1 t+1 t+1 t+1

t+1

(B.5)

REFERENCES C. Andrieu, A. Doucet, S. S. Singh, and V. B. Tadic. Particle methods for change detection, system identification, and control. Proceedings of IEEE, 92(3):423–438, 2004. M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp. A tutorial on particle filters for online nonlinear/nongaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50(2):174–188, 2002. N. Bergman. Sequential Monte Carlo Methods in Practice, chapter Posterior Cramer-Rao ´ Bounds for Sequential Estimation. New York: Springer–Verlag, 2001. A. Doucet, J. F. G. de Freitas, and N.J N. J. Gordon. Sequential Monte Carlo Methods in Practice, chapter On Sequential Monte Carlo Methods. New York: Springer– Verlag, 2001. R. B. Gopaluni. A particle filter approach to identification of nonlinear processes under missing observations. The Canadian Journal of Chemical Engineering, 86(6):1081– 1092, 2008. R. A. Horn and C. R. Johnson. Matrix Analysis. New York: Cambridge Univ. Press, 1985. B. Ristic, S. Arulampalam, and N. Gordon. Beyond the Kalman Filter: Particle Filters for Tracking Applications. Boston: Artech House, 2004. C. P. Robert and G. Casella. Monte Carlo Statistical Methods. New York: Springer–Verlag, 1999. T. B. Sch¨on, A. Wills, and B. Ninness. System identification of nonlinear state-space models. Automatica, 47: 39–49, 2011. P. Tichavsk´ y, C. Muravchik, and A. Nehorai. Posterior Cram´er-Rao bounds for discrete-time nonlinear filtering. IEEE Transactions on Signal Processing, 46:1386–1396, 1998. H. L. V. Trees. Detection, Estimation and Modulation Theory. New York: Wiley, 1968. ˇ M. Simandl, J. Kr´alovec, and P. Tichavsk´ y. Filtering, predictive and smoothing Cr´amer-Rao bounds for discretetime nonlinear dynamic systems. Automatica, 37:1703– 1716, 2001.