LMMSE Filtering in Feedback Systems with White Random Modes: Application to Tracking in Clutter
aa r X i v : . [ c s . I T ] A p r SUBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 1
LMMSE Filtering in Feedback Systems with WhiteRandom Modes: Application to Tracking in Clutter
Daniel Sigalov, Tomer Michaeli, andYaakov Oshman,
Fellow, IEEE
Abstract —A generalized state space representation of dynamical sys-tems with random modes switching according to a white random processis presented. The new formulation includes a term, in the dynamicsequation, that depends on the most recent linear minimum mean squarederror (LMMSE) estimate of the state. This can model the behavior ofa feedback control system featuring a state estimator. The measurementequation is allowed to depend on the previous LMMSE estimate of thestate, which can represent the fact that measurements are obtained froma validation window centered about the predicted measurement and notfrom the entire surveillance region. The LMMSE filter is derived forthe considered problem. The approach is demonstrated in the context oftarget tracking in clutter and is shown to be competitive with severalpopular nonlinear methods.
I. I
NTRODUCTION
State estimation in dynamical systems with randomly switchingcoefficients is an important problem in many applications. Naturalexamples are maneuvering target tracking and fault detection andisolation algorithms, featured, e.g., in aerospace navigation systems.In the standard modeling the dynamics of the continuously-valuedstate, and, possibly, its measurement equation, are controlled by adiscrete evolving mode. This is the well known concept of hybridsystems [1].Various problems have been formulated using the hybrid systemsframework. In problems involving uncertain observations, such as [2],[3], the mode affects the matrices of the measurement equation. Intarget tracking applications, considered in, e.g., [4]–[6], the modeusually affects the dynamics equation.We consider a state space representation of dynamical systemswith random coefficients that constitute a white stochastic sequence,accompanied by the following feedback terms. First, we allow thesystem input to depend on the latest estimate of the state, as iscommon practice in closed loop control systems. In this work, thestate estimate is taken to be the linear minimum mean squared error(LMMSE) estimate. In addition, the measurement equation is also setto depend on the latest LMMSE state estimate. This can representthe fact that observations are not taken in the entire feasible space,but, rather, in a small validation window set about the predictedmeasurement of the state.It is well known [5] that, even for the case of independentlyswitching modes, the optimal estimate of the state cannot be obtainedwithout resorting to exhaustive enumeration. Therefore, significantefforts have been dedicated to developing suboptimal approaches forstate estimation in hybrid systems and especially for the importantsubclass of jump linear systems (JLS). The most popular nonlinearmethods include the generalized pseudo-Bayesian (GPB) filter [5]and the interacting multiple model (IMM) algorithm [6]. Alterna-tively, one may consider optimality within the narrower family oflinear filters. Among these we mention [2] and [3] that consideredestimation with uncertain observations, [7] that derived a Kalmanfilter-like (KF) algorithm for a JLS with independently switchingmodes and uncorrelated matrices within each time step, and [8] thatderived an LMMSE scheme for a Markov JLS by means of stateaugmentation. In addition, in some cases, parts of the state may beestimated optimally while others in a linear optimal manner, as wasshown in [9].In this paper we concentrate on feedback JLS with independentmode transitions and consider optimal estimation within the family of linear filters. We derive a recursive LMMSE algorithm that maybe conveniently implemented in a recursive form, eliminating theneed for unbounded memory. Unlike [7], we do not assume thatthe matrices within each time step are uncorrelated. This allowstackling a wider variety of problems, such as tracking in clutter, whichcannot be modeled directly within the framework of [7]. On the otherhand, since we still treat the easier case of independent, rather thanMarkov, mode transitions, we do not require state augmentation, asdoes the algorithm of [8]. Our filter reduces to several previouslyreported results when the parameters of the underlying problem areappropriately adjusted. As an illustration, we formulate the problemof target tracking in clutter within the proposed framework and showthat the resulting filter is competitive with several classical nonlinearmethods.The paper is organized as follows. In Sec. II we describe theproposed modeling and survey some related work. The recursiveLMMSE algorithm is derived in Sec. III. An application to targettracking in clutter, followed by a numerical study, is presented inSec. IV. Concluding remarks are given in Sec. V.II. S
YSTEM M ODEL AND R ELATED W ORK
We consider the dynamical system x k +1 = A k x k + B k u k + C k w k (1a) y k = H k x k + G k v k + F k ˆ x k − , (1b)where x k ∈ R n and y k ∈ R m are the state and measurement vectorsat time k , respectively. The processes { w k } and { v k } constitute zero-mean unity-covariance strictly white sequences, and x is a randomvector (RV) with mean ¯ x and second-order moment P .We consider two variants for the modeling of u k . In the first case, u k is a known deterministic input. However, because in some cases u k serves as a closed loop control signal, it is common practice tolet it depend on the most recent estimate of the state. Thus, in thesecond variant we set u k = ˆ x k , where ˆ x k is the LMMSE estimateof x k using the measurement history Y k , { y , . . . , y k } .Likewise, the term ˆ x k − in the measurement equation is theLMMSE estimate of x k − based on the measurement history Y k − .Affecting the measurement at time k , the term F k ˆ x k − can be used torepresent the fact that observations are not taken in the entire space,but, rather, in a small validation window, set about the predictedmeasurement.The system mode, M k , { A k , B k , C k , H k , G k , F k } , is a strictlywhite random process with known distribution. The quantities { w k } , { v k } , {M k } , and x are assumed to be independent.We seek to obtain the LMMSE estimate ˆ x k +1 using the measure-ments Y k +1 . It will be shown in the sequel that, in our setting, ˆ x k +1 conveniently possesses the recursive form ˆ x k +1 = L k ˆ x k + K k y k +1 + J k u k (2)thus avoiding the need to store the entire measurement sequence.When u k = ˆ x k , the terms L k ˆ x k and J k ˆ x k in (2) may be groupedtogether.Note that the described problem does not require the system modeto assume values in a discrete domain as opposed to, e.g. [2], [3],[8]. In addition, the above formulation allows evolution not only ofthe entries of the mode matrices, but also of their dimensions [10].This observation allows treatment of problems that, to the best ofour knowledge, have not been previously considered in the contextof LMMSE algorithms. One such example is given in Section IV.For the setting without feedback terms, several variants and specialcases of the presented problem have been considered in the past. Inde-pendent measurement faults were treated, in an LMMSE sense, in [2]. UBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 2
De Koning [7] considered a more general case of independentlyswitching modes where, however, the mode elements are assumeduncorrelated, and Costa [8] developed, by means of state augmen-tation, a recursive LMMSE filter for systems with discrete modesobeying Markov dynamics. Additional contributions include [3], thatconsidered correlated faults, [11], that allowed correlations betweensubsequent fault variables, and [4], that proposed an LMMSE filter forthe static multiple model problem [12]. Related nonlinear solutionswere proposed in [5], [6], [13] and references therein.Besides the novel introduction of the feedback terms, this papercontains several additional contributions. First, we derive a recursiveLMMSE algorithm without assuming uncorrelatedness of the modeelements, as done in [7]. This assumption precludes the utilizationof the algorithm of [7] even for the simple problem of uncertainobservations where measurement noise has a higher variance whenfaults occur, not to mention more involved settings, such as trackingin clutter. In addition, our algorithm is derived without state augmen-tation and without assuming discrete modes, as done in [8]. Finally,the approach allows a broader class of problem to be formulatedwithin a single state-space model. Specifically, the new feedbackterms allow the application of the idea to the problem of trackingin clutter.III. L
INEAR O PTIMAL R ECURSIVE E STIMATION
We begin the derivation with deterministic u k . The stochastic case istreated in Section III-E.Let Y k be the RV obtained by concatenating the elements of Y k . We derive the result using the following lemma, which followsfrom [14, p. 190] and the linearity of the MMSE estimator in theGaussian case. Lemma.
Let x , y and z be RVs and let ˆ x ( z ) and ˆ x ( y, z ) denote,respectively, the LMMSE estimates of x using z , and using both y and z . Let ˆ y ( z ) be the LMMSE estimate of y using z . Then ˆ x ( y, z ) =ˆ x ( z )+Γ x ˜ y Γ − y ˜ y ˜ y, where ˜ y = y − ˆ y ( z ) and Γ ab is the cross-covariancematrix between the RVs a and b . Letting z , Y k , y , y k +1 and using the lemma, the LMMSEestimate of x k +1 using Y k +1 is ˆ x k +1 = ˆ x − k +1 + Γ x k +1 ˜ y k +1 Γ − y k +1 ˜ y k +1 ˜ y k +1 , (3)where ˆ x − k +1 is the LMMSE estimate of x k +1 using Y k , ˜ y k +1 , y k +1 − ˆ y − k +1 , and ˆ y − k +1 is the LMMSE estimate of y k +1 using Y k . If Γ ˜ y k +1 ˜ y k +1 is singular the lemma still holds with the inverse replacedby the Moore-Penrose pseudo-inverse. It is easily verified that ˆ x − k +1 = E [ A k ] ˆ x k + E [ B k ] u k (4) ˆ y − k +1 = E [ H k +1 ] ˆ x − k +1 + E [ F k +1 ] ˆ x k = ( E [ H k +1 ] E [ A k ]+ E [ F k +1 ])ˆ x k + E [ H k +1 ] E [ B k ] u k . (5)Plugging (4) in (3) we identify the desired matrix coefficients K k , L k , and J k of (2) as follows: K k = Γ x k +1 ˜ y k +1 Γ − y k +1 ˜ y k +1 (6) L k = ( I − K k E [ H k +1 ]) E [ A k ] − K k E [ F k +1 ] (7) J k = ( I − K k E [ H k +1 ]) E [ B k ] . (8)We now compute the covariance terms Γ x k +1 ˜ y k +1 and Γ ˜ y k +1 ˜ y k +1 . A. Computation of Γ x k +1 ˜ y k +1 Since ˆ y − k +1 is unbiased, and using (1b) and (5), Γ x k +1 ˜ y k +1 = E h x k +1 ( y k +1 − ˆ y − k +1 ) ⊤ i = E h x k +1 ( H k +1 x k +1 + G k +1 v k +1 + F k +1 ˆ x k ) ⊤ i − E h x k +1 (( E [ H k +1 ] E [ A k ] + E [ F k +1 ])ˆ x k ) ⊤ i − E h x k +1 ( E [ H k +1 ] E [ B k ] u k ) ⊤ i . (9)Using the independence of x k +1 and v k +1 , and canceling outidentical terms, (9) becomes Γ x k +1 ˜ y k +1 = E [ x k +1 x ⊤ k +1 ] E [ H ⊤ k +1 ] − E [ x k +1 ˆ x ⊤ k ] E [ A ⊤ k ] E [ H ⊤ k +1 ] − E [ x k +1 ] u ⊤ k E [ B ⊤ k ] E [ H ⊤ k +1 ] . (10)Before proceeding, we define Σ k , E [ x k x ⊤ k ] , ∆ k , u k u ⊤ k and, inaddition, Λ k , E [ˆ x k ˆ x ⊤ k ] = E [ˆ x k x ⊤ k ] (11) Υ k , E [ x k ] u ⊤ k = E [ˆ x k ] u ⊤ k , (12)where the RHS of (11) and (12) follow from the orthogonalityprinciple and from the unbiasedness of ˆ x k , respectively. Note that Σ k , Λ k , and ∆ k are symmetric.Using the independence of ˆ x k and w k , E h x k +1 ˆ x ⊤ k i = E h ( A k x k + B k u k + C k w k )ˆ x ⊤ k i = E [ A k ] Λ k + E [ B k ] Υ ⊤ k , (13)which yields for (10) Γ x k +1 ˜ y k +1 = (cid:0) Σ k +1 − ( E [ A k ]Λ k + E [ B k ]Υ ⊤ k ) E [ A ⊤ k ] − E [ x k +1 ] u ⊤ k E [ B ⊤ k ] (cid:1) E [ H ⊤ k +1 ] . (14)From (1a) we have E [ x k +1 ] = E [ A k x k + B k u k + C k w k ]= E [ A k ] E [ x k ] + E [ B k ] u k , (15)which, when substituted in (14), leads to Γ x k +1 ˜ y k +1 = (cid:0) Σ k +1 − ( E [ A k ](Λ k E [ A ⊤ k ] + Υ k E [ B ⊤ k ])+ E [ B k ](Υ ⊤ k E [ A ⊤ k ] + ∆ k E [ B ⊤ k ])) (cid:1) E [ H ⊤ k +1 ] . (16) B. Computation of Γ ˜ y k +1 ˜ y k +1 Since ˆ y − k +1 is the LMMSE estimate of y k +1 using Y k , ˜ y k +1 isorthogonal to ˆ y − k +1 and, using (5), Γ ˜ y k +1 ˜ y k +1 = E [( y k +1 − ˆ y − k +1 ) y ⊤ k +1 ]= E [ y k +1 y ⊤ k +1 ] − E [ˆ y − k +1 y ⊤ k +1 ]= E [ y k +1 y ⊤ k +1 ] − ( E [ H k +1 ] E [ A k ] + E [ F k +1 ]) E [ˆ x k y ⊤ k +1 ] − E [ H k +1 ] E [ B k ] u k E [ y ⊤ k +1 ] . (17)Using (1b) and the independence of { ˆ x k , x k +1 } , { H k +1 , G k +1 , F k +1 } and v k +1 , we have E [ˆ x k y ⊤ k +1 ] = E h ˆ x k ( H k +1 x k +1 + F k +1 ˆ x k ) ⊤ i = E [ˆ x k x ⊤ k +1 ] E [ H ⊤ k +1 ] + Λ k E [ F ⊤ k +1 ] , (18)which, using (13), becomes E [ˆ x k y ⊤ k +1 ] = Λ k ( E [ A ⊤ k ] E [ H ⊤ k +1 ] + E [ F ⊤ k +1 ])+ Υ k E [ B ⊤ k ] E [ H ⊤ k +1 ] . (19) UBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 3
Due to the independence of { x k +1 , ˆ x k } , v k +1 , and { H k +1 , G k +1 } E [ y k +1 y ⊤ k +1 ] = E [ H k +1 x k +1 x ⊤ k +1 H ⊤ k +1 ] + E [ G k +1 v k +1 v ⊤ k +1 G ⊤ k +1 ]+ E [ F k +1 ˆ x k ˆ x ⊤ k F ⊤ k +1 ] + E [ H k +1 x k +1 ˆ x ⊤ k F ⊤ k +1 ]+ E [ F k +1 ˆ x k x ⊤ k +1 H ⊤ k +1 ] . (20)Consider the last summand. From the smoothing property of theconditional expectation, E [ F k +1 ˆ x k x ⊤ k +1 H ⊤ k +1 ] = E h E [ F k +1 ˆ x k x ⊤ k +1 H ⊤ k +1 | F k +1 , H k +1 ] i = E h F k +1 E [ˆ x k x ⊤ k +1 ] H ⊤ k +1 i , (21)where we utilized the independence of { H k +1 , F k +1 } and { x k +1 , ˆ x k } .Similarly, since E (cid:2) x k +1 x ⊤ k +1 (cid:3) = Σ k +1 , E (cid:2) v k +1 v ⊤ k +1 (cid:3) = I , and E (cid:2) ˆ x k ˆ x ⊤ k (cid:3) = Λ k , we obtain: E [ H k +1 x k +1 x ⊤ k +1 H ⊤ k +1 ] = E [ H k +1 Σ k +1 H ⊤ k +1 ] (22) E [ G k +1 v k +1 v ⊤ k +1 G ⊤ k +1 ] = E [ G k +1 G ⊤ k +1 ] (23) E [ F k +1 ˆ x k ˆ x ⊤ k F ⊤ k +1 ] = E [ F k +1 Λ k F ⊤ k +1 ] . (24)For future reference, we also note that E [ A k x k x ⊤ k A ⊤ k ] = E [ A k Σ k A ⊤ k ] (25) E [ A k x k u ⊤ k B ⊤ k ] = E [ A k Υ k B ⊤ k ] (26) E [ B k u k u ⊤ k B ⊤ k ] = E [ B k ∆ k B ⊤ k ] (27) E [ C k w k w ⊤ k C ⊤ k ] = E [ C k C ⊤ k ] . (28)Substituting (13) in (21), and using (21)-(24) in (20), E [ y k +1 y ⊤ k +1 ] = E [ H k +1 Σ k +1 H ⊤ k +1 ] + E [ G k +1 G ⊤ k +1 ]+ E [ F k +1 Λ k F ⊤ k +1 ]+ E h H k +1 ( E [ A k ] Λ k + E [ B k ] Υ ⊤ k ) F ⊤ k +1 i + E h F k +1 (Λ k E [ A ⊤ k ] + Υ k E [ B ⊤ k ]) H ⊤ k +1 i . (29)In addition, we obtain, in a straightforward manner, E [ y k +1 ] = ( E [ H k +1 ] E [ A k ] + E [ F k +1 ]) E [ x k ]+ E [ H k +1 ] E [ B k ] u k . (30)Using (22), (23), and (24) in (29), and substituting (19), (29), and (30)in (17) we finally obtain Γ ˜ y k +1 ˜ y k +1 = E h H k +1 Σ k +1 H ⊤ k +1 i + E h G k +1 G ⊤ k +1 i + E h F k +1 Λ k F ⊤ k +1 i − E [ F k +1 ] Λ k E [ F k +1 ] ⊤ − E [ H k +1 ] E [ A k ] Λ k E h A ⊤ k i E h H ⊤ k +1 i + E h H k +1 ( E [ A k ] Λ k + E [ B k ] Υ ⊤ k ) F ⊤ k +1 i + E h F k +1 (Λ k E h A ⊤ k i + Υ k E h B ⊤ k i ) H ⊤ k +1 i − E [ H k +1 ] E [ A k ] Λ k E h F ⊤ k +1 i − E [ F k +1 ] Λ k E h A ⊤ k i E h H ⊤ k +1 i − E [ H k +1 ] E [ A k ] Υ k E h B ⊤ k i E h H ⊤ k +1 i − E [ F k +1 ] Υ k E h B ⊤ k i E h H ⊤ k +1 i − E [ H k +1 ] E [ B k ] u k E h y ⊤ k +1 i . (31)Notice, that a sufficient condition for the nonsingularity of Γ ˜ y k +1 ˜ y k +1 is E [ G k +1 G ⊤ k +1 ] ≻ . To see this, recall that, by definition, Γ ˜ y k +1 ˜ y k +1 is positive semi-definite for any choice of E (cid:2) G k +1 G Tk +1 (cid:3) and, in particular, for G k +1 = 0 . But this means that the matrix on theRHS of (31) without E (cid:2) G k G Tk (cid:3) is positive semi-definite, rendering E (cid:2) G k G Tk (cid:3) ≻ a sufficient condition for the non-singularity of Γ ˜ y k +1 ˜ y k +1 . C. Computation of the Second-Order Moments
Utilizing the independence of x k , w k and { A k , B k , C k } , and (25)–(28), Σ k +1 is given by Σ k +1 = E [ x k +1 x ⊤ k +1 ]= E h ( A k x k + B k u k + C k w k )( A k x k + B k u k + C k w k ) ⊤ i = E [ A k Σ k A ⊤ k ] + E [ A k Υ k B ⊤ k ] + E [ B k Υ ⊤ k A ⊤ k ]+ E [ B k ∆ k B ⊤ k ] + E [ C k C ⊤ k ] , (32)Next consider Λ k +1 . Direct computation yields: Λ k +1 = E h ˆ x k +1 x ⊤ k +1 i = ( L k + K k E [ F k +1 ]) E [ˆ x k x ⊤ k +1 ]+ K k E [ H k +1 ] Σ k +1 + J k u k E [ x ⊤ k +1 ] . (33)Using (13) the latter becomes Λ k +1 = ( L k + K k E [ F k +1 ])(Λ k E [ A ⊤ k ] + Υ k E [ B ⊤ k ])+ J k (Υ ⊤ k E [ A ⊤ k ] + ∆ k E [ B ⊤ k ]) + K k E [ H k +1 ]Σ k +1 . (34)Finally, Υ k +1 = E [ x k +1 ] u ⊤ k +1 . Note that ∆ k is known for all k . D. Algorithm Summary a) Initialization: ˆ x = ¯ x , Σ = P + ¯ x ¯ x ⊤ , Λ = ¯ x ¯ x ⊤ , Υ =¯ x u ⊤ , ∆ = u u ⊤ .b) Recursion: For k = 1 , , . . . perform the routine of Alg. 1. Algorithm 1Input: y k +1 , u k +1 , ˆ x k , E [ x k ] , Σ k , Λ k , Υ k , ∆ k Compute E [ A k ] , E [ B k ] , E (cid:2) C k C ⊤ k (cid:3) , E (cid:2) A k Σ k A ⊤ k (cid:3) , E (cid:2) A k Υ k B ⊤ k (cid:3) , and E (cid:2) B k ∆ k B ⊤ k (cid:3) . Compute E [ x k +1 ] and Σ k +1 using Eqs. (15) and (32). Compute E [ H k +1 ] , E (cid:2) G k +1 G ⊤ k +1 (cid:3) , E [ F k +1 ] , E (cid:2) F k +1 Λ k F ⊤ k +1 (cid:3) , E (cid:2) H k +1 Σ k +1 H ⊤ k +1 (cid:3) , and E (cid:2) H k +1 ( E [ A k ] Λ k + E [ B k ] Υ ⊤ k ) F ⊤ k +1 (cid:3) . Compute Γ x k +1 ˜ y k +1 and Γ ˜ y k +1 ˜ y k +1 using Eqs. (16) and (31). Compute K k , L k , and J k using Eqs. (6), (7), and (8), and ˆ x k +1 using Eq. (2). Compute Λ k +1 using Eq. (34) and Υ k +1 by plugging E [ x k +1 ] into (12). Output: ˆ x k +1 , E [ x k +1 ] , Σ k +1 , Λ k +1 , Υ k +1 Since the distribution of M k is known, the expectations of steps 1and 3 of Alg. 1 may be calculated by, e.g., direct summations in caseof discrete modes. In some cases, as demonstrated in Section IV,closed form expressions exist for the above expectations.We note that the standard KF for a system with no inputs shouldbe obtained when {M k } is a deterministic sequence with B k = 0 , F k = 0 . In this setting we have Γ x k +1 ˜ y k +1 = (Σ k +1 − A k Λ k A ⊤ k ) H ⊤ k +1 and Γ ˜ y k +1 ˜ y k +1 = H k +1 (Σ k +1 − A k Λ k A ⊤ k ) H ⊤ k +1 + G k +1 G ⊤ k +1 . Substituting these in (3) we indeed obtain the standard KF in the formwhere the time and measurement updates are combined together. Theerror covariances follow in a similar manner.
UBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 4
E. Random Inputs
In the second variant of (1a), in which u k = ˆ x k , it turns out that theroles played by A k and B k are identical. Specifically, after replacing u k with ˆ x k , at each step of the derivation of Section III, A k and B k are multiplied by the same quantities. Thus, the filter for the modifiedproblem is obtained from the one described in Alg. 1 by replacing A k with A k + B k and nullifying u k and Υ k . An alternative derivation,based on the orthogonality principle, may be found in [15].IV. A PPLICATION TO T ARGET T RACKING IN C LUTTER
In this section we demonstrate the proposed concept by casting theclassical problem of tracking in clutter within our formulation, andapplying the LMMSE filter of Section III.
A. System and Clutter Models
Consider a single target obeying a linear model. Setting A k = A , B k = 0 , and C k = C in (1a) x k +1 = Ax k + Cw k . (35)Here A and C are deterministic matrices, accounting for the statedynamics and process noise covariance, respectively, and { w k } is ascalar process noise sequence. The target state is observed via thethe equation y true k = H nom x k + G nom v true k , (36)where v true k represents measurement noise. In addition, at each time,a number of clutter detections are obtained. These will be denotedas { y cl k,i } N − i =1 , where N is the total number of detections. Cluttermeasurements do not carry any information about the target ofinterest. They are, however, indistinguishable from true detections inthe sense that they carry information of the same type (say, position).At each time, the clutter measurements are assumed to be independentof each other, of the clutter measurements at other times, and ofthe true state and observation. In addition, we assume that they areuniformly distributed in space. To correctly model the distributionof the clutter detections, we note that, typically, at each scan, thesensor initiates a validation window centered about the predictedtarget position, and the algorithm processes only those measurementsobtained within the window. Since the clutter detections are uniformlydistributed in space, they are also uniformly distributed within thevalidation window.We define the measurement vector y k to be the concatenation ofall measurements from time k , N − of which correspond to clutter,and one originating from the true target. The location of the truemeasurement within this concatenated vector is, of course, unknownto the algorithm. This setting can be modeled using (1b) by lettingthe mode M k be distributed as M k = { H k , G k , F k } = H nom ... ! , diag G nom G cl ...G cl , H nom A...H nom A ! , w.p. N ... ... ... ... H nom ! , diag G cl ...G cl G nom , H nom A...H nom A ! , w.p. N , (37)where G cl is the square-root of the covariance matrix associated withthe clutter.For example, the first realization of { H k , G k , F k } in (37) corre-sponds to the scenario in which the first of the N observations is the true target measurement, y true k , generated according to (36), while theother N − measurements are clutter, each of which is generatedaccording to y cl k,i = H nom A ˆ x k − + G cl v cl k,i , i = 2 , . . . , N. (38)Here, H nom A ˆ x k − is the predicted true measurement at time k ,which is also the center of the validation window, so that cluttermeasurements at time k are uniformly distributed around this quan-tity. Namely, v cl k,i has a uniform distribution. The overall number ofmeasurements in the validation window, N , is assumed to be known,but may vary in time. Thus, the dimensions of H k , G k , and F k maydepend on k .It is readily observed that the matrices { H k , G k , F k } are correlatedin this setting. This renders the approach of [7] inapplicable inthe current scenario. Furthermore, it can be seen that without thefeedback term in the measurement equation, it is impossible toaccount for the fact that clutter is uniformly distributed in a windowcentered about the predicted measurement. In fact, any linear methoddisregarding this term, such as [7], [8], must assume that cluttermeasurements are distributed about .Notice that we assumed, for simplicity, that the true measurement isalways present in the validation window. To account for the possibilitythat the true measurement does not fall in the validation window, theoption { H k , G k , F k } = { , I N ⊗ G cl , N ⊗ H nom A } needs to be added to the set of possible realizations in (37). Here, ⊗ stands for the Kronecker product, N is an N × vector comprisingall ones, and I N is the N × N identity matrix. The probability ofthis outcome is (1 − P D )(1 − P G ) where P D is the probabilty oftarget detection, assumed known, and P G is the probability that,upon target detection, the true measurement falls in the validationwindow. This parameter is defined by the user and, typically, it affectsthe window size as discussed in the sequel. Note that, when nomeasurements are available, N = 0 , and (2) becomes (at the absenceof u k ) ˆ x k +1 = L k ˆ x k , which corresponds to a simple prediction (timeupdate) without consecutive measurement update, as expected. B. Matrix Computations
To invoke the algorithm presented in Section III we need tocompute the expectations of Steps 1 and 3 of Alg. 1. Althoughthese may be evaluated numerically, via direct summations, in thepresent example closed-form expressions exist, as we show nextfor the simple setting in which the true measurement is alwayspresent in the validation window (extensions are straightforward.)As the matrices of the dynamics equation are deterministic, E [ A k ] = A , E [ B k ] = 0 , E (cid:2) C k C ⊤ k (cid:3) = CC ⊤ , E (cid:2) A k Υ k B ⊤ k (cid:3) = 0 , E (cid:2) B k ∆ k B ⊤ k (cid:3) = 0 , and E (cid:2) A k Σ k A ⊤ k (cid:3) = A Σ k A ⊤ . Also, accordingto the distribution defined in (37), E [ H k +1 ] = 1 N N ⊗ H nom (39) E [ F k +1 ] = N − N N ⊗ H nom A. (40)The remaining terms read E [ H k +1 Σ k +1 H ⊤ k +1 ] = 1 N I N ⊗ H nom Σ k +1 H ⊤ nom (41) E [ G k +1 G ⊤ k +1 ] = 1 N I N ⊗ (cid:16) G nom G ⊤ nom + ( N − G cl G ⊤ cl (cid:17) (42) E h F k +1 Λ k F ⊤ k +1 i = Ξ ⊗ (cid:16) H nom A Λ k A ⊤ H ⊤ nom (cid:17) , (43) UBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 5 where
Ξ = ( N (cid:0) ( N − N ⊤ N + I N (cid:1) , N > , N = 1 . (44)Finally, E h H k +1 ( E [ A k ] Λ k + E [ B k ] Υ ⊤ k ) F ⊤ k +1 i = 1 N ( N ⊤ N − I N ) ⊗ (cid:16) H nom A Λ k A ⊤ H ⊤ nom (cid:17) . (45)The spatial distribution of clutter is uniform in the validationwindow, whose size determines G cl G ⊤ cl . C. Discussion
It is easy to see that, in the present case, Γ ˜ y k +1 ˜ y k +1 = I N ⊗ D where D = 1 N H nom A Λ k A ⊤ H ⊤ nom + 1 N H nom Σ k +1 H ⊤ nom + 1 N G nom G ⊤ nom + N − N G cl G ⊤ cl . Moreover, Γ x k +1 ˜ y k +1 = (Σ k +1 − A Λ k A ⊤ ) E h H ⊤ k +1 i = 1 N (Σ k +1 − A Λ k A ⊤ ) (cid:16) H ⊤ nom · · · H ⊤ nom (cid:17) ⊤ , (46)and K k = Γ x k +1 ˜ y k +1 Γ − y k +1 ˜ y k +1 = 1 N ⊤ N ⊗ (cid:16) (Σ k +1 − A Λ k A ⊤ ) H ⊤ nom D − (cid:17) . (47)Since y k +1 is a concatenation of all the observations from time k +1 ,the product K k y k +1 in (2) is the average of these measurements,pre-multiplied by (Σ k +1 − A Λ k A ⊤ ) H ⊤ nom D − . Consequently, theLMMSE estimator for tracking a target in clutter is a KF-likealgorithm, operating on the average of all detections in the validationwindow. In this respect, its mode of operation resembles classicalmethods. For example, the probabilistic data association (PDA) [16]method implements a KF driven by the weighted average of all mea-surements in the window, and the nearest neighbor (NN) filter [17] isa KF driven by the measurement nearest to the prediction assigningit a weight of and assigning to the rest of the measurements. D. Numerical Study
We consider a one-dimensional tracking scenario, in which thestate comprises position and velocity information, x k = ( p k v k ) ⊤ .Starting at x ∼ N (¯ x , P ) with ¯ x = (0 0) ⊤ and P = 30 I , thetarget is simulated for time units using (35) with A = ( .
20 0 . ) and C = (cid:0) / (cid:1) . The process and measurement noises are takento be Gaussian. The true measurement is generated using (36) with H nom = (1 0) and G nom = √ . The target is detected withprobability P D = 0 . and the probability that the true observationfalls in the validation window is taken to be P G = 0 . . A validationwindow is set about the predicted measurement position. Its size, d ,is determined to comply with P G (see [17, p.130] for details). Oncethe window is determined, the clutter variance of (38) is G cl G ⊤ cl = d / .The derived algorithm is compared with NN and PDA filters,that are equipped with the same windowing logic and parameters.All algorithms are initialized with ˆ x = ¯ x and the initial errorcovariance matrix is taken to be P . When dealing with trackingin clutter, using the MSE as the only performance measure mayresult in misleading conclusions, since, eventually the estimate will draw away from the true measurement and follow the clutter, and theerrors will become meaninglessly large. We thus use two measuresof performance to evaluate the algorithms. The first is the time untilthe target is lost, defined as the third consecutive time when themeasurement of a detected target falls outside the validation window.The second measure is the root MSE (RMSE) calculated over thetime interval until the first of the three algorithms loses track.We test the algorithms at a range of clutter densities. Let ρ to bethe average number of clutter measurements falling in an intervalof one standard deviation of the (true) measurement noise. Averagedover independent Monte Carlo runs, the average position RMSEand track loss times are plotted, versus ρ , in Fig. 1. It is readily seenthat the LMMSE filter attains competitive performance relatively tothe nonlinear algorithms. Specifically, for heavy clutter regimes itmaintains longest track loss times. It is not very surprising that theerrors of PDA are better, since these are calculated before the firstof the three algorithms has lost track (NN in all cases). During thisperiod the PDA performs a more efficient, nonlinear manipulation onthe measurements. However, for high clutter rates, it is probable thatclutter measurements will be assigned higher weights than the truedetection, eventually leading to a track loss. In this case, it is betterto simply average the measurements, as the linear filter does.V. C ONCLUSION
We proposed a new formulation of JLS, where the dynamics andmeasurement equations are allowed to depend on previous estimatesof the state representing closed-loop control input and measurementvalidation window. We derived an LMMSE recursive algorithm forthis setting, and illustrated the approach in the context of tracking inclutter. In this case, our filter demonstrates competitive performance,when compared with classical, nonlinear methods.R
EFERENCES [1] E.-K. Boukas and Z.-K. Liu,
Deterministic and Stochastic Time-DelaySystems . Boston: Birkh¨auser, 2002.[2] N. Nahi, “Optimal recursive estimation with uncertain observation,”
IEEE Trans. Inf. Theory , vol. IT-15, no. 4, pp. 457–462, 1969.[3] M. Hadidi and S. Schwartz, “Linear recursive state estimators underuncertain observations,”
IEEE Trans. Autom. Control , vol. AC-24, no. 6,pp. 944–948, December 1979.[4] N. Nahi and E. Knobbe, “Optimal linear recursive estimation withuncertain system parameters,”
IEEE Trans. Autom. Control , vol. 21,no. 2, pp. 263–266, 1976.[5] G. Ackerson and K. Fu, “On state estimation in switching environments,”
IEEE Trans. Autom. Control , vol. 15, no. 1, pp. 10–17, 1970.[6] H. Blom and Y. Bar-Shalom, “The interacting multiple model algorithmfor systems with Markovian switching coefficients,”
IEEE Trans. Autom.Control , vol. 33, no. 8, pp. 780–783, 1988.[7] W. De Koning, “Optimal estimation of linear discrete-time systems withstochastic parameters,”
Automatica , vol. 20, no. 1, pp. 113–115, 1984.[8] O. Costa, “Linear minimum mean square error estimation for discrete-time Markovian jump linear systems,”
IEEE Trans. Autom. Control ,vol. 39, no. 8, pp. 1685–1689, 1994.[9] T. Michaeli, D. Sigalov, and Y. Eldar, “Partially linear estimation withapplication to sparse signal recovery from measurement pairs,”
IEEETrans. Signal Process. , vol. 60, no. 5, pp. 2125–2137, 2012.[10] T. Yuan, Y. Bar-Shalom, P. Willett, E. Mozeson, S. Pollak, and D. Hardi-man, “A multiple IMM estimation approach with unbiased mixing forthrusting projectiles,”
IEEE Trans. Aerosp. Electron. Syst. , vol. 48, no. 4,pp. 3250–3267, 2012.[11] R. Jackson and D. Murthy, “Optimal linear estimation with uncertainobservations (Corresp.),”
IEEE Trans. Inf. Theory , vol. 22, no. 3, pp.376–378, 1976.[12] D. Magill, “Optimal adaptive estimation of sampled stochastic pro-cesses,”
IEEE Trans. Autom. Control , vol. 10, no. 4, pp. 434–439, 1965.[13] D. Sigalov and Y. Oshman, “State estimation in hybrid systems witha bounded number of mode transitions,” in
Proc. Fusion 2010, 13thInternational Conference on Information Fusion , 2010.
UBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 6 ρ P o s i t i o n R M S E LMMSEPDANN ρ T r a c k L o ss T i m e s LMMSEPDANN
Fig. 1: Position RMSE (left) and track loss time (right) vs. clutter density. [14] J. Mendel,
Lessons in digital estimation theory . Prentice-Hall, Inc.,1986.[15] D. Sigalov, T. Michaeli, and Y. Oshman, “Linear optimal state estimationin systems with independent mode transitions,” in
Proc. CDC 2011, 50thConf. on Decision and Control . IEEE, 2011.[16] Y. Bar-Shalom and E. Tse, “Tracking in a cluttered environment withprobabilistic data association,”
Automatica , vol. 11, no. 5, pp. 451–460,1975.[17] Y. Bar-Shalom and X. Li,