1

Bayesian Filtering with Unknown Sensor
Measurement Losses

arXiv:1801.07945v1 [] 24 Jan 2018

Jiaqi Zhang, Keyou You, Senior Member, IEEE, and Lihua Xie, Fellow, IEEE

Abstract—This work studies the state estimation problem of a
stochastic nonlinear system with unknown sensor measurement
losses. If the estimator knows the sensor measurement losses of a
linear Gaussian system, the minimum variance estimate is easily
computed by the celebrated intermittent Kalman filter (IKF).
However, this will no longer be the case when the measurement
losses are unknown and/or the system is nonlinear or nonGaussian. By exploiting the binary property of the measurement
loss process and the IKF, we design three suboptimal filters for
the state estimation, i.e., BKF-I, BKF-II and RBPF. The BKF-I
is based on the MAP estimator of the measurement loss process
and the BKF-II is derived by estimating the conditional loss
probability. The RBPF is a particle filter based algorithm which
marginalizes out the loss process to increase the efficiency of
particles. All the proposed filters can be easily implemented in
recursive forms. Finally, a linear system, a target tracking system
and a quadrotor’s path control problem are included to illustrate
their effectiveness, and show the tradeoff between computational
complexity and estimation accuracy of the proposed filters.
Index Terms—Stochastic systems, networked estimation, intermittent Kalman filter, sensor measurement losses, particle filter.

I. I NTRODUCTION

T

HE Kalman filter (KF) has a simple structure in optimally correcting propagated state estimates with the
sensor measurement of the observed system, and has been
successfully applied in countless guidance, navigation and
control (GNC) related applications. While in some practical
applications, the estimator may not always access the true
sensor measurement [1]–[9]. For example, a sensor fault
results in that the estimator only receives a pure noise [4]–
[6], which does not contain any information of the estimated
system. In networked systems, the channel traffic congestion
will also lead to sensor measurement losses [7]. Similar issues
commonly exist in the following important applications.
• In a target tracking system, the position, distance, and
velocity of the target to the base station is measured
by a radar or Global Positioning System (GPS). The
sensor information is sent to a missile or an unmanned
aircraft via resource limited wireless channels where
communications between devices are power constrained
and therefore limited in range and reliability [7]–[9]. If
the true sensor measurement is lost, the base station only
receives error data, e.g., channel noises.
J. Zhang and K. You are with the Department of Automation, and TNList, Tsinghua University, Beijing 100084, China. e-mail:
zjq16@mails.tsinghua.edu.cn, youky@tsinghua.edu.cn.
L. Xie is with the School of Electrical and Electronic Engineering,
Nanyang Technological University, Singapore 639798, Singapore. e-mail:
elhxie@ntu.edu.sg.

GPS is a common way for positioning, which however
requires signals to be detected precisely. Thus any slight
inaccuracy in the signal’s reception or some external
disturbances can lead to an error in the measured location [2]. Environmental factors such as trees, valleys,
buildings, or even heavy cloud cover can impact or even
interrupt the transmission of the GPS signal between the
receiver and the satellites. Besides, urban environments
surrounded by tall buildings also cause inaccurate GPS
signals.
• A vehicle or aircraft is usually equipped with many
sensors. For instance, a quadrotor may have an inertial
measurement unit (IMU), a GPS receiver, barometers,
and magnetometers. Large aircrafts even use more sensor
units. It is possible that some error measurements occur
from time to time due to unexpected faults in the system
[3]–[5].
In the above situations, the KF or extended KF can only
intermittently access the true sensor measurements. If an
wrong measurement or pure noise is used for updating the
state estimate, the estimation performance will significantly
degrade, and the filter may even diverge.
Then, a natural question is how to handle issues induced by
occasional measurement errors or losses. If the measurement
loss is known to the estimator, the minimum variance estimate
(MVE) of a linear Gaussian system is provided by an intermittent Kalman filter (IKF) [10], and a vast body of literature, see
e.g. [11], [12] and references therein, focuses on the stability
of the IKF. For instance, it is shown in [10] that there exists
a critical measurement loss probability beyond which the IKF
may diverge. In [13], this critical value for certain type of
systems is explicitly obtained under Markovian measurement
losses. More cases can be found in [14].
The problem is further complicated if the measurement
loss is unknown to the estimator. In real applications, the
estimator receives a data packet which may be a pure noise or
fake sensor measurement. A faulty sensor might also return
“wrong” measurement data to the estimator. Thus, it is of
practical importance to study the state estimation problem
under unknown sensor measurement losses, which is the
main focus of this work. Clearly, the unknown measurement
loss results in a non-Gaussian and nonlinear system where
the IKF is no longer applicable. Although there are some
generic estimation methods for nonlinear and/or non-Gaussian
systems, e.g., extended KF (EKF), unscented KF [15], [16]
and particle filter (PF) [17], they do not particularly explore
the unique feature of the current problem, in which the measurement loss or error results in a discontinuous measurement
•

2

equation. To resolve this issue, we model the process of sensor
measurement losses as a binary sequence {γk }. Particularly,
γk = 1 means that the true sensor measurement is received
and γk = 0 indicates the occurrence of a sensor measurement
loss. Our objective is to design practical suboptimal filters to
accommodate the unknown intermittent measurement losses in
a nonlinear stochastic system.
Motivated by the optimality of the IKF, our idea is that we
first estimate the binary sequence {γk } under the maximum
a posteriori probability (MAP) estimation criterion. Following
the principle of certainty equivalence, we then use the MAP
estimate to replace the unknown measurement loss γk in the
IKF for the linearized systems. This results in our Bayesian
Kalman filter I (BKF-I). We also derive a Bayesian Kalman
filter II (BKF-II) by estimating the conditional measurement
loss probability in the linearized system, which is a compromise between the KF and the IKF. Clearly, both filters are very
different from the existing nonlinear extended KF, unscented
KF or PF, and will reduce to the intermittent extended Kalman
filter (IEKF) if the measurement loss process {γk } is known
to the estimator.
Another method to address the non-Gaussianity and/or nonlinearity is the PF which samples particles to approximate
the conditional density of state. However, the amount of
computations required for the PF in the high-dimensional state
space is usually large. To increase the sampling efficiency, one
can marginalize out some of the states and use the standard
algorithms such as the KF to estimate them. Then, the PF is
applied to estimate the rest of state variables, which is called
Rao-Blackwellised particle filter (RBPF). The implementation
and comparison between the standard PF and RBPF are well
documented in [15], [18]–[20]. In this work, we adopt this
idea by using the PF to estimate the conditional distribution of
the measurement loss process {γk }, which is a binary process
and requires only a few number of particles to approximate
its posterior distribution. As in the BKF-I, the state of the
linearized system is then estimated by the IKF. Compared to
the standard PF, the efficiency of using particles is significantly
improved, no matter how large the dimension of the system
state is.
We compare the computational complexities and the estimation performance of the proposed suboptimal filters. Both the
BKF-I and BKF-II have similar computational complexities
to the IEKF. The computational complexity of the RBPF
is essentially proportional to the number of particles. The
larger number of particles is adopted in the RBPF, the better
estimation accuracy is expected in general. Moreover, the
RBPF weakly converges to an optimal MVE for a linear
stochastic system with unknown measurement losses if the
number of particles tends to infinity. Although the estimation
accuracies of the BKF-I and BKF-II cannot be ensured, they
are usually easier to implement. Thus, the choice of the three
filters depends on the problem on hand.
Finally, we use the proposed filters to solve the state
estimation problem of a linear stochastic system, a target
tracking problem and a quadrotor’s path control problem,
all of which are subject to unknown measurement losses.
In these examples, the proposed filters well complete their

estimation tasks under different levels of measurement losses,
and their performance even comes close to the case with
known measurement losses, showing that the measurement
loss process {γk } is correctly estimated. For the RBPF, we
also illustrate how the estimation accuracy is improved by
increasing the number of particles in the linear system. This
explains the tradeoff between computational complexity and
the estimation accuracy of the RBPF.
It should be noted that the conference version of this paper
[21] only studies linear stochastic systems. In comparison,
this work focuses on nonlinear systems and provides detailed
comparisons of the proposed three filters in terms of estimation
accuracy and computational complexity, see Section V and
Section VI-A. Besides, this paper includes the quadrotor’s path
control problem in Section VI-C to validate the effectiveness
of our filters and provides a fast algorithm for implementing
the RBPF in Section IV-C.
The rest of this paper is organized as follows. In Section II,
we formulate the estimation problem with unknown measurement losses. In Section III, we design the BKF-I and the BKFII based on the Bayes’ theorem and the technique of EKF. In
Section IV, we derive the RBPF to address the estimate of
the process of unknown measurement losses. We compare the
three filters both in terms of computational complexity and
estimation accuracy in Section V. Simulation is performed in
Section VI to show the effectiveness of the proposed filters and
compare their performance. Finally, we draw some concluding
remarks in Section VII.
II. P ROBLEM FORMULATION
In this section, we formulate the state estimation problem
with unknown measurement losses and introduce the celebrated intermittent KF, whose optimality is ensured for a linear
Gaussian system with known measurement losses.
A. State Estimation with Unknown Measurement Losses
Consider a discrete dynamical system with measurement
losses and the additive noise of the form:
xk+1 = f (xk , uk ) + wk
yk = ~γk · h(xk ) + vk

(1)

where xk ∈ Rn and yk ∈ Rm are the vector state and
measurement, respectively. The random vectors wk ∈ Rn and
vk ∈ Rm are independent white Gaussian noises with zero
means and covariance matrices Q ≥ 0 and R ≥ 0, respectively.
The initial state x0 is assumed to be a random Gaussian vector
with mean x̄0 and covariance matrix Σ0 > 0. ~γk is a diagonal
matrix with its i-th diagonal element γik , i ∈ {1, 2, .., m},
representing the i-th sensor measurement loss at time step k. In
particular, γik = 1 indicates that the true sensor measurement
yik is successfully received by the estimator while γik = 0
means that the estimator only receives a fake measurement of
yik , i.e., the pure noise vik .
In many real systems, all entries of yk are obtained from
the same sensor and are transmitted as a single packet, e.g.,
the GPS signal usually consists of position and velocity
measurements. Then, all γik , i ∈ {1, ..., m} are identical, and

3

~γk is reduced to a random variable. For brevity, this paper
mainly studies this case and focuses on the following nonlinear
stochastic system
xk+1 = f (xk , uk ) + wk

(2)

yk = γk · h(xk ) + vk

where γk is a binary random variable that represents the sensor
measurement loss at time step k. Thus, γk = 1 indicates that
the true sensor measurement yk is contained in the arrival data
packet while γk = 0 means that the estimator only receives
pure noise.
The goal of this work is to propose recursive filters to
estimate the state of the system in (2) under unknown measurement loss process {γk }.
B. Intermittent Extended KF (IEKF)
For a linear Gaussian model, the IKF in [10] assumes that
{γk } is known to the estimator. To revise it for the nonlinear
system in (2), define Γk = {γ0 , ..., γk }, Uk = {u0 , ..., uk }
and Yk = {y0 , ..., yk }, and the conditional minimum variance
estimate and error covariance matrices are given by
x̂k|k−1 = E[xk |Yk−1 , Γk−1 , Uk−1 ]

A. Bayesian Kalman Filter I
We design a nonlinear filter called Bayesian Kalman filter
I (BKF-I) to recursively compute the state estimate with
unknown measurement losses. An intuitive idea, which is
motivated by the principle of certainty equivalence, is that we
first estimate the measurement losses Γk , based on which the
IEKF (4) is then applied to compute the state estimate. We
shall elaborate it in this subsection.
Since γk is binary, it is natural to adopt the maximum a
posteriori probability (MAP) estimate, and the MAP estimate
of Γk is given as follows:
b k = argmax p(Γk |Yk )
Γ
Γk
b k = {γ̂0 , ..., γ̂k } and for notional simplicity, we
where Γ
directly use p(x) to denote either the probability density or
mass function of a random vector X.
b k into (4), i.e., use the estimate γ̂k to replace
Substitute Γ
unknown measurement loss γk in (4), we obtain the BKF-I, see
Algorithm 1. The remaining problem reduces to the derivation
of the MAP estimate of Γk . To solve it, we use the Bayes’
formulas and obtain that
p(Γk |Yk )

x̂k|k = E[xk |Yk , Γk , Uk−1 ]
T

Σk|k−1 = E[(xk − x̂k|k−1 )(xk − x̂k|k−1 ) |Yk−1 , Γk−1 , Uk−1 ]

=

Σk|k = E[(xk − x̂k|k )(xk − x̂k|k )T |Yk , Γk , Uk−1 ]
,u )

p(γk , Γk−1 |Yk )
p(γk , Γk−1 , yk |Yk−1 )p(Yk−1 )
.
p(Yk )

To recursively compute the above, we note that

ŷk|k−1 = E[yk |Yk−1 , Γk−1 , Uk−1 ].
∂f (x̂

=

∂h(x̂

)

k
k|k−1
k|k
and Ck =
, this enables
Denote Ak =
∂x
∂x
us to approximate (2) by a linearized system

xk+1 = Ak xk + wk + bk
yk = γk · Ck xk + vk + zk

p(γk , Γk−1 , yk |Yk−1 )
= p(yk |γk , Γk−1 , Yk−1 )p(γk , Γk−1 |Yk−1 )
= p(yk |γk , Γk−1 , Yk−1 )p(γk |Γk−1 , Yk−1 )p(Γk−1 |Yk−1 ).

(3)
This implies that

where bk and zk are computed online from the equations
bk = f (x̂k|k , uk ) − Ak x̂k|k , zk = γk (h(x̂k|k−1 ) − Ck x̂k|k−1 ).
If γk is known to the estimator, the measurement update
for the nonlinear system (2) can be obtained by applying
the IKF [10] to the linearized system (3), which leads to the
Intermittent Extended KF (IEKF), i.e.,
x̂k|k = x̂k|k−1 + γk Kk (yk − h(x̂k|k−1 ))
Σk|k = Σk|k−1 − γk Kk Ck Σk|k−1

(4)

and the time update is the same as the EKF, i.e.,
x̂k+1|k = f (x̂k|k , uk ),
Σk+1|k = Ak Σk|k AT
k + Q,

(5)

where the Kalman gain Kk = Σk|k−1 CkT (Ck Σk|k−1 CkT +
R)−1 and x̂0|−1 = x̄0 , Σ0|−1 = Σ0 .
In the present situation, {γk } is unknown to the estimator,
which renders the above IEKF inapplicable. However, it is still
very helpful in designing an effective filter in a recursive form.
III. BAYESIAN K ALMAN F ILTERS
In this section, we propose two suboptimal recursive filters
to solve the filtering problem with unknown sensor measurement losses.

p(Γk |Yk ) =

p(yk |γk , Γk−1 , Yk−1 )p(γk |Γk−1 , Yk−1 )
p(yk |Yk−1 )
× p(Γk−1 |Yk−1 ).

(6)

To compute p(Γk |Yk ), it requires to consider all possible
values of Γk , which grows unboundedly with respect to the
time steps. Therefore, it is impossible to recursively obtain the
MAP estimate of Γk . In real applications, recursive algorithms
are essential. Thus, we consider to approximately compute
p(Γk |Yk ) in a recursive way. This is achieved by using the
b k−1 rather than the unknown Γk−1 to estimate γk .
estimate Γ
b k−1 , it follows from (6) that
By substituting Γk−1 with Γ
b k−1 |Yk )
p(Γk |Yk ) ≈ p(γk , Γ
(7)
b
b
b
∝ p(yk |γk , Γk−1 , Yk−1 )p(γk |Γk−1 , Yk−1 )p(Γk−1 |Yk−1 ).
b k−1 is already obtained at time step k, our objective
Since Γ
is to find γk to maximize the approximated posterior probab k−1 |Yk ) in (7). As γk is a binary variable, this
bility p(γk , Γ
is easily solved by letting
b k−1 |Yk )
γ̂k = argmaxγk p(γk , Γ
(8)
(
b k−1 |Yk ) > p(γk = 0, Γ
b k−1 |Yk ),
1, if p(γk = 1, Γ
=
0, otherwise.

4

In view of (7), we further obtain that
b k−1 |Yk )
p(γk = 1, Γ

(9)
b k−1 |Yk )
p(γk = 0, Γ
b k−1 , Yk−1 )p(γk = 1|Γ
b k−1 , Yk−1 )
p(yk |γk = 1, Γ
=
.
b k−1 , Yk−1 )p(γk = 0|Γ
b k−1 , Yk−1 )
p(yk |γk = 0, Γ

b k−1 , Yk−1 ) and
It remains to compute p(yk |γk , Γ
b
b k−1 , Yk−1 ), we
p(γk |Γk−1 , Yk−1 ). To obtain p(yk |γk , Γ
denote the probability density function of the Gaussian
distribution with mean µ and covariance matrix σ 2 by
b k is known, it follows from the IEKF that
N (µ, σ 2 ). Once Γ
b k−1 )
p(yk |γk , Yk−1 , Γ
(
N (h(x̂k|k−1 ), Ck Σk|k−1 CkT + R),
≈
N (0, R),

(10)
γk = 1,
γk = 0,

where x̂k|k−1 and Σk|k−1 are computed in (4) by replacing
b k−1 .
Γk−1 with Γ
If the prior probability distribution of γk is also available,
b k−1 ). Two common
we are able to easily compute p(γk |Yk−1 , Γ
cases are illustrated below.
(a) {γk } is a Bernoulli process with parameter θ, then
b k−1 ) = p(γk = 1) = θ.
p(γk = 1|Yk−1 , Γ
(b) {γk } is a Markov process with the transition probability
matrix [22]


1−q
q
p(γk = j|γk−1 = i) =
, i, j ∈ {0, 1},
p
1−p
b k−1 ) = p(γk |γ̂k−1 ).
then p(γk |Yk−1 , Γ
b k−1 |Yk ), and thus obtain
Finally, we can compute p(γk , Γ
the MAP estimate γ̂k by (8) and (9). The estimate x̂k|k is
then updated from (4) by replacing γk with γ̂k .
Overall, the BKF-I is summarized in Algorithm 1. Compared with the IEKF, we need to further compute the MAP
estimate of γk . The good news is that its computational
complexity is still comparable to that of the IEKF.
Algorithm 1 Bayesian Kalman filter I
1) Initialization: Let x̂0|−1 = x̄0 , Σ0|−1 = Σ0 .
2) MAP estimate of γk : Given yk , x̂k|k−1 , Σk|k−1 and
b k−1 , compute (9) and approximately obtain the MAP
Γ
estimate of γk by (8).
3) Measurement update: Given yk and γ̂k , update the state
estimate and its estimation error covariance matrix via (4)
by replacing γk with γ̂k .
4) Time update: Update the state prediction via (5).

B. Bayesian Kalman Filter II
In this subsection, we design the second suboptimal filter
by approximately computing the posterior distribution of the
measurement loss process {γk }.
Consider the minimum variance estimate, which is given by
the conditional expectation
Z
x̂k|k = E[xk |Yk ] = xk p(xk |Yk )dxk .

Since the system is nonlinear and non-Gaussian, the EKF
cannot be applied to compute the optimal estimate. Since the
posterior distribution p(xk |Yk ) is not tractable, the optimal
filter is unavailable or cannot be implemented recursively. Note
that a recursive suboptimal filter is necessary in practice. To
this end, we adopt the following approximation
p(xk |Yk ) ≈ p(xk |yk , x̂k|k−1 ).
That is, we use x̂k|k−1 to synthesize the information contained in Yk−1 to compute the posterior distribution. Then, it
follows that
x̂k|k = E[xk |Yk ] ≈ E[xk |yk , x̂k|k−1 ],
Σk|k = E[(xk − E(xk |Yk ))(xk − E(xk |Yk ))T |Yk ]
≈ E[(xk − x̂k|k )(xk − x̂k|k )T |x̂k|k−1 , yk ].
Besides, one can derive that
p(xk |x̂k|k−1 , yk )
=

1
X

p(xk |x̂k|k−1 , γk = i, yk )p(γk = i|x̂k|k−1 , yk )

i=0

= λk p(xk |γk = 1, x̂k|k−1 , yk )
+ (1 − λk )p(xk |γk = 0, x̂k|k−1 , yk )
where λk denotes the probability of event {γk = 1} conditioned on x̂k|k−1 and yk , and is given by
λk = p(γk = 1|x̂k|k−1 , yk )
p(yk |γk = 1, x̂k|k−1 )p(γk = 1|x̂k|k−1 )
,
= P1
j=0 p(yk |γk = j, x̂k|k−1 )p(γk = j|x̂k|k−1 )

(11)

where p(γk = 1|x̂k|k−1 ) is approximated by λk−1 .
It remains to compute p(yk |γk , x̂k|k−1 ). Consider a Gaussian approximation p(xk |x̂k|k−1 ) ≈ N (x̂k|k−1 , Σk|k−1 ), it
yields that
p(yk |γk , x̂k|k−1 )
(
N (h(x̂k|k−1 ), Ck Σk|k−1 CkT + R), γk = 1,
≈
N (0, R),
γk = 0.

(12)

Adopting the extended Kalman filtering technique [23], we
obtain that
(
N (µ1 , σ12 ), γk = 1,
p(xk |γk , x̂k|k−1 , yk ) ≈
N (µ2 , σ22 ), γk = 0,
where the mean and covariance are respectively given by
µ1 = x̂k|k−1 + Kk (yk − h(x̂k|k−1 )),
σ12 = Σk|k−1 − Σk|k−1 CkT (Ck Σk|k−1 CkT + R)−1 Ck Σk|k−1 ,
µ2 = x̂k|k−1 ,
σ22 = Σk|k−1 .
(13)

5

Combining the above, the measurement update for the
minimum variance estimate is approximately given as
Z
x̂k|k ≈ xk p(xk |yk , x̂k|k−1 )dxk
Z
= λk xk p(xk |γk = 1, x̂k|k−1 , yk )dxk
Z
+ (1 − λk ) xk p(xk |γk = 0, x̂k|k−1 , yk )dxk
= λk µ1 + (1 − λk )µ2
= x̂k|k−1 + λk Kk (yk − h(x̂k|k−1 )).
where the last equality follows from (13).
To run the algorithm recursively, we still need to update
Σk|k as well, which is derived below.
Σk|k
≈ E[(xk − E(xk ))(xk − E(xk ))T |x̂k|k−1 , yk ]
T
= E[xk xT
k |x̂k|k−1 , yk ] − E[xk |x̂k|k−1 , yk ]E[xk |x̂k|k−1 , yk ]
T
≈ E[xk xT
k |x̂k|k−1 , yk ] − x̂k|k x̂k|k .

The first term of Σk|k is further written as
Z
T
E[xk xk |x̂k|k−1 , yk ] = xk xT
k p(xk |x̂k|k−1 , yk )dxk
Z
= λ k xk xT
k p(xk |γk = 1, x̂k|k−1 , yk )dxk +
Z
(1 − λk ) xk xT
k p(xk |γk = 0, x̂k|k−1 , yk )dxk
= λk E[xk xT
k |γk = 1, x̂k|k−1 , yk ]
+ (1 − λk )E[xk xT
k |γk = 0, x̂k|k−1 , yk ]
2
T
2
= λk (µ1 µT
1 + σ1 ) + (1 − λk )(µ2 µ2 + σ2 ).

Finally, Σk|k is updated by the following formula
Σk|k
T
≈ E[xk xT
k |x̂k|k−1 , yk ] − x̂k|k x̂k|k
2
T
2
= λk (µ1 µT
1 + σ1 ) + (1 − λk )(µ2 µ2 + σ2 )

− [λk µ1 + (1 − λk )µ2 )][λk µ1 + (1 − λk )µ2 )]T
= λk (1 − λk )(µ1 − µ2 )(µ1 − µ2 )T + λk σ12 + (1 − λk )σ22
= Σk|k−1 − λk Kk Ck Σk|k−1 + λk (1 − λk )
× [Kk (yk − h(x̂k|k−1 ))(yk − h(x̂k|k−1 ))T KkT ].
Overall, the BKF-II is summarized in Algorithm 2. In the
BKF-II, we use λk to represent the uncertainty induced by the
measurement losses. By (11), it is clear that if γk is known to
the estimator, then λk = γk , which reduces to the IEKF. As
γk is unknown in our case, λk plays the role of estimating γk .
Note again that approximation is adopted to derive the
recursive filter. However, from the storage and computation
points of view, this suboptimal filter may be better than
an optimal one in applications, since an optimal minimum
variance estimate conditioned on Yk is generally not tractable.

Algorithm 2 Bayesian Kalman filter II
1) Initialization: Let x̂0|−1 = x̄0 , Σ0|−1 = Σ0 .
2) Measurement update: The state estimate and its approximated error covariance matrix are updated as follows.
x̂k|k = x̂k|k−1 + λk Kk (yk − h(x̂k|k−1 )),
Σk|k = Σk|k−1 − λk Kk Ck Σk|k−1 + λk (1 − λk )
[Kk (yk − h(x̂k|k−1 ))(yk − h(x̂k|k−1 ))T KkT ]
where λk is given in (11).
3) Time update: Update the state prediction via (5).

IV. R AO -B LACKWELLISED PARTICLE F ILTER
Clearly, both the BKF-I and BKF-II are not the minimum
variance estimate and their performance cannot be guaranteed
due to the use of approximation. In this section, a numerical method called Rao-Blackwellised particle filter (RBPF)
is applied to approximately compute the minimum variance
estimate. While the RBPF is computationally more demanding
than the BKF-I and BKF-II, its estimation performance can be
improved by increasing the number of particles.
A. RB Particle Filter
Recall that a minimum variance filter of xk is expressed by
Z
ZZ
x̂k|k = xk p(xk |Yk )dxk =
xk p(xk , Γk |Yk )dxk dΓk .
(14)
Since p(xk , Γk |Yk ) is not Gaussian, it is impossible to
be analytically obtained, even for linear systems. Then, the
integral is not computable and we have to resort to a numerical
approach.
In this section, we choose a particle based algorithm to
approximate the integral in (14). The particle filter (PF)
is a powerful sampling based method to approximate any
probability distribution, and address the non-linearity and nonGaussianity problem. However, it is acknowledged that the
number of particles grows dramatically with the dimension of
the underlying random vector, which substantially increases
the computation load.
In our case, if the particles are directly used to approximate p(xk , Γk |Yk ), it will result in a significant waste of
particles since the unconditional distribution of the state for
the linearized model in (3) is deemed to be Gaussian, and
the non-linearity and/or non-Gaussianity only appear in the
measurement equation, which results from the unknown measurement losses. To well exploit this observation, we shall
adopt the Rao-Blackwellised particle filter (RBPF) and use all
the particles to approximate the conditional distribution of the
binary variable γk . Intuitively, a few particles are sufficient to
accomplish it. Furthermore, it follows from [19] that the RBPF
leads to better estimation performance than the standard PF.
To exposit it, it follows from the conditional probability
definition that
p(xk , Γk |Yk ) = p(xk |Γk , Yk )p(Γk |Yk )
where p(xk |Γk , Yk ) is an approximated Gaussian density and
is recursively computed via the IEKF.

6

However, p(Γk |Yk ) is difficult to obtain. As Γk is a binary
sequence, we use N particles {Γik }N
i=1 that are drawn from an
importance density q(Γk |Yk ) to approximate p(Γk |Yk ), i.e.,
p(Γk |Yk ) ≈

N
X

ωki δ(Γk − Γik ),

(15)

i=1

where δ(·) is the standard Dirac delta function and ωki is the
normalized particle weight associated with Γik , i.e.,
p(Γik |Yk )
ωki ∝
.
q(Γik |Yk )
Inserting (15) into (14), we obtain that
ZZ
x̂k|k =
xk p(xk , Γk |Yk )dxk dΓk
ZZ
=
xk p(xk |Γk , Yk )dxk p(Γk |Yk )dΓk
≈

=

N
X
i=1
N
X

ωki

Z

xk p(xk |Γik , Yk )dxk

ωki E[xk |Γik , Yk ] =

i=1

N
X

A nice feature of the above is that p(yk |Yk−1 , Γik ) is approximately a conditional Gaussian density and is recursively
computed by IEKF, i.e.,
p(yk |Yk−1 , Γik )
(
N (h(x̂ik|k−1 ), Ck Σik|k−1 CkT + R),
≈
N (0, R),

ωki x̂ik|k ,

i=1

where
= E[xk |Γik , Yk ].
Similarly, let Σik|k = E[(xk

Σk|k ≈
=

N
X
i=1
N
X

ωki E[(xk − x̂k|k )(xk − x̂k|k )T |Γik , Yk ]
(17)
ωki Σik|k .

It should be noted from (3) that both (16) and (17) can
be recursively computed by using the IEKF. Specifically, we
obtain that
x̂ik|k = x̂ik|k−1 + γki Kki (yk − h(x̂ik|k−1 )),
Σik|k = Σik|k−1 − γki Kki Ck Σik|k−1 .

B. Importance Density
If an importance density is chosen to factorize such that
(18)

one can obtain particles Γik ∼ q(Γk |Yk ) by augmenting each
existing particle Γik−1 ∼ q(Γk−1 |Yk−1 ) with the new state
γki ∼ q(γk |Γk−1 , Yk ) and obtain a recursive filter algorithm.
To elaborate it, we express p(Γk |Yk ) in the following form
p(Γk |Yk ) ∝ p(yk |Yk−1 , Γk )p(γk |Γk−1 , Yk−1 )p(Γk−1 |Yk−1 ).
Jointly with (18), it implies that
∝

p(yk |Yk−1 , Γik )p(γki |Γik−1 , Yk−1 )p(Γik−1 |Yk−1 )
q(γki |Γik−1 , Yk )q(Γik−1 |Yk−1 )

i
= wk−1

p(yk |Yk−1 , Γik )p(γki |Γik−1 , Yk−1 )
.
q(γki |Γik−1 , Yk )

If {γk } a Bernoulli process or Markov process, the update
equation for the particle weight in (19) becomes particularly
simple. Specifically,
ωki ∝

i
p(yk |Yk−1 , Γik )p(γki |γk−1
)
p(Γik |Yk )
i
=
ω
k−1
i
i
i
q(Γk |Yk )
q(γk |Γk−1 , Yk )

In the current situation, the number of particles is very
small since γk is binary. The RBPF with a resampling step
is summarized in Algorithm 3.
C. Fast Implementation

The remaining problem is how to recursively generate
particles {Γik } and compute their associated weights {ωki }.
This is resolved in next subsection.

wki

i
γki ∼ p(γk |γk−1
).

i
= ωk−1
p(yk |Yk−1 , Γik ).

i=1

q(Γk |Yk ) = q(γk |Γk−1 , Yk )q(Γk−1 |Yk−1 ),

N
1 + var(wk∗i )

where wk∗i = p(γki |Yk )/q(γk |Γik−1 , Yk ) is referred to as the
true importance weights.
i
The other is p(γk |γk−1
), which makes it easy to draw particles and compute the important weights. Since p(γk |Γk−1 , Yk )
is usually difficult to access, we choose q(γk |Γk−1 , Yk ) =
p(γk |γk−1 ) as the importance density in this work. That is,
the new particle is generated from the following distribution

x̂ik|k

− x̂k|k )(xk − x̂k|k )T |Γik , Yk ],
then the estimation error covariance matrix is given as

γki = 0.

To alleviate the degeneracy problem, which is key to the
success of PFs, there are two good choices of the importance density q(γk |Γik−1 , Yk ) in the literature [17]. One is
p(γk |Γik−1 , Yk ), which minimizes a suitable measure of the
degeneracy of the algorithm, i.e.,
Neff =

(16)

γki = 1,

(19)

The major computation of the RBPF is taken in the Output
step in Algorithm 3, which requires an IEKF iteration for
each triple (γki , x̂ik|k−1 , Σik|k−1 ). Thus, the computational complexity essentially is proportional to the number of particles.
When the number of particles is very large, one may observe
that there are many duplicate triples after the Resampling c),
i.e., there may be only m  N triples that take different
values. We can take advantage of this observation to reduce
the computational complexity of RBPF.
Instead of directly executing an IEKF iteration for each
triple (γki , x̂ik|k−1 , Σik|k−1 ), we first remove the duplicate
triples from the list and obtain a set of triples that take different
values, see step a) in Algorithm 4, and then do an IEKF
iteration for each of them, see step b) in Algorithm 4. Finally,
we set the updated value for the whole list, see step c).
Clearly, Algorithm 4 returns the same estimate and prediction as those of RBPF in Algorithm 3. It may significantly
reduce the computational complexity of the RBPF, especially
when the number of particles is large. Besides, the reduction
of the computational complexity reflects the effectiveness

7

Algorithm 3 Rao-Blackwellised particle filter
1) Initialization: For i = 1, ..., N. draw N particles γ0i from
the prior p(γ0 ) and let x̂i−1|0 = x0 , Σi−1|0 = Σ0 .
2) Importance sampling:
a) For i ∈ {1, ..., N } and k > 0, draw N particles from
i
the importance distribution γki ∼ p(γk |γk−1
).
b) For i ∈ {1, ..., N } and given yk , update the normalized
i
importance weights ωki ∝ ωk−1
p(yk |Yk−1 , Γik ).
3) Resampling:
a) Compute an estimate of the effective number of particles
beff = P 1
.
N
N
i 2
i=1 (ωk )
beff < NT , which is a prescribed threshold, perform
b) If N
i∗
resampling. Take N new samples (γki∗ , x̂i∗
k|k−1 , Σk|k−1 )
i
i
i
with replacement from the list (γk , x̂k|k−1 , Σk|k−1 ),
i ∈ {1, . . . , N } according to the probability distribution that Pr{γki∗ = γki } = ωki .
i∗
c) Let (γki , x̂ik|k−1 , Σik|k−1 ) = (γki∗ , x̂i∗
k|k−1 , Σk|k−1 ) and
i
ωk = 1/N for all i ∈ {1, . . . , N }.
4) Output
a) For each triple (γki , x̂ik|k−1 , Σik|k−1 ) and yk , do
Kki
x̂ik|k
Σik|k

=
=
=

Σik|k−1 CkT (Ck Σik|k−1 CkT + R)−1
x̂ik|k−1 + γki Kki (yk − h(x̂ik|k−1 ))
Σik|k−1 − γki Kk Ck Σik|k−1 .

b) Measurement update
PN
PN
x̂k|k = i=1 ωki x̂ik|k , Σk|k = i=1 ωki Σik|k .
c) Update the state prediction by (5), and set
x̂ik+1|k = f (x̂ik|k , uk ), Σik+1|k = Ak Σik|k ATk + Q.
5) Let k = k + 1.

of particles to some extent. However, when the number of
particles is small, the fast implementation may be not efficient
as it further requires to remove duplicate triples. This will be
illustrated in Section VI-A by an numerical example.
V. C OMPARISONS OF THE THREE FILTERS
This section provides some comparisons of the proposed
filters in terms of the estimation accuracy and computational
complexity. We shall further verify the major results of this
section by numerical examples in the next section.
A. Estimation Accuracy
If both f and h in (2) are linear, any filter cannot outperform
the optimal IKF, which however requires {γk } is known. The
performance of the RBPF can approach that of the IKF arbitrarily well by increasing the number of particles, and hence
has a predictable performance. This cannot be guaranteed for
the BKF-I and BKF-II due to approximations in deriving the
recursive filters. However, if yk − h(x̂k|k−1 ) is very different

Algorithm 4 Fast RBPF
The Output step in Algorithm 3 is replaced as follows.
a) For the list of triples (γkj , x̂jk|k−1 , Σjk|k−1 ), j ∈
{1, ..., N }, remove duplicate triples and obtain an index
set I such that there is no duplicate triple in the set
{(γki , x̂ik|k−1 , Σik|k−1 ), i ∈ I}.
b) For each triple (γki , x̂ik|k−1 , Σik|k−1 ), i ∈ I and yk , do
Kki = Σik|k−1 CkT (Ck Σik|k−1 CkT + R)−1
x̂ik|k = x̂ik|k−1 + γki Kki (yk − h(x̂ik|k−1 ))
Σik|k = Σik|k−1 − γki Kk Ck Σik|k−1
x̂ik+1|k = f (x̂ik|k , uk ),
Σik+1|k = Ak Σik|k ATk + Q.
c) For each j ∈ {1, ..., N }, set
(x̂jk+1|k , x̂jk|k , Σ̂jk+1|k , Σ̂jk|k )
=(x̂ik+1|k , x̂ik|k , Σ̂ik+1|k , Σ̂ik|k )
if (γkj , x̂jk|k−1 , Σjk|k−1 ) = (γki , x̂ik|k−1 , Σik|k−1 ) for some
i ∈ I.
d) Do measurement update
x̂k|k =

N
X

ωki x̂ik|k , Σk|k =

i=1

N
X

ωki Σik|k ,

i=1

and update the state prediction by (5).

from 0, the measurement loss γk can be estimated correctly
with high probability. Then, both the BKF-I and BKF-II are
likely to approach the IKF. Compared to the BKF-II, the
BKF-I is more intuitive and easier to implement. The BKFII appears to be more stable as λk retains the probability
information of {γk = 1}. Since the system in (2) is nonlinear
under the unknown measurement losses, these observations
are impossible to prove in theory. We shall provide numerical
examples for validation in the next section.
If either f or h in (2) is nonlinear, our filters first linearize
the nonlinear function and then applies the corresponding
algorithms to the system in (3). Due to the approximation
in the linearization step, the IEKF is not optimal even it uses
the true Γk , and all our filters might have better performance
than the IEKF. Note that even when there is no measurement
loss, i.e., γk = 1 for all k, the estimation accuracy of the EKF
cannot be guaranteed. By using a sufficiently large number
of particles, the estimation accuracy of the RBPF may even
outperform that of the IEKF. Since approximation is used
twice in both the BKF-I and BKF-II, it is expected that the
IEKF and RBPF usually outperform the BKF-I and BKF-II.

B. Computational Complexity
We provides Table I to compare the computational complexities of the three proposed filters. The computation for
linearization in deriving the system in (3) is omitted since all
the proposed filters require it.

8

TABLE I
C OMPUTATIONAL COMPLEXITIES IN ONE ITERATION
Filters
Number of times of
Gaussian PDF evaluations
Number of times of
matrix multiplications
Computational complexity
a
b

IEKF

RBPF

BKF-I

BKF-II

0

2N

2

2

7

O(N )

7

9

O(n3 )

O(N n3 )

O(n3 )

O(n3 )

N is the number of particles in RBPF.
n is the dimension of state vectors.

The computational complexity can be quantified by the
number of times of evaluating the Gaussian probability density
function (PDF) and the matrix multiplications in one iteration.
It is apparent that the RBPF has the highest computational
complexity which is proportional to the number of particles
and is far more than that of BKF-I and BKF-II. When
N is large, the fast RBPF might significantly reduce the
computational complexity, and will be illustrated by numerical
examples in the next section.
The BKF-I and BKF-II have similar computational complexities, and both require additionally two Gaussian PDF
evaluations than that of the IEKF. Generally speaking, the
computational complexity of evaluating a Gaussian PDF is
O(n3 ), where n is the dimension of the random vector
and is the dimension of xk here. Since the computational
complexity of multiplication of two n × n matrices is also
O(n3 ) in general, both the BKF-I and BKF-II have the same
computational complexity O(n3 ), and are similar to that of
the IEKF and EKF.
C. Filter Selection
The filter selection depends on many factors, and usually
requires to find a good tradeoff between the estimation accuracy and the computational complexity. The tradeoff between
the estimation accuracy and the computational complexity is
illustrated in Section VI-A. If the computation resource is
abundant, it is recommended to use the RBPF with a sufficient
number of particles to obtain a good estimation accuracy. If
not, we may consider to use the BKF-I or BKF-II. Note that
whether one filter is always superior to the other two filters is
not conclusive.
VI. A PPLICATION E XAMPLES
In this section we adopt the proposed filters to solve three
estimation problems. The first one is a state estimation problem of a linear system with randomly unknown measurement
losses. In this example, it does not need to do linearization
as in (3). Then, the approximation in both BKF-I and BKF-II
is solely used to handle the measurement loss process {γk }
in a recursive way. We compare the performance of proposed
filters and show the relations between the number of particles,
the computation time, and the root-mean-square error (RMSE)
of the RBPF.
The second one is a target tracking problem, which is a
nonlinear system with measurement losses, and illustrates the

effectiveness of our filters. We also provide a case where our
filters do not work very well and give some explanations of
it in this example. The last example is a quadrotor’s path
control problem with random measurement losses, where the
quadrotor’s states estimated by our filters are used in feedback
to control a quadrotor. This requires precise state estimates,
otherwise the system may diverge quickly. The result shows
the quadrotor can complete the task well under all our filters.
A. Linear System with Random Measurement Losses
Consider the following system


0.6 0.4
xk+1 =
x + wk ,
0.1 0.9 k


yk = γk 1 −2 xk + vk ,

(20)

where wk and vk are independent Gaussian noise with zero
means, and the covariance matrices are both identity matrices.
{γk } is an i.i.d. Bernoulli process with P r{γk = 0} = p. The
initial state is also a Gaussian random vector with zero mean
and identity covariance matrix. Since the open-loop poles of
the system (20) are λ1 = 0.5 and λ2 = 1, it follows from [13]
that the optimal IKF is expected to be mean square stable if
the measurement loss level p is strictly less than one.
To validate the proposed filters, we use the same sequence
of noisy observations Y200 for each filter, and the RBPF is
implemented by using 20 particles. We also compare them
with the optimal IKF, which requires to know the measurement
loss, and the filter that does not consider the measurement
loss, i.e., we simply set γk = 1 for all k in (4) and (5),
which is denoted as ‘KF’ in Fig. 1. To compute the RMSE
of each estimate where the RMSE of an estimate x̂k|k at
time k is defined as (E[kxk − x̂k|k k2 ])1/2 , we adopt the
Monte Carlo method with 500 independent experiments. Fig.
1 illustrates that the estimation accuracy of the proposed three
filters comes close to that of the optimal IKF, and is much
better than that of the ‘KF’. Moreover, the larger probability
of the measurement losses, the better of the improvement of
the estimation accuracy over the ‘KF’. This clearly shows the
advantage of considering the measurement losses in the filter
design.
In Fig. 2, we illustrate how the number of particles affects
the computational complexity and the estimation accuracy of
the RBPF. One can easily observe that the computational
complexity per iteration of the RBPF is proportional to the
number of particles, which validates the result for the RBPF
in Table I. It also verifies that the fast implementation of
RBPF in Algorithm 4 can greatly reduce the computation load
when N is large, while it is slightly slower than the RBPF
for a small number of particles. The reason is that the fast
RBPF further requires ‘find and set’ steps, see steps a) and
c) in Algorithm 4, which are dominated by the computation
for executing the IEKF iteration if the number of particles is
large. Moreover, Fig. 2 reveals that the sum of RMSE reduces
quickly when N is smaller, and slowly when N is large. This
suggests the existence of a critical number of particles that
is vital to the tradeoff between the estimation accuracy and
the computational complexity. Note that the critical number is
also related to the levels of random measurement losses.

9

TABLE II
PARAMETERS IN THE TARGET TRACKING SYSTEM

12
p=0.8
p=0.6
p=0.2

Sum of RMSE

10

Parameter

Description

Value

6

τ
σm

0.01s
4m/s2

4

α

Sampling period
Variance of target acceleration
Reciprocal of the
maneuver (acceleration) time constant
Standard deviation of
angular measurement noise
Standard deviation of
distance measurement noise
Probability of measurement losses

8

2

σφ

0
KF

IKF

BKF-I

BKF-II

RBPF

σr

Time per iteration (ms)

Fig. 1. The sum of RMSE over all the 200 time steps under different levels
of measurement losses.

5◦
5m
0.1, 0.3, 0.5, 0.7

2
1.5
1
0.5

RBPF
Fast RBPF

0
0

5

10

15

20

25

30

35

40

45

50

12
p=0.8

10

where τ is the sampling period.
Here the measurement loss process {γk } is also a Bernoulli
process, which represents the loss of the radar measurement.
(i)
The input noise wk and the output noise vk = [vkr , vkφ ]T are
independently Gaussian distributed with zero mean. As in [24]
and [20], the covariance matrix of wk is set to

 5
τ /20 τ 4 /8 τ 3 /6
2  4
τ /8 τ 3 /3 τ 2 /2
Q = 2ασm
τ 3 /6 τ 2 /2
τ
and the covariance matrix of vk is given as
 2

σ
0
R= r
.
0 σφ2

p=0.6

Sum of RMSE

p

1

p=0.9

(i)

8

The initial state x0 is a Gaussian random vector with mean
= [10, 0, 0]T , i ∈ {1, 2} and covariance
 2

σr
σr2 /τ
0
Σ0 = σr2 /τ 2σr2 /τ 2 0 .
0
0
0

(i)
x̄0
6
0

5

10

15

20

25

30

35

40

45

50

Number of Particles

Fig. 2. Computation time per iteration and RMSE of the RBPF with respect
to the number of particles. The time is measured using Matlab.

As in (3), we obtain the linearized measurement equation



(1)

pk|k−1

B. Target Tracking
In this example, we test our filters on a target tracking problem, which further requires the linearization in (3). Consider
a 2D target tracking problem where a radar is positioned on
the ground to measure the distance and the angle of a moving
target to the base station. The radar measurements are then
sent to an unmanned tracking aircraft via a lossy network,
which is subject to unknown measurement losses.
The motion and measurement equations are described by
(i)

(i)

(i)

(i)

xk+1 = F xk + uk + wk , i ∈ {1, 2}
   r
v
r
yk = γk
+ φk
φ
vk
(i)

(i)

(i)

(i)

(21)

where xk = [pk , ṗk , p̈k ]T denotes the target state at time
step k, including the target position, speed and acceleration in
(i)
the i-th dimension. The system input uk is given as a prior,
and the system matrices are


q

(1) 2
(2) 2
1 τ τ 2 /2  
(p
)
+
(p
)
r
k
k
,
τ ,
F = 0 1
=
(2)
pk
φ
arctan
(1)
0 0
1
p
k

 rk|k−1
yk = γk 
 −p(2)

k|k−1

2
rk|k−1

(2)

pk|k−1
rk|k−1

(1)

 pk 
 + vk + zk ,

(1)
−pk|k−1   (2) 
pk
r2
k|k−1

where zk is a function of pk|k−1 . All the parameters are of
the same as those in [24], [25], see Table II, and the initial
(i)
estimate is set to x̂0|0 = [10, 0, 0]T , i ∈ {1, 2}.
To track the system, we use the proposed three filters and
the IEKF, which requires to know the measurement loss, and
run 1500 independent experiments for each filter to compute
the RMSE of the position per time step. The RBPF adopts
200 particles.
From Fig. 3 and Fig. 4, one can observe that the RMSE
of the RBPF is generally smaller than that of the BKF-I
and BKF-II, and the BKF-I appears to outperform the BKFII in this example. However, we cannot conclude that the
BKF-I is better than the BKF-II. Note that the computational
complexities of these two filters come close to each other.
Clearly, the larger the probability of random measurement
losses, the larger the RMSE as expected.
It is interesting to observe that the estimation accuracy of the
RBPF is better than that of the IEKF in many cases, especially

10

140
120
100

RMSE

TABLE III
PARAMETERS IN THE QUADROTOR ’ S PATH CONTROL

RBPF
BKF-I
BKF-II
IEKF

80
60

Parameter

Description

Value

m
Jxx , Jyy

Quadrotor mass
Moments of inertia
Standard deviation of
angular measurement noise
Standard deviation of
position measurement noise
Control frequency
Measurements’ error probability

1.023kg
9.5g·m2

σφ , σθ , σψ

40

σX , σY , σZ

20

T
p

0
p=0.1

p=0.3

p=0.5

p=0.7

Fig. 3. The RMSE of position under different levels of random measurement
losses.
30
IEKF

BKF-I

BKF-II

RBPF

20

0.1◦
2m
100Hz
0.2

from different systems because of the non-linearity of the
systems. One can not conclude that one filter always performs
better than another. Usually, the RBPF outperforms the other
filters by increasing the number of particles.

10
0
0

2

4

6

8

10

12

14

16

18

20

RMSE (m)

100

50

0
0

2

4

6

8

0

2

4

6

8

10

12

14

16

18

20

10

12

14

16

18

20

400

200

0

Time (s)

Fig. 4. The RMSE of position under different levels of random measurement
losses, e.g., p = 0.1, p = 0.5 and p = 0.8.

when the probability of measurement losses p is large. It is
known that IKF is optimal when both f and h in (2) are linear
functions and γk is known. Under this case, the proposed filter
cannot perform better than the IKF. However, the IEKF is not
optimal if either f or h is nonlinear. In such cases, the RBPF
may achieve a better performance than the IEKF due to the
non-linearity in the measurement equation (21).
Fig. 5a and Fig. 5b show the true and estimated trajectories
in one experiment when p = 0.5 and p = 0.7 respectively. We
can observe that our filters track the target well, especially
under the low probability of measurement losses.
Fig. 5c indicates that the proposed filters perform unsatisfactorily under p = 0.9 and the poor initial estimate, i.e.,
(i)
x̂0|0 = [200, 0, 0]T , i ∈ {1, 2}. In this case, all the filters track
the trajectory badly, including the IEKF. This is due to the
non-linearity of the measurement equation and the poor initial
state estimate. Once the filters have an inaccurate estimate at
time k, the linearization step will further introduce estimation
errors and tends to a bad estimate. Similar to the standard EKF,
a good initial estimate is essential for nonlinear systems.
To summarize, the performance of the proposed filters varies

C. Quadrotor’s Path Control
Finally, we test our filters on a quadrotor’s path control
problem where the state estimate is used for the feedback
design. A quadrotor is a very maneuverable unmanned aircraft
with four rotors, and has been popular both in research and
real applications over last years. The structure of quadrotors is
shown in Fig. 6 [26], where [X, Y, Z] and [φ, θ, ψ] are position
and angular of the quadrotor in the inertia frame, respectively.
[U, V, W ] and [P, Q, R] are position and angular velocities in
the body frame. Due to page limitation, we omit the dynamical
model of quadrotors, see [27]–[31] for more details.
The control architecture is illustrated in Fig. 7. The IMU
and GPS block provides the noisy measurements of [φ, θ, ψ],
[P, Q, R] and [X, Y, Z], [U, V, W ] respectively. Due to bad
GPS data and/or the unstable IMU, an estimator may only
receive the true measurements intermittently. To this end,
the proposed filters are applied to this situation and the
estimated state are further sent to the outer and inner loop PID
controller to control both the quadrotor’s position and attitude
respectively. This is realized by building a Matlab Simulink
model based on a Matlab package called Quadcopter Dynamic
Modeling and Simulation [26].
Our goal is to control the quadrotor to follow a given
triangle path in space, see Fig. 8 (black line), and the projected
trajectories in X,Y and Z directions are shown in black dashed
lines in Fig. 9. We consider that the IMU and GPS measurements are sent as a single packet, and they will be either
simultaneously lost or received by the estimator. Moreover,
the probability of measurement losses is set to be p = 0.2.
Parameters of the PID controller are selected to be the same
as those in [26], which however does not consider the problem
of measurement losses, and are provided in Table III.
The trajectories of the closed-loop system are shown in Fig.
8 and Fig. 9, where the reference path represents the desired
trajectory. The trajectory of ‘EKF’ denotes the controlled path
of using the standard EKF without any measurement loss, and
‘IEKF’ is that of using the IEKF to estimate the state, which
relies on the true value of γk at each time step. This is different
from ‘BKF-I’, ‘BKF-II’, and ‘RBPF’ where the packet loss

11

IEKF
BKF-I
BKF-II
RBPF
real trajectory
start point

0

Y-coordinate (m)

Y-coordinate (m)

0
-20
-40
IEKF
BKF-I
BKF-II
RBPF
real trajectory
start point

-60
-80

-200

-150

-100

-50

-50

IEKF
BKF-I
BKF-II
RBPF
real trajectory
start point

6000
4000

Y-coordinate (m)

20

-100

-150

2000
0
-2000
-4000
-6000

0

-150

-100

X-coordinate (m)

-50

0

50

0

2000

X-coordinate (m)

(a)

4000

6000

8000

X-coordinate (m)

(b)

(c)

Fig. 5. One run trajectory.

Reference path
EKF
IKF
BKF-I
BKF-II
RBPF

Z (m)

3

2

1

0
3
2
1

Fig. 6. The structure of a quadrotor [26].

0
Task

Given Outer loop
PID
path controller

Inner loop
PID
controller

Y (m)
Motors

-3

-2

Quadrotor
Dynamics

-1

0

1

2

3

X (m)

Fig. 8. The trajectory of the controlled quadrotor in space under different
filters. The red triangle mark is the end point of the trajectories.

Estimators

IMU & GPS

Fig. 7. The control architecture of the quadrotor’s path control system.

process {γk } is unknown to the estimator. For comparison, all
filters are tested under the same input noise, measurement loss
process {γk }, initial conditions and covariance matrices.
In view of Fig. 8 and Fig. 9, it is not difficult to observe that
all the proposed filters can fulfill the tracking task and their
performance comes close to that of the standard EKF without
any measurement loss. In contrast, we are unable to control
the quadrotor to follow the desired trajectory if we directly
apply the EKF update to all time steps without considering
the measurement losses, see Fig. 10, where the tracking error
diverges.
VII. CONCLUSION
In this paper, we have designed three suboptimal filters
to deal with the state estimation problem with unknown

measurement losses. Among these filters, the BKF-I and the
BKF-II were established by using the Bayesian point of view
and the IEKF. Some approximations were made to obtain
recursive forms. The RBPF is a particle filter based numerical
method with a relatively small number of particles. All these
proposed filters were applied to three different application
problems, showing their effectiveness. Future work will focus
on the theoretically convergence of the proposed filters.
R EFERENCES
[1] R. Teo, D. M. Stipanovic, and C. J. Tomlin, “Decentralized spacing
control of a string of multiple vehicles over lossy datalinks,” IEEE
Transactions on Control Systems Technology, vol. 18, no. 2, pp. 469–
473, 2010.
[2] M. N. Soltani, R. Izadi-Zamanabadi, and R. Wisniewski, “Reliable
control of ship-mounted satellite tracking antenna,” IEEE Transactions
on Control Systems Technology, vol. 19, no. 1, pp. 221–228, 2011.
[3] D. Choukroun and J. L. Speyer, “Mode estimation via conditionally
linear filtering: Application to gyro failure monitoring,” Journal of
Guidance, Control, and Dynamics, vol. 35, no. 2, pp. 632–644, 2012.
[4] X. Yu and J. Jiang, “Hybrid fault-tolerant flight control system design
against partial actuator failures,” IEEE Transactions on Control Systems
Technology, vol. 20, no. 4, pp. 871–886, 2012.

12

X (m)

2

[11]
0
-2

[12]
0

5

10

15

20

[13]

Y (m)

3
2

[14]

1
0

[15]
0

5

10

15

20

[16]

Z (m)

3

[17]

2
Reference path

1
IKF

BKF-I

EKF
BKF-II

RBPF

0
0

5

10

15

20

[18]

Time (s)
[19]
Fig. 9. The trajectory of the controlled quadrotor under different filters in X,
Y, and Z directions.

X (m)

[20]
400

[21]

200
0
0

2

4

6

8

10

12

14

16

18

20

[23]

100

Y (m)

[22]

[24]
50

[25]

0
0

2

4

6

8

10

12

14

16

18

20

[26]
Z (m)

-2.5

[27]
[28]

-3
0

2

4

6

8

10

12

14

16

18

20

Time (s)

(a)

[29]

Fig. 10. Tracking error of the EKF without considering the measurement
losses.
[30]
[5] M. Joerger and B. Pervan, “Kalman filter-based integrity monitoring
against sensor faults,” Journal of Guidance, Control, and Dynamics,
vol. 36, no. 2, pp. 349–361, 2013.
[6] R. C. Avram, X. Zhang, and J. Muse, “Nonlinear adaptive fault-tolerant
quadrotor altitude and attitude tracking with multiple actuator faults,”
IEEE Transactions on Control Systems Technology, 2017.
[7] D. Zhang and D. Ionescu, “A new method for measuring packet loss
probability using a Kalman filter,” IEEE Transactions on Instrumentation
and Measurement, vol. 58, no. 2, pp. 488–499, 2009.
[8] N. Kottenstette, J. F. Hall, X. Koutsoukos, J. Sztipanovits, and P. Antsaklis, “Design of networked control systems using passivity,” IEEE Transactions on Control Systems Technology, vol. 21, no. 3, pp. 649–665,
2013.
[9] B. L. Reed and F. S. Hover, “JLS-PPC: A jump linear system framework for networked control,” IEEE Transactions on Control Systems
Technology, vol. 25, no. 3, pp. 924–939, 2017.
[10] B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M. I. Jordan,

[31]

and S. S. Sastry, “Kalman filtering with intermittent observations,” IEEE
Transactions on Automatic Control, vol. 49, no. 9, pp. 1453–1464, 2004.
J. Wu, G. Shi, B. D. Anderson, and K. H. Johansson, “Kalman filtering
over Gilbert-Elliott channels: Stability conditions and critical curve,”
IEEE Transactions on Automatic Control, 2017.
K. You, N. Xiao, and L. Xie, Analysis and Design of Networked Control
Systems. Springer, 2015.
K. You, M. Fu, and L. Xie, “Mean square stability for Kalman filtering
with Markovian packet losses,” Automatica, vol. 47, no. 12, pp. 2647–
2657, 2011.
E. R. Rohr, D. Marelli, and M. Fu, “Kalman filtering with intermittent
observations: On the boundedness of the expected error covariance,”
IEEE Transactions on Automatic Control, vol. 59, no. 10, pp. 2724–
2738, 2014.
R. Van Der Merwe, A. Doucet, N. De Freitas, and E. Wan, “The unscented particle filter,” in Conference on Neural Information Processing
Systems, vol. 2000. Denver, CO, USA, 2000, pp. 584–590.
S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear
estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, 2004.
M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial
on particle filters for online nonlinear/non-Gaussian Bayesian tracking,”
IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188,
2002.
A. Doucet, N. De Freitas, and N. Gordon, “An introduction to sequential
Monte Carlo methods,” in Sequential Monte Carlo Methods in Practice.
Springer, 2001, pp. 3–14.
A. Doucet, N. De Freitas, K. Murphy, and S. Russell, “RaoBlackwellised particle filtering for dynamic Bayesian networks,” in
Proceedings of the Sixteenth Conference on Uncertainty in Artificial
Intelligence. Morgan Kaufmann Publishers Inc., 2000, pp. 176–183.
F. Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell, J. Jansson,
R. Karlsson, and P.-J. Nordlund, “Particle filters for positioning, navigation, and tracking,” IEEE Transactions on Signal Processing, vol. 50,
no. 2, pp. 425–437, 2002.
J. Zhang and K. You, “Kalman filtering with unknown sensor measurement losses,” IFAC-PapersOnLine, vol. 49, no. 22, pp. 315–320, 2016.
M. Huang and S. Dey, “Stability of Kalman filtering with Markovian
packet losses,” Automatica, vol. 43, no. 4, pp. 598–607, 2007.
B. D. Anderson and J. B. Moore, Optimal Filtering. Courier Corporation, 2012.
R. A. Singer, “Estimating optimal tracking filter performance for manned
maneuvering targets,” IEEE Transactions on Aerospace and Electronic
Systems, vol. 6, no. 4, pp. 473–483, 1970.
G. R. Curry, “Radar system performance modeling,” Artech House,
2005.
“Quadcopter dynamic modeling and simulation,” https://github.com/
dch33/Quad-Sim/blob/master/README.md.
D. W. Kun and I. Hwang, “Linear matrix inequality-based nonlinear
adaptive robust control of quadrotor,” Journal of Guidance, Control,
and Dynamics, pp. 996–1008, 2015.
J. L. Crassidis and F. L. Markley, “Three-axis attitude estimation
using rate-integrating gyroscopes,” Journal of Guidance, Control, and
Dynamics, pp. 1–14, 2016.
J. J. Dougherty, H. El-Sherief, and D. S. Hohman, “Use of the Global
Positioning System for evaluating inertial measurement unit errors,”
Journal of Guidance, Control, and Dynamics, vol. 17, no. 3, pp. 435–
441, 1994.
R. Mahony, V. Kumar, and P. Corke, “Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor,” IEEE Robotics Automation
Magazine, vol. 19, no. 3, pp. 20–32, Sept 2012.
T. Bresciani, “Modelling, identification and control of a quadrotor
helicopter,” MSc Theses, 2008.

13

PLACE
PHOTO
HERE

Jiaqi Zhang received the B.S. degree from the
School of Electronic and Information Engineering,
Beijing Jiaotong University, Beijing, China, in 2016.
He is currently pursuing the Ph.D. degree at the
Department of Automation, Tsinghua University,
Beijing, China. His research interests include networked control systems, distributed optimization and
their applications.

Keyou You received the B.S. degree in statistical
science from Sun Yat-sen University, Guangzhou,
China, in 2007 and the Ph.D. degree in electrical and
electronic engineering from Nanyang Technological
PLACE
University (NTU), Singapore, in 2012. After briefly
PHOTO
working as a Research Fellow at NTU, he joined
HERE
Tsinghua University, Beijing, China in 2012 where
he is currently an Associate Professor with the Department of Automation. He held visiting positions
at Politecnico di Torino, Turin, Italy, the Hong Kong
University of Science and Technology, Hong Kong,
and the University of Melbourne, Parkville, VIC, Australia. His current
research interests include networked control systems, distributed algorithms,
and their applications.
Dr. You was a recipient of the Guan Zhaozhi Award at the 29th Chinese
Control Conference in 2010, the CSC-IBM China Faculty Award in 2014,
and the National Science Fund for Excellent Young Scholars in 2017. He was
nominated for the National 1000-Youth Talent Program of China in 2014.

Lihua Xie received the B.E. and M.E. degrees in
electrical engineering from Nanjing University of
Science and Technology in 1983 and 1986, respectively, and the Ph.D. degree in electrical engineerPLACE
ing from the University of Newcastle, Australia, in
PHOTO
1992. Since 1992, he has been with the School
HERE
of Electrical and Electronic Engineering, Nanyang
Technological University, Singapore, where he is
currently a professor and Director, Delta-NTU Corporate Laboratory for Cyber-Physical Systems. He
served as the Head of Division of Control and
Instrumentation from July 2011 to June 2014. He held teaching appointments
in the Department of Automatic Control, Nanjing University of Science and
Technology from 1986 to 1989 and Changjiang Visiting Professorship with
South China University of Technology from 2006 to 2011.
Dr. Xie’s research interests include robust control and estimation, networked
control systems, multi-agent networks, localization and unmanned systems.
He is an Editor-in-Chief for Unmanned Systems and an Associate Editor
for IEEE Transactions on Network Control Systems. He has served as an
editor of IET Book Series in Control and an Associate Editor of a number
of journals including IEEE Transactions on Automatic Control, Automatica,
IEEE Transactions on Control Systems Technology, and IEEE Transactions
on Circuits and Systems-II. He is an elected member of Board of Governors,
IEEE Control System Society (Jan 2016-Dec 2018). Dr. Xie is a Fellow of
IEEE and Fellow of IFAC.

