arXiv:1503.02878v1  [cs.IT]  10 Mar 2015
Mobile Node Localization via Pareto Optimization:
Algorithm and Fundamental Performance
Limitations
Alessio De Angelis, Member, IEEE, Carlo Fischione, Member, IEEE
Abstract—Accurate estimation of the position of network
nodes is essential, e.g., in localization, geographic routing, and
vehicular networks. Unfortunately, typical positioning techniques
based on ranging or on velocity and angular measurements
are inherently limited. To overcome the limitations of specific
positioning techniques, the fusion of multiple and heterogeneous
sensor information is an appealing strategy. In this paper, we
investigate the fundamental performance of linear fusion of mul-
tiple measurements of the position of mobile nodes, and propose
a new distributed recursive position estimator. The Cramér-
Rao lower bounds for the parametric and a-posteriori cases
are investigated. The proposed estimator combines information
coming from ranging, speed, and angular measurements, which is
jointly fused by a Pareto optimization problem where the mean
and the variance of the localization error are simultaneously
minimized. A distinguished feature of the method is that it
assumes a very simple dynamical model of the mobility and
therefore it is applicable to a large number of scenarios providing
good performance. The main challenge is the characterization of
the statistical information needed to model the Fisher information
matrix and the Pareto optimization problem. The proposed analy-
sis is validated by Monte Carlo simulations, and the performance
is compared to several Kalman-based filters, commonly employed
for localization and sensor fusion. Simulation results show that
the proposed estimator outperforms the traditional approaches
that are based on the extended Kalman filter when no assumption
on the model of motion is used. In such a scenario, better
performance is achieved by the proposed method, but at the
price of an increased computational complexity.
Index Terms—Sensor Fusion, Positioning, Distributed Models,
Networks, Optimization, Cramér Rao Lower Bound.
I. INTRODUCTION
In many sensor network applications it is desirable to
accurately estimate the position of mobile wireless nodes [1].
Relevant applications are vehicular traffic monitoring, asset
tracking, process monitoring, and control of autonomous
agents. As an example, accurate position information is crucial
for emergency personnel and first responders, see, e.g., [2].
The commonly used Global Navigation Satellite Systems
(GNSSs) provide position information with an accuracy of
about 3 to 10 m in outdoor scenarios. However, in indoor and
electromagnetically-challenged environments, such as urban
canyons or forests, the coverage provided by GNSSs is not
sufficient. Therefore, for such environments, several alternative
Alessio De Angelis is with the Department of Engineering, University of
Perugia, Italy. E-mail: alessio.deangelis@unipg.it. Carlo Fischione is with
ACCESS Linnaeus Centre, Department of Electrical Engineering, KTH Royal
Institute of Technology, Stockholm, Sweden. E-mail: carlofi@kth.se
positioning technologies have been developed, particularly for
wireless indoor localization; for an extensive review see [3]
and references therein.
The topic of localization using wireless signals has received
much attention in many different research areas. Numerous
solutions and applications have been proposed, including
distributed algorithms for wireless sensor networks [4], the
robust distributed method proposed in [5], and the technique
for passive device-free sensor localization proposed in [6]. The
key aspect of cooperation between the nodes of a wireless
network for localization purposes has been investigated in
[7], and experimentally evaluated in [8]. The fundamental
importance of position information for nodes in wireless
sensor networks is highlighted in [9], where topology-aware
estimation algorithms are developed, using techniques from
the spectral graph theory field. In this context, distributed
gossip algorithms for sensor networks localization have been
studied in [10] and [11]. In particular, [10] shows that kernel
averaging techniques for localizing an isotropic source based
on measurements by distributed sensors may provide im-
proved robustness and accuracy than least squares estimators.
Moreover, in [11], a distributed algorithm is proposed and
characterized, which estimates the positions of nodes given
a minimal number of anchors and using only data-exchange
between neighboring nodes.
The robustness with respect to range-measurement outliers
due to non-line of sight conditions has been analyzed in [12],
which introduces a robust extended Kalman filter for locating
a mobile terminal node in wireless networks. Moreover, in
[13], a Bayesian algorithm is proposed for self-localization
and tracking of a mobile node through ranging with known-
position anchors. The algorithm is robust to NLOS propagation
by combining a Markov model for sight-state estimation
and a particle filter for location estimation, with a simple
general motion model. Experiments using the IEEE 802.15.4a
chirp-spread-spectrum ranging technology show accuracy of
approximately 2 m in an indoor environment with varying
NLOS and LOS conditions. In [14], the concept of schedule-
based network localization is introduced, which enables self-
localization of mobile nodes, as well as localization of the
entire network, without communication overhead.
In the context of wireless short range indoor positioning,
Pulse-based Ultra-Wideband (UWB) techniques are of particu-
lar interest due to numerous good qualities, such as a fine time
resolution (which allows for a ranging accuracy of the order
of centimeters when used in conjunction with time-of-arrival
2
methods) a resilience to multipath propagation effects of the
wireless channel, and a low-power device implementation [15].
Recent results on the fundamental limits for wideband radio
localization of static nodes have been presented in [16]. For
the case of mobile cellular systems, a performance analysis
has been presented in [17]. The fundamental limitations of
mobile localization have been extensively studied in [18] for
bearings-only measurements in the target tracking field, and
in [19] in the terrain-aided navigation field. In the context
of sensor networks node localization, the Cramér-Rao lower
bound (CRLB) has been evaluated for the case of received-
signal-strength (RSS) range measurements in [20] and in [21].
Therein, the RSS model has been linearized with only ranging
measurements and without using any velocity and orientation
information.
To improve the performance of a positioning system in
terms of reliability of the estimates (integrity), accuracy, and
availability, it is appealing to process information obtained
from a number of sensors by means of fusion techniques [22],
[23]. These techniques involve the processing of different
information sources to overcome the fundamental limitations
of sensor measurements such as GNSSs, inertial navigation
systems (INSs), odometry, and local radio technologies. An
extensive survey of the most common information sources and
sensor fusion approaches is provided in [24] in the context
of automotive positioning. Concerning information fusion for
UWB, the commonly employed techniques are based on the
extended Kalman filter (EKF), [25], or other non-linear filter-
ing approaches such as particle filters methods [26]. In [27],
a 6 degrees-of-freedom tracking system is presented, which
performs sensor fusion on a UWB distance-measuring device
and an inertial sensor consisting of triaxial accelerometers and
gyroscopes. In [28], an UWB sensor fusion technique based
on the complementary filter approach is presented, where the
error states are estimated by an EKF using the UWB ranging
measurements. Subsequently, the estimated errors are fed back
to correct for the inertial navigation system biases, which
would otherwise grow unbounded.
In this paper, we investigate a novel sensor fusion technique,
based on Pareto optimization, for node self-localization using
heterogeneous information sources. It is assumed that a node
wishing to estimate its own position has available ranging,
speed and orientation information. We stress here that the
problem that we consider in this paper is self-localization,
which is fundamentally different from the target tracking
problem that is widely studied in the signal processing litera-
ture, see e.g. [29]. Furthermore, the fundamental performance
limitation of the linear fusion is investigated in this paper.
This investigation is a substantial extension of our earlier
work in [30] and [31], where the performance of the Pareto
optimization were not fully investigated and no derivation of
the CRLB were considered. The posterior CRLB [32] and the
parametric CRLB [19] are characterized for both any generic
trajectory of the mobile node and for a specific trajectory.
Since ranging information gives an error with low bias and
a high variance, and measurements of the speed and absolute
orientation (dead-reckoning) give an error with a low variance
and a high bias, the overall information is jointly processed
to overcome the individual limitations of ranging, speed and
orientation measurements. By numerical simulations and by
evaluating the CRLB, we show that our new method may
outperform existing solutions, including several commonly
employed techniques based on the Kalman filter. However,
the price payed is a higher computational complexity.
This paper has several distinguished features compared to
the related work mentioned above. Specifically, our system
model is original because the mobile node position estima-
tion is performed by ranging measurements, with respect to
fixed position nodes, together with velocity and orientation
measurements without any linearization or simplification. For
the derivation of the estimator, we assume a very simple
dynamical model of the movement of the mobile node. This
model allows for flexibility, robustness, low-complexity imple-
mentation and works well in several different scenarios, as we
show later. Only local processing at such a node is used, thus
enabling a completely distributed strategy. Furthermore, the
estimation is performed in a cooperative fashion, in the sense
that operation of the estimator at the mobile node relies on
the cooperation of fixed-position nodes acting as responders
during ranging measurements.
Therefore, our method is of a general nature and can
be applied to any motion scenario. Compared to, e.g., [20]
and [21], we investigate the CRLB without any linearization
of the ranging model and we consider simultaneously ranging,
velocity and orientation measurements, which were not con-
sidered therein. Moreover, in this paper, we address one of the
main challenges for developing sensor fusion methods, namely
the statistical characterization of the estimation error. Based
on this characterization, a sensor fusion method is derived
by solving a Pareto optimization problem, where a tradeoff
between the variance of the estimation error and its bias is
exploited. The Pareto optimization approach was first proposed
in our initial study in this context in [31].
The remainder of this paper is organized as follows: The
problem is formulated in Section II. Then, the proposed sensor
fusion method is derived in Section III. Furthermore, the
fundamental limits of an estimator that fuses ranging, velocity
and orientation measurements are investigated in Section IV
by means of the CRLB. Numerical simulation results are
presented in Section V. Finally, conclusions are drawn in
Section VI.
A. Notation
We denote real n-dimensional vectors with lowercase bold-
face letters, such as x, and Rn×m matrices with uppercase
boldface letters, such as A. The superscript (·)T indicates
the transpose of a matrix. Given two matrices A and B,
the inequality A ⪯B denotes that the matrix B −A is
positive semidefinite. Given a scalar function f(x) : Rn →R,
∇xf(x)
=
h
df(x)
dx1 , . . . , df(x)
dxn
iT
. Given a vector function
f(x) : Rn →Rn, we use the gradient matrix definition
∇xf(x) =
∇f1(x)
. . .
∇fn(x)
, which is the transpose
3
Slave 3
Slave M
0
Slave 2
Slave 1
Master
x1
x2
ݏଵ,ଵ, ݏଵ,ଶ 
ݏଶ,ଵ, ݏଶ,ଶ 
ݏଷ,ଵ, ݏଷ,ଶ 
ݏெ,ଵ, ݏெ,ଶ 
ݔଵ,௞, ݔଶ,௞ 
ܸ௞ 
߶௞ 
ݎଵ,௞ 
ݎଶ,௞ 
ݎଷ,௞ 
ݎெ,௞ 
Fig. 1.
Model of the localization system. The slave nodes are placed at fixed
and known positions. At time k, the unknown-position master node measures
its distance ri,k with respect to each slave i by measuring the round-trip-time
of UWB pulses. Speed Vk and orientation φk of the master are measured using
on-board sensors.
of the Jacobian matrix. We indicate with ˜x a noisy mea-
surement of x and with ˆx an estimator of x. The statistical
expectation operator is denoted by E {·}. p(x) denotes a
probability density function (pdf) of the random vector x.
II. PROBLEM FORMULATION
Consider a system where a mobile node, which we denote
as master, needs to estimate its own planar position xk ∈R2,
with xk = [x1,k
x2,k]T , where k is the discrete time index.
To do so, the master measures its distance with respect to
M devices, which we denote as slaves, with M ≥4. We
do not assume any a-priori information on the mobility of
the master. It is allowed to move by following a linear
trajectory with constant velocity as well as a random trajectory
with random acceleration. An example of such a system is
provided by the in-house developed experimental platform
that has been characterized in [28], [33], where the distance
measurement is obtained by means of the round-trip-time
of UWB signals. This distance measurement approach, also
known in the literature as two-way time-of-arrival, avoids the
need for synchronization between mobile node and slaves,
therefore allowing for asynchronous operation and reduced
infrastructure requirements. Alternative approaches, such as
one-way time-of-arrival and time-difference-of-arrival, suffer
from clock skew and require additional infrastructure for
precise clock synchronization [34].
We assume that the master is equipped with UWB trans-
mission and sensors that measure its speed and absolute
orientation. As an example, speed sensors might be imple-
mented by means of odometry in vehicular or mobile-robotics
applications, such as those modeled in [35]. Furthermore,
information about absolute orientation could be obtained from
a heading sensor, for example a compass or magnetometer,
as in [36]. The diagram in Fig. 1 shows a representation of
the system model. Furthermore, we assume that the slaves
are located in fixed and known positions that are denoted as
si,j, with i = 1, . . . , M and j = 1, 2. We denote the ranges,
namely the distances between the master and each slave at
time instant k, by ri,k, i = 1, . . . , M. The measurements of
the ranges, denoted as ˜ri,k, are affected by additive zero-mean
Gaussian noise wr,i,k:
˜ri,k = ri,k + wr,i,k .
(1)
The variance of wr,i,k depends on the range, according to the
following exponential model that has been derived from real-
world round-trip-time UWB measurements in [28]:
σ2
r,i,k = σ2
0 exp(κσri,k) ,
(2)
where σ2
0 and κσ are known constants.
Further, we denote the speed and orientation of the master
by Vk and φk, respectively, whereas ˜Vk is the measurement of
the speed and ˜φk the measurement of the orientation. These
measurements are expressed as
˜Vk = Vk + wV,k ,
˜φk = φk + wφ,k
where the noise terms wV,k and wφ,k are zero-mean Gaussian
random variables with known variances σ2
V and σ2
φ, respec-
tively.
In this paper, assuming the measurement models described
above, we propose a novel method to combine and improve
existing estimators which are well known and commonly used.
Specifically, we use a so called loosely-coupled approach,
in which we fuse estimates already available from rang-
ing measurements, with the Weighted Least Squares (WLS)
method, and from velocity and orientation measurements, by
performing dead-reckoning. Note that, in order to achieve the
best possible accuracy, a tightly-coupled approach is necessary,
where information from the ranging, velocity and orientation
sensors are fused directly without preprocessing. However, the
loosely-coupled approach is motivated by the reduced compu-
tational complexity as well as by the increased robustness and
modularity that it offers.
This paper has two goals: First, to propose an highly
accurate estimator of the position of the master by processing
the results provided by existing ranging-only WLS and dead
reckoning estimators. Second, to understand the fundamental
limitations when estimating the position of the master by the
fusion of information from these estimators. In the remainder
of this section, we start by formulating the estimation problem
that leads to the derivation of the proposed estimator.
The estimate of the master’s position at time k +1 by using
measurements available at that time is denoted by ˆxk+1|k+1 =

ˆx1,k+1|k+1
ˆx2,k+1|k+1
T . We propose that it is derived as
follows:
ˆx1,k+1|k+1 = α1,kˆx(r)
1,k+1 + β1,kˆx(v)
1,k+1 ,
(3)
ˆx2,k+1|k+1 = α2,kˆx(r)
2,k+1 + β2,kˆx(v)
2,k+1 ,
(4)
where ˆx(r)
k+1 =
h
ˆx(r)
1,k+1
ˆx(r)
2,k+1
iT
is the position estimate
based only on the ranging measurements, and ˆx(v)
k+1
=
h
ˆx(v)
1,k+1
ˆx(v)
2,k+1
iT
is the position estimate based on the dead
reckoning block, which processes the position estimate at the
previous time step and the on-board speed and orientation
sensors. The terms α1,k, α2,k, β1,k and β2,k are the sensor
fusion design parameters that need to be optimally chosen.
4
Our estimation problem is separable on the x1 and x2 axes, in
the sense that the estimates on the x2 axis does not affect the
x1 axis, and vice versa. Therefore in the following we will
provide derivations only for the x1 component, because the
derivations for the x2 component are similar to those for the
x1 component.
Typically, the dead reckoning block provides us with an
estimate having the following expression:
ˆx(v)
1,k+1 = ˆx1,k|k + ˜VkT cos ˜φk .
(5)
Note that we do not make the assumption of a linear motion
model for the dead-reckoning block. Furthermore, the dead-
reckoning estimate is biased because the orientation measure-
ment appears as the argument of a cosinus. It follows that
the estimate (3) is biased. This motivates us to model our
estimator (3) as
ˆx1,k+1|k+1 ≜x1,k+1 + wx1,k+1 ,
(6)
where wx1,k+1 is the error in the position estimate. This error
has a non zero average, since the dead reckoning block gives
a biased estimate.
It is possible to write a simplified form of Equations (3),
by first noting that
E

ˆx1,k+1|k+1
	
= α1,kx1,k+1 + β1,kx1,k+1
+ α1,k E {wr,k+1}
+ β1,k E {wv,k+1} + β1,k E {wx,k}
= x1,k+1 + E {wx,k+1}
where wr,k+1 and wv,k+1 are the errors in the estimate
obtained from the ranging block and from the dead-reckoning
block, respectively. Now, in order to have α1,kx1,k+1 +
β1,kx1,k+1 = x1,k+1, it must be α1,k +β1,k = 1 and therefore
α1,k = 1 −β1,k .
(7)
Hence, by substituting (7) in (3), we obtain
ˆx1,k+1|k+1 = (1 −β1,k) ˆx(r)
1,k+1 + β1,kˆx(v)
1,k+1 ,
(8)
In the following, we will use the simplified expression (8)
instead of (3).
From (6) it follows that the error is
wx1,k+1 = ˆx1,k+1|k+1 −x1,k+1
(9)
= (1 −β1,k) w(r)
x1,k+1 + β1,kwx1,k
+ β1,kT ˜Vk cos ˜φk −β1,kT Vk cos φk ,
(10)
where w(r)
x1,k+1 is the error of the estimates obtained from
the ranging. To develop an estimator for the master’s po-
sition, we define a cost function that takes into account
the estimator variance and bias simultaneously. We de-
fine the bias term of the estimation error as µw1,k+1 ≜
E

wx1,k+1
	
, and the variance term of the estimation er-
ror as σ2
w1,k+1 ≜E
nˆx1,k+1|k+1 −E

ˆx1,k+1|k+1
	
2o
=
E
nwx1,k+1 −E

wx1,k+1
	
2o
.
Dead 
Reckoning
Least 
Squares
Information 
Fusion
Speed 
sensor
Absolute 
orientation 
sensor
Range 
measurement 1
Range 
measurement M
)
(
1
ˆ r
k
x
)
(
1
ˆ v
k
x
1
1
ˆ

 k
k
x
k
k
xˆ
k
V~
kI~
1
,1~

k
r
1
,
~

k
M
r
Fig. 2. Sensor fusion system architecture used at a master node for estimating
its position at time k + 1 given available information. The dead reckoning
block gives speed and absolute orientation information, whereas the least
square block gives ranging information as computed with respect to slave
nodes. The information provided by these two blocks is then fused to provide
an accurate estimation ˆxk+1|k+1.
We are now ready to formulate a Pareto optimization
problem to select the fusion coefficients at each time k:
min
β1,k
ρ1,kµ2
w1,k + (1 −ρ1,k) σ2
w1,k
(11)
s.t.
β1,k ∈B = [−1, 1] ,
ρ1,k ∈R = [0, 1] .
This problem is motivated by that we would like to reduce
as much as possible both the average of the estimation error,
namely the bias due to the dead reckoning block, and the
variance of the estimation error. The Pareto weighting factor
(or scalarization coefficient) ρ1,k must be chosen for each time
k so to trade off the low bias and high variance of the error of
the ranging estimate with the high bias and low variance of the
error of the dead-reckoning estimate [37, Section 4.7.5]. Based
on the statistical properties of the bias, which we characterize
below, the bias itself may grow unstable with time, which
motivates the constraint β1,k ∈B. The intuition is that the
minimization of both the average and the variance of the
estimation error at every time k does not ensure that the bias is
stable or decreases over time. Notice that in the special case of
ρ1,k = 0 ∀k, the cost function is the variance of the estimation
error, and when ρ1,k = 0.5 ∀k the cost function is the mean
square error (MSE).
The challenge for such an optimization is the analytical
characterization of the cost function (11), which we study in
the following section.
III. DERIVATION OF A SENSOR FUSION ESTIMATOR
The architecture of the proposed estimator is shown in
Fig. 2. In the following subsections, we provide a statistical
characterization of the ranging estimate ˆx(r)
k+1 and the dead-
reckoning estimate ˆx(v)
k+1. Then, in subsection III-D we solve
the optimization problem, thus achieving the optimal weight-
ing coefficients β1,k and β2,k to use in Eq. (8). The original
contribution of this section is the statistical characterization of
the estimation errors, and the optimal weighting derivation in
Proposition 3.7.
5
A. Ranging estimate statistical characterization
In this subsection we give a new statistical characterization
of the estimation error of the ranging estimate. We consider the
problem of finding the first and second order moments of the
ranging estimates of the x1 and x2 components of the position.
This requires first a derivation of the ranging estimates. In the
following, for the sake of simple notation, we drop the k + 1
indices.
Since it must be that M ≥4 for linear-based ranging
estimates, the master’s position estimate ˆx(r) may be obtained
by following a commonly employed strategy (see, e.g., [38])
as the solution of Aˆx(r) = ˜s, where A and ˜s are easily seen
to be
A =


2 (s1,1 −sM,1)
2 (s1,2 −sM,2)
...
...
2 (sM−1,1 −sM,1)
2 (sM−1,2 −sM,2)

,
(12)
˜s =
 r2
M −r2
1, . . . , r2
M −r2
M−1
T + a + ws ,
(13)
where the constant vector a is given by the following function
of the known coordinates of the slaves
a =


s2
1,1 −s2
M,1 + s2
1,2 −s2
M,2
. . .
s2
M−1,1 −s2
M,1 + s2
M−1,2 −s2
M,2

,
and ws is a random vector given by
ws ≜


w2
rM + 2rMwrM −w2
r1 −2r1wr1
. . .
w2
rM + 2rMwrM −w2
rM−1 −2rM−1wrM−1

.
(14)
We derive the ranging estimate by the weighted least squares
(WLS) technique, since ranging measurements are affected by
unequal statistical errors across the various links. Then, the
solution to the following classical WLS optimization prob-
lem gives the ranging estimates: minˆx(r)



W
1
2 ˜s −Aˆx(r)



,
where W is a positive-definite weighting matrix introduced
for the purpose of emphasizing the contribution of those range
measurements that are deemed to be more reliable [39], and
∥·∥is the Euclidean norm. A common choice for W is R−1,
namely the inverse of the noise covariance matrix:
R ≜E
n
(ws −E {ws}) (ws −E {ws})T o
.
(15)
Such a choice is motivated by that it would give a maximum
likelihood estimator in the case of Gaussian errors. Thus the
well-known solution to the WLS problem is
ˆx(r) =
AT WA

−1 AT W˜s ,
(16)
where W = R−1.
We are now ready to introduce the new part of this subsec-
tion. From (16) we can characterize the first and second order
moments of the error in the ranging estimate. First, however,
we need the expression of the inverse of the error covariance
matrix. This inverse matrix has a simple expression and is
provided by the following lemma:
Lemma 3.1: Consider the matrix R as defined in (15). The
matrix has full rank, and the inverse is
R−1 =

I −
1
1 + q G11T

D−1 ,
where I is the (M −1) × (M −1) identity matrix, D is a
(M −1) × (M −1) diagonal matrix whose entries are Dll =
4r2
l σ2
wrl +2σ4
wrl, l = 1, . . . , M −1, G = pD−1 is a (M −1)×
(M −1) matrix with p = 4r2
Mσ2
wrM +2σ4
wrM , q =
M−1
P
i=1
p/Dii,
and 1 is the all ones vector.
Proof: The lemma follows from tedious algebraic com-
putations, the matrix inversion lemma, and the Woodbury
identity [39].
Based on Eq. (16) and the previous lemma, we are now in
the position of deriving the expectation and the variance of
the ranging estimation error wr, defined as wr ≜ˆx(r) −x =
h
w(r)
x1
w(r)
x2
iT
.
Lemma 3.2: The expectation of the ranging estimation
error is given by:
E {wr} =
AT WA

−1 AT W
×

1σ2
wrM −
h
σ2
wr1 , . . . , σ2
wrM−1
iT 
.
Proof: See Appendix A.
Furthermore, we have the following result:
Lemma 3.3: The correlation of the range estimation error
is given by
E

w2
r
	
=
AT WA

−1AT WCWT A
AT WA

−1 , (17)
where C is a matrix whose diagonal elements are Cll =
3σ4
wrM + 4r2
Mσ2
wrM + 3σ4
wrl + 4r2
l σ2
wrl −2σ2
wrM σ2
wrl , with
l = 1, . . . , M, and the off-diagonal elements are Clj =
3σ4
wrM −σ2
wrM σ2
wrj + 4r2
Mσ2
wrM −σ2
wrl σ2
wrM + σ2
wrl σ2
wrj ,
with l ̸= j.
Proof: See Appendix B.
This concludes the efforts to characterize the statistics of
the ranging estimation error. In the following subsection, we
focus on the statistics of the dead-reckoning estimation.
B. Dead-reckoning estimate statistical characterization
In this subsection, we give a new characterization of the
estimation errors of the dead-reckoning block. We have the
following results for the first and second order moments of
the dead-reckoning estimate:
Lemma 3.4: Let ˜φ(k) be Gaussian, and assume that ˜V (k)
and ˜φ(k) are statistically independent. Then
E
n
˜Vk cos

˜φk
o
= Vk cos (φk) e−
σ2
φ
2 ,
(18)
E
n
˜V 2
k cos2 
˜φk
o
= σ2
V
1
2 + 1
2 cos (2φk) e−2σ2
φ

. (19)
Proof: See Appendix C.
In the following subsection, we use these results to character-
ize the estimation bias.
6
C. Estimation Bias
We can now put together the statistical characterization
of the ranging and dead-reckoning estimations. We have the
following results that give the terms µw1,k and σ2
w1,k in
Eq. (11):
Lemma 3.5: Consider the estimation error (9). The bias is
µw1,k+1 = E {wx1,k+1} = (1 −β1,k) E
n
w(r)
x1,k+1
o
+ β1,k E {wx1,k} + β1,kT Vk cos (φk)

e−
σ2
φ
2 −1

.
(20)
Proof: See Appendix D.
Remark From the previous lemma, we see that the average
of the estimation error at time k + 1 depends on the bias of
time k. To avoid an accumulation of the bias, we need to
impose a condition on β1,k. Thus, the absolute value of the
average of the estimation error must be non-expansive, which
can be easily achieved when |β1,k|∈[0, 1], and thus we have
B = [−1, 1] in the Pareto optimization problem (11).
Lemma 3.6: Consider the estimation error (9). Then, the
second order moment is
E
n
w2
x1,k+1
o
= β2
1,kak + 2β1,kbk + E

w(r)
x1,k+1
2
(21)
where
ak = E

w(r)
x1,k+1
2
+ E
n
w2
x1,k
o
+ T 2 E { ˜V 2
k cos2(˜φk)}
+ T 2V 2
k cos2(φk)

1 −2e−
σ2
φ
2

−2 E {wx1,k} E {w(r)
x1,k+1}
−2 E {w(r)
x1,k+1}T Vk cos (φk)

e−
σ2
φ
2 −1

+ 2 E {wx1,k}T Vk cos (φk)

e−
σ2
φ
2 −1

,
(22)
bk = −E

w(r)
x1,k+1
2
+ E {wx1,k} E {w(r)
x1,k+1}
+ E {w(r)
x1,k+1}T Vk cos (φk)

e−
σ2
φ
2 −1

.
(23)
Proof: The result follows from Lemma 3.2, Eq. (17), and
Proposition 3.4.
Finally, to derive an expression for the variance of the
estimation error σ2
w1,k that we need in (11), we define the
variables vr, vx1 and vv as
vr ≜w(r)
x1,k+1 −E
n
w(r)
x1,k+1
o
, vx1 ≜wx1,k −E

wx1,k
	
,
vv ≜T ˜Vk cos

˜φk

−T Vk cos (φk) e−
σ2
φ
2 .
Notice that vr, vx1 and vv are independent; we denote their
variances as σ2
vr, σ2
vx1 and σ2
vv respectively. Also, note that
E {vr} = E {vx1} = 0. By substituting the expression for
the error bias (20) in Eq. (21) we obtain
σ2
w1,k+1 = E
n
(1 −β1,k) w(r)
x1,k+1 + β1,kwx1,k
+β1,kT ˜Vk cos

˜φk

−β1,kT Vk cos (φk)
−(1 −β1,k) E
n
w(r)
x1,k+1
o
−β1,k E

wx1,k
	
+β1,kT Vk cos (φk) e−
σ2
φ
2
2)
= E
n
(1 −β1,k)2 v2
r + β2
1,kv2
x1 + β2
1,kv2
v
+2 (1 −β1,k) β1,kvrvx1 + 2 (1 −β1,k) β1,kvrvv
+2β2
1,kvx1vv
	
= (1 −β1,k)2 σ2
r + β2
1,kσ2
x + β2
1,kσ2
v ,
(24)
where we have used that vr, vx1 and vv are independent and
E {vr} = E {vx1} = 0.
D. Pareto Optimization
In this section, we put together the results of the previous
sections to derive a position estimator by solving the opti-
mization problem in (11). Lemmata 3.5 and 3.6 give us the
analytical expression of the cost function of problem (11).
Thus, we have the following result, which is one of the core
contributions of this paper:
Proposition 3.7: Consider optimization problem (11). Let
the mean and the variance of the estimation error (9) be
computed by Lemmata 3.5 and 3.6. Then, the optimal solution
for a fixed Pareto weighting factor ρ1,k is
β∗
1,k(ρ1,k) = max (−1, min (ξ, 1)) ,
(25)
with
ξ = 2 (1 −ρ1,k) σ2
vr −2ρ1,kγ E {wr,k+1}
2 (1 −ρ1,k) η + 2ρ1,kγ2
,
where
η ≜σ2
vr + σ2
vx + σ2
vv ,
γ ≜−E {wr,k+1} + E {wx,k} + T Vk cos (φk)

e−
σ2
φ
2 −1

.
Proof: See Appendix E.
As a particular case of the previous proposition is given by the
MSE, when the optimal solution of problem (11) is obtained
with the Pareto weighting factor ρ1,k = 0.5 ∀k and results
β∗
1,k = max

−1, min

−bk
ak
, 1

,
(26)
where ak and bk are given by Eq. (22) and Eq. (23), respec-
tively. In general, the optimal value of β1,k in (25) depends
on the scalarization coefficient ρ1,k (see [37, Section 4.7]).
The best value of such a coefficient is found by building the
Pareto trade-off curve and selecting the “knee-point” on this
curve [37]. Thus, we choose ρ∗
1,k such that µw1,k and σ2
w1,k
computed in β∗
1,k(ρ∗
1,k) are σ2
w1,k ≃µ2
w1,k. This is given by
the solution to the following further optimization problem
min
ρ1,k
σ2
w1,k
β∗
1,k (ρ1,k)


−µ2
w1,k
β∗
1,k (ρ1,k)


2 ,
(27)
7
where we have evidenced that µw1,k and σ2
w1,k are computed
in β∗
1,k (ρ1,k), according to Eqs. (20) and (24). This problem
is highly non-linear, but simple and standard numerical pro-
cedures, based on a discrimination of ρ1,k, can be employed
to compute approximately the optimal Pareto coefficient ρ∗
1,k.
The standard numerical procedure is described in [37, Section
4.7.5]. Such a procedure consists of first defining a set of
values that ρ1,k can assume, then for each such value the
procedure does the following steps:
1) compute β∗(ρ1,k) via (25) for the current value of ρ1,k,
2) use the value of β∗(ρ1,k) computed in 1) and plug it
in (20) and in (24) to evaluate µ2
w1,k

β∗
1,k (ρ1,k)

and
σ2
w1,k

β∗
1,k (ρ1,k)

;
3) use
the
values
µ2
w1,k

β∗
1,k (ρ1,k)

and
σ2
w1,k

β∗
1,k (ρ1,k)

computed at step 2) to compute
σ2
w1,k
β∗
1,k (ρ1,k)


−µ2
w1,k
β∗
1,k (ρ1,k)


2 ,
(28)
and store such a value.
The approximate optimum ρ∗
1,k is the one which corresponds
to the minimum of the stored values of (28). Such a ρ1,k cor-
responds to the "knee point" of the Pareto trade off curve [37,
Section 4.7.5].
Finally, recall that the values of the optimal β2,k and
ρx2,k yielding the x2 component of the position estimate are
obtained similarly to what done for the x1 component.
IV. CRAMÉR-RAO LOWER BOUND
In this section, we are interested in investigating the funda-
mental performance limitation of our estimator proposed in the
previous section. To do so, given any estimation problem, the
CRLB is a fundamental tool for performance analysis, since it
quantifies the best achievable mean-square error performance
of estimators [39]. We wish to derive such a bound to
compare it to the performance of our estimator. There are two
fundamental approaches: the parametric CRLB and posterior
CRLB.
The parametric CRLB, that we denote as ParCRLB, is the
CRLB evaluated for a given state-space trajectory [19]. There-
fore it is equivalent to the problem of estimating unknown
deterministic states considering a zero process noise in the
dynamical system model; for the proof see, e.g., the derivations
in [40]. In other words, the state is treated as a deterministic
but unknown parameter for the purpose of the evaluation of the
ParCRLB. Therefore, the ParCRLB is simple to evaluate in a
numerical simulation environment where the true trajectory is
known, because the equations required for its calculation do
not contain any statistical expectation operator; see, e. g., [41].
However, due to its nature, the ParCRLB does not provide
information about the fundamental performance limitations in
the general case. For such a case, we need to investigate the
posterior CRLB, which we do in the following.
A. Posterior CRLB
For recursive estimation problems in the Bayesian frame-
work, where the state is treated as a random variable, the
bound is denoted as the posterior CRLB (PCRLB). A recursive
formulation of the PCRLB has been first introduced in [32].
In the following, we derive such a bound for our localization
problem. The core contribution of this section is summarized
in Proposition 4.5 below.
Differently from the case of the ParCRLB, where one is
interested in calculating the performance bound for a given
trajectory, here we aim at characterizing the fundamental
performance limits, on the average, for all the trajectories
belonging to a class, which is defined by a dynamical model.
Therefore we introduce the state pk ∈R4 to model the bi-
dimensional position of the master, its absolute speed and
its orientation, that is pk = [x1,k
x2,k
Vk
φk]T , and we
consider the following non-linear state-space model:
pk+1 = f (pk) + vk
(29)
yk = h (pk) + ek ,
(30)
where
vk
=
[v1,k
v2,k
v3,k
v4,k]T
∈
R4
∼
N (0, Q) is the white Gaussian process noise with Q =
diag
σ2
1, σ2
2, σ2
3, σ2
4


, where it is natural to assume that
σ1 = σ2, namely that the process noise on the two coordinates
has same variance, and ek ∈RM+2 ∼N (0, Rk) is the white
Gaussian measurement noise, which is independent of vk,
with Rk = diag

σ2
r,1,k, . . . , σ2
r,M,k, σ2
V , σ2
φ

. Simple algebra
gives that the state update and measurement functions are,
respectively,
f (pk) =


x1,k + T Vk cos φk
x2,k + T Vk sin φk
Vk
φk

,
h (pk) =
 h1, . . . , hM, hM+1, hM+2
T
=


p
(s1,1 −x1,k)2 + (s1,2 −x2,k)2
...
p
(sM,1 −x1,k)2 + (sM,2 −x2,k)2
Vk
φk


,
where we recall that T is the sampling interval.
Consider the dynamical system of (29) and (30), and denote
as Pk the covariance matrix of any unbiased estimator of the
system’s state at time k. Then, Pk is lower bounded by the
inverse of the Fisher information matrix Jk, that is J−1
k
⪯Pk.
Such a matrix can be recursively computed as [32]
Jk+1 =D22
k −D21
k
Jk + D11
k

−1 D12
k
(31)
with
D11
k = −E
n
∇xk [∇xk log p (xk+1|xk)]T o
,
(32)
D21
k =
D12
k

T = −E
n
∇xk

∇xk+1 log p (xk+1|xk)
T o
,
(33)
D22
k = −E
n
∇xk+1

∇xk+1 log p (xk+1|xk)
T o
−E
n
∇xk+1

∇xk+1 log p (yk+1|xk)
T o
,
(34)
where the initial condition is given by the information matrix
J0, which is computed based on the initial prior density p (x0).
8
The system model in (29) and (30) has additive Gaussian
noise. Hence, the terms in (32), (33) and (34) may be ex-
pressed as (see [32, eq. (34) – (36)])
D11
k = E
n
[∇pkf(pk)] Q−1 [∇pkf(pk)]T o
,
(35)
D12
k = −E {[∇pkf(pk)]} Q−1 ,
(36)
D22
k = Q−1 + E
n
[∇pkh(pk)] Rk
−1 [∇pkh(pk)]T o
.
(37)
In general, there is no guarantee that the expectations in (35)-
(37) exist. In other words, the CRLB may not exist. Therefore,
in the following we establish the existence of the CRLB
by deriving closed form expressions for the first two terms,
namely D11
k
and D12
k , while we derive bounds for the third
term D22
k . We start by proving the following useful Lemma:
Lemma 4.1: Consider the dynamical system in (29) – (30),
denote the speed and the orientation at time k = 0 as
V0 and φ0, respectively, and define the following quanti-
ties: c0 ≜cos φ0,
s0 ≜sin φ0,
c′
0 ≜cos 2φ0,
s′
0 ≜
sin 2φ0,
εk ≜e−
(k−1)σ2
4
2
. Then
E {Vk} = V0 ,
E

V 2
k
	
= (k −1) σ2
3 + V 2
0
E {cos φk} = c0εk ,
E {sin φk} = s0εk
E {sin φk cos φk} = 1
2s′
0ε2
k ,
E

cos2 φk
	
= 1
2 + 1
2c′
0ε4
k
Proof: See Appendix F.
Then, we note that the gradient of f is easily evaluated as
∇pkf(pk) =


1
0
0
0
0
1
0
0
T cos φk
T sin φk
1
0
−T Vk sin φk
T Vk cos φk
0
1

. (38)
By substituting Eq. (38) in Eq.(35) and using Lemma 4.1 we
obtain Eq. (39) (see top of next page). Note that, when k tends
to infinity, this matrix tends to a diagonal matrix, because εk
tends to zero, but the right-lower element tends to infinity.
This implies that, as k grows, the most relevant contribution
to the Fisher information matrix in Eq. (31) is given by the
D22
k term. Therefore, in the limit case, only the measurements
influence the PCRLB, not the model.
Similarly, by substituting Eq. (38) in Eq.(36), we obtain:
D12
k = −


σ−2
1
0
0
0
0
σ−2
2
0
0
T σ−2
1 c0εk
T σ−2
2 s0εk
σ−2
3
0
−T σ−2
1 V0s0εk
T σ−2
2 V0c0εk
0
σ−2
4

.
We now turn our attention to the D22
k
term in (31). For
the sake of simple notation, we drop the k subscript in the
remainder of this section and we define the following quantity:
di,j ≜
xi −sj,i
q
(x1 −sj,1)2 + (x2 −sj,2)2 .
After some algebra, the gradient of h is
∇ph(p) ≜


d1,1
. . .
d1,M
0
0
d2,1
. . .
d2,M
0
0
0
. . .
0
1
0
0
. . .
0
0
1

.
(41)
Let us introduce the following definition for the matrix Π in
Eq. (40) (see top of next page). By substituting Eq. (41) in
Eq. (37), and using the above definitions, we obtain
D22 = Q−1 +
 Π
0
0
R−1
2

,
(42)
where 0 is the 2 × 2 zero matrix and R2 is defined as the
2 × 2 right-lower block of the measurement noise covariance
matrix, that is R2 ≜diag

σ2
V , σ2
φ

.
In Eq. (42), there are two main issues: first, one has to
show that the expectations in the submatrix Π exist, and
then one has to compute them. These expectations are taken
over non-linear functions of the state random variables, which
makes it hard to find a closed-form expression of the PCRLB.
This is a well known issue, as pointed out in [42], where
after arguing the difficulty of deriving a PCRLB in cartesian
coordinates, as we do, the PCRLB has been derived by
using logarithmic polar coordinates. However, we note that
the model adopted in [42] does not not include speed and
orientation measurements by on-board sensors, and that the
derivation of the PCRLB in a coordinate system does not give
insight on which is the best estimator on another coordinate
system [37], [41], which greatly limits the model of [42] for
our case. In the cartesian coordinate system one could resolve
the issue by a Monte Carlo approach to calculate numerically
the PCRLB. However, analytical expressions of the PCRLB
are of fundamental interest for many localization applications,
such as when trying to understand where to place the slaves
in order to maximize the information content given by the
Fisher information matrix [43], [44], or when planning the
path of a robot to minimize the uncertainty in its location.
This further motivates the derivation of the analytical results
that we present below.
We must first prove that the expectations in the submatrix
Π exist. In the following, we establish the existence of such
expectations by deriving them in closed-form, and we derive
closed-form expressions of upper and lower bounds for the
entries of Π for which the exact closed form expression is
too complex for practical numerical evaluations. Then, we
calculate them numerically in Section V. To compute the
expectations on the diagonal elements of Π, we establish the
following result.
Proposition 4.2: Let q and z be independent Gaussian
random variables with average and standard deviation µq, σq
and µz and σz, respectively, with σz = σq 1. Then,
E

q2
q2 + z2

=
1
1 + Υ
"
1 +
∞
X
k=1
(−1)k
µF,k
(1 + Υ)k
#
, (43)
where
Υ ≜µ2
z + σ2
z
µ2
q + σ2
q
"
1 +
∞
X
k=1
(−1)k
σ2k
q µk
µ2q + σ2q

k
#
,
and µk denotes the k-th central moment of the noncentral
chi-squared distribution of 1 degree of freedom, and µF,k
1Recall that these are the variances of the process noise on the coordinates,
see Eq.(29).
9
D11
k = E















σ−2
1
0
T
σ2
1 cos φk
−T
σ2
1 Vk sin φk
0
σ−2
2
T
σ2
2 sin φk
T
σ2
2 Vk cos φk
T
σ2
1 cos φk
T
σ2
2 sin φk
T 2
σ2
1 cos2 φk + T 2
σ2
2 sin2 φk + σ−2
3
T 2 
σ2
1−σ2
2
σ2
1σ2
2

Vk cos φk sin φk
−T
σ2
1 Vk sin φk
T
σ2
2 Vk cos φk
T 2 
σ2
1−σ2
2
σ2
1σ2
2

Vk cos φk sin φk
T 2
σ2
1 V 2
k sin2 φk + T 2
σ2
2 V 2
k cos2 φk + σ−2
4















=


σ−2
1
0
T c0
σ2
1 εk
−T V0s0
σ2
1
εk
0
σ−2
2
T s0
σ2
2 εk
T V0c0
σ2
2
εk
T c0
σ2
1 εk
T s0
σ2
2 εk
T 2 
σ2
1+σ2
2
2σ2
1σ2
2 + σ2
2−σ2
1
2σ2
1σ2
2 c′
0ε4
k

+ σ−2
3
T 2 
σ2
1−σ2
2
σ2
1σ2
2

V0s′
0
2 ε4
k
−T V0s0
σ2
1
εk
T V0c0
σ2
2
εk
T 2 
σ2
1−σ2
2
σ2
1σ2
2

V0s′
0
2 ε4
k
T 2 (k −1) σ2
3

σ2
1+σ2
2
2σ2
1σ2
2 + σ2
1−σ2
2
2σ2
1σ2
2 c′
0ε4
k

+ σ−2
4


. (39)
Π = E

σ−2
r1 d2
1,1 + . . . + σ−2
rM d2
1,M
σ−2
r1 d1,1d2,1 + . . . + σ−2
rM d1,Md2,M
σ−2
r1 d1,1d2,1 + . . . + σ−2
rM d1,Md2,M
σ−2
r1 d2
2,1 + . . . + σ−2
rM d2
2,M

.
(40)
denotes the k-th central moment of the doubly noncentral F-
distribution.
Proof: See Appendix G.
The previous proposition allows us to establish the exis-
tence of the expectations on the diagonal elements of Π. A
similar result can be derived for the off-diagonal elements.
Unfortunately, these closed form expressions are of limited
practical usage due to the computational complexity of the
many summations and function evaluations therein involved.
Accordingly, we derive the following upper and lower bounds,
which are much more easy to use in practice:
Proposition 4.3: Let q and z be independent Gaussian
random variables of average and standard deviation µq, σq
and µz and σz, respectively. Then
eα−ln(σ2
q+µ2
q+σ2
z+σ2
z) ≤E

q2
q2 + z2

≤1 ,
where α = −γe −ln (2)−2 ln (σq)−M
0, 1/2, −µ2
q/2σ2
q


,
the symbol γe denotes the Euler Gamma constant, and the
symbol M(a, b, z) denotes Kummer’s confluent hypergeomet-
ric function [45].
Proof: See Appendix H.
To compute bounds on the off-diagonal elements of Π, we
give the following proposition:
Proposition 4.4: Let q and z be independent Gaussian
random variables with average and standard deviation µq, σq,
and µz and σz, respectively. Let the expectation be computed
in the Cauchy principal value sense. Then
−1
2 ≤E

qz
q2 + z2

≤1
2 ,
Proof: The bounds follow immediately from that −1/2 ≤
qz/(q2 + z2) ≤1/2, ∀q ∈R, ∀z ∈R.
We are now in the position to compute the upper and
lower bounds on the Fisher information matrix, in the positive
semidefinite sense, by the following result.
Proposition 4.5: Let J(ub) and J(lb) be two matrices con-
sisting of element-wise upper and lower bounds, respectively,
of the Fisher information matrix J in Eq. (31) as obtained by
Propositions 4.3 and 4.4 for the elements of Π of Eq. (42).
Let
J(ub,G)
ii
≜









J(ub)
ii
,
if
J(ub)
ii
> min

u(row)
i
, u(col)
i

min

u(row)
i
, u(col)
i

+ ǫ ,
if
J(ub)
ii
≤min

u(row)
i
, u(col)
i

(44)
J(ub,G)
ij
≜J(ub)
ij
,
∀i, j = 1 . . . n , with i ̸= j ,
(45)
J(lb,G)
ii
≜J(lb)
ii
,
∀i = 1 . . . n ,
(46)
J(lb,G)
ij
≜









J(ub)
ij
if
J(lb)
ii
> min

l(row)
i
, l(col)
i


min

J(lb)
ii /l(row)
i
, J(lb)
ii /l(col)
i

−ǫ

J(lb)
ij
if
J(lb)
ii
≤min

l(row)
i
, l(col)
i

∀i, j = 1 . . . n , with i ̸= j ,
(47)
where u(row)
i
=
X
j̸=i
J(ub)
ij
 , u(col)
i
=
X
j̸=i
J(ub)
ji
 , l(row)
i
=
X
j̸=i
J(lb)
ij
 , l(col)
i
=
X
j̸=i
J(lb)
ji
 and ǫ is an arbitrary constant
with ǫ > 0. Then
0 ⪯J(lb,G) ⪯J ⪯J(ub,G) .
(48)
Proof: See Appendix I.
This concludes the steps to establish the existence and
upper and lower bounds of the Fisher information matrix.
In the following section, we use these results to compare
the performance of our proposed localization method to the
CRLB.
V. SIMULATION RESULTS
In this section we present simulation results to validate
both the PCRLB derivations and the estimator we presented
in the previous sections. In particular, we compare our new
estimator based on Pareto optimization that we have developed
in Section III to some solutions from the literature and to the
PCRLB we derived in Section IV.
10
A. Performance of the proposed sensor fusion method
In this section we present extensive simulation results of our
new estimator (8) based on the optimal fusion coefficient of
Eq. (25) and Pareto weighting factor of Eq. (27). The values
of the parameters used in the numerical simulations are the
following:
• Speed measurement noise standard deviation σv
=
0.05 m/s, to model the worst-case performance of an
odometry sensor, [35];
• Orientation measurement noise standard deviation σφ =
π/8 rad, to model the worst-case performance of a mag-
netometer subject to disturbances due to the environment,
see, e.g., [36], [46];
• Ranging noise model parameters (in Eq. (2)): σ0 = 0.25,
κσ = 0.25;
• Sampling time interval T = 0.1 s.
The following two scenarios have been considered:
(A) Linear trajectory with constant speed of 0.1 m/s. This is a
deterministic trajectory, that is, no process noise is added.
As an example, this scenario emulates the real-world
situation in which the mobile node is on a vehicle moving
along a straight line without performing any maneuver.
(B) Piece-wise-linear (PWL)-acceleration trajectory, which is
more complex and realistic than the linear one, with
a maximum acceleration of 0.5 m/s2, see Fig. 5. This
scenario emulates a mobile node performing maneuvers
with an acceleration range consistent with the indoor
environment, e.g. a walking person carrying the mobile
node.
The absolute value of β∗, provided by the algorithm,
has been clipped to 0.99 to ensure that the average of the
estimation error is non-expansive, as discussed in Section III.
Note that, in 3.7, in the expression for γ the true speed
Vk and orientation φk terms are present. Since these are not
available, we use these approximations:
Vapprox (k) ≃1
T
ˆx1,k|k −ˆx1,k−1|k−1

2
+
ˆx2,k|k −ˆx2,k−1|k−1

2 1
2 ,
(49)
φapprox (k) ≃arctan
 ˆx2,k|k −ˆx2,k−1|k−1
ˆx1,k|k −ˆx1,k−1|k−1

.
(50)
Also, note that the second-order moment of the ranging error,
given by (17), depends on the true ranges ri. Therefore, for
the calculation of the term E

w(r)
x1,k+1
2
in bk and ak,
we use the approximations for the ranges
r(approx)
i,k+1
=
q
(si,1 −ˆxk)2 + (si,2 −ˆyk)2 ,
i = 1, . . . , M .
(51)
The absolute error and empirical Cumulative Distribution
Function (CDF) obtained by our new method for scenario A
are shown in Fig. 3 and 4. In this scenario, a root-mean-
squared error (RMSE) of 4.0 cm in the position estimate
over the entire trajectory has been obtained by our method.
Furthermore, in Fig. 5 the more realistic PWL-acceleration
case is shown, where an RMSE of 5.5 cm is obtained.
0
10
20
30
40
50
60
0
0.1
0.2
0.3
0.4
0.5
0.6
x1 error [m]
Time [s]
 
 
ranging error
dead reckoning error
proposed sensor fusion error
0
10
20
30
40
50
60
0
0.2
0.4
0.6
x2 error [m]
Time [s]
 
 
ranging error
dead reckoning error
proposed sensor fusion error
Fig. 3.
Absolute errors in scenario A. The dead-reckoning estimates,
asterisks, have a high bias, while the WLS ranging estimates, dashed line,
present a higher variance, but a reduced bias. The proposed sensor fusion
technique is able to reach a good tradeoff between estimator bias and variance.
0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Absolute position error [m]
Empirical CDF
 
 
Ranging, RMSE = 0.113 m
Dead−reckoning, RMSE = 0.223 m
Proposed sensor fusion, RMSE = 0.035 m
Fig. 4. Empirical CDF of the absolute positioning error in 2D for scenario A.
Using the proposed method (blue line), about 95% of the position estimates
have an error smaller than 7 cm, and the performance is improved with respect
to both the ranging and the dead-reckoning stand-alone systems.
Numerical simulation results show that the approximations
(49)-(51) give good performance and that the proposed method
yields a good trade-off between variance and bias of the esti-
mated position. Therefore, the proposed estimator overcomes
the limitation of dead reckoning (that is the slowly accumu-
lating bias) while also reducing the relatively higher variance
of the ranging estimator. The result is an overall smooth
and accurate estimate of the trajectory. In the following, we
compare the proposed estimator to other solutions that are
commonly adopted in similar localization problems.
B. Performance Comparison
Here we compare the sensor fusion based on Pareto op-
timization proposed in this paper to several other methods
in the two mobility scenarios considered in the previous
section. To provide an extensive comparison, the simulations
11
0
1
2
3
4
5
0
1
2
3
4
5
x1 position [m]
x2 position [m]
 
 
Ranging position estimates
Dead reckoning
Estimated trajectory
True trajectory
Slaves
Fig. 5.
PWL-acceleration trajectory case, scenario B. The proposed sensor
fusion technique closely tracks the true trajectory with RMSE = 5.5 cm.
TABLE I
AVERAGE EXECUTION TIME PER ITERATION [s].
Matlab® implementation on an Intel® Core2 Duo,
2.40 GHz PC with 3 GB RAM.
EKF
1.3 × 10−4
MSE Eq. (26)
1.7 × 10−4
LC-KF
1.9 × 10−4
UKF
4.2 × 10−4
Proposed sensor fusion, Eq. (25)
4.5 × 10−4
are performed for different values of the sampling period T
and of the speed V for scenario A, whereas for scenario B the
simulations are performed for different values of T and of the
maximum acceleration. Moreover, for each configuration, the
simulations are performed for 10 realizations of the random
noises and the resulting RMSE values are averaged.
The sensor fusion methods available from the literature that
we have considered are the following:
• Extended Kalman Filter (EKF) [25], where the measured
speed and orientation data are employed as inputs, and the
ranges as measurements. No motion model is assumed.
• Unscented Kalman Filter (UKF) [47], with the same
state-space model as the EKF above.
• Loosely Coupled Kalman Filter (LC-KF), where range
measurements are preprocessed with the WLS algorithm
to obtain preliminary position measurements, thus using
a linear measurement equation. The usual linearization of
the state function is performed for calculating the input
noise matrix [25].
Notice that the methods based on the Kalman filter theory
have been implemented without assumptions about the motion
model, to provide a fair comparison with the developed method
of this paper. If prior information about motion models is
included in the state space formulation, the Kalman approaches
might be able to provide better performance, as we will
show in the example in the next section. Notice also that we
provide results of the MSE special case when ρ1,k = 0.5
∀k in Proposition 3.7. This method does not provide good
10
−2
10
−1
10
0
0
0.1
0.2
0.3
0.4
0.5
Speed [m/s]
RMSE [m]
 
 
MSE special case
Proposed sensor fusion
EKF
UKF
LC−KF
Fig. 6. RMSE vs speed of the different methods for Scenario A with T = 0.5 s.
10
−1
10
0
10
1
0
0.05
0.1
0.15
0.2
0.25
Maximum acceleration [m/s2]
RMSE [m]
 
 
MSE special case
Proposed sensor fusion
EKF
UKF
LC−KF
Fig. 7.
RMSE vs maximum acceleration of several methods for Scenario B
with T = 0.1 s.
performance compared to the more general case of choosing
the best ρ1,k by problem (27). However, it features a much
smaller computational complexity.
Comparisons for two selected simulation configurations are
provided in Fig. 6 and 7, as obtained in the two scenarios
with different values of T . The proposed sensor fusion based
on Pareto optimization outperforms the other methods, while
the KF-based methods have similar performance among them.
However, the proposed sensor fusion based on Pareto opti-
mization requires the highest computational load among the
implemented techniques, as shown in the average simulation
execution times listed in Table I, although it is close to the
execution time of UKF.
C. CRLB evaluation
In this section, we provide a comparison of the RMSE
performance of the proposed sensor fusion method based on
Pareto optimization and the Kalman filter based methods we
used in the previous section to the fundamental limits given
by the Cramér Rao lower bound. In particular, we first provide
comparisons with the ParCRLB, and then with the PCRLB as
derived in Section IV.
The performance of the filters was compared based on the
RMSE of the position estimate, where the error is defined
as the Euclidean distance between the true position and the
estimate. The results compared to the ParCRLB for several
specific trajectories are presented in Fig. 8, where the error
curves were obtained by averaging over 100 Monte-Carlo runs.
We notice that the proposed sensor fusion method based on
Pareto optimization provides better performance than the other
considered KF methods as it is the closest to the ParCRLB. We
12
 
 
ParCRLB
Proposed fusion
Special case (MSE)
LC−KF
EKF
UKF
0
1
2
3
4
5
6
7
8
9
0
0.02
0.04
0.06
0.08
0.1
0.12
Time [s]
RMSE [m]
(a)
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
0
0.05
0.1
0.15
Time [s]
RMSE [m]
(b)
0
10
20
30
40
50
0
0.02
0.04
0.06
0.08
0.1
0.12
Time [s]
RMSE [m]
(c)
0
5
10
15
20
25
0
0.05
0.1
0.15
Time [s]
RMSE [m]
(d)
Fig. 8.
Error performance of the proposed sensor fusion method based on Pareto optimization of Eq. (25) and other methods against the theoretical bound
given by the square root of the trace of the ParCRLB (dashed line): (a) linear trajectory, constant velocity of 0.5 [m/s]; (b) linear trajectory, constant velocity
of 1 [m/s]; (c) spline trajectory, maximum acceleration of 0.3 [m/s2]; (d) spline, maximum acceleration of 1 [m/s2].
13
0
20
40
60
80
100
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
trajectory step
RMSE [m]
 
 
PCRLB
Proposed sensor fusion
Special case (MSE)
KF loosely coupled
EKF CV
Fig. 9. Numerical simulation results on the class of trajectories defined by the
model in (29)-(30). The behavior of the proposed sensor fusion method based
on Pareto optimization is compared to other estimators from the literature and
to the PCRLB derived in Section IV. Error curves are obtained as the average
of 1000 Monte Carlo iterations.
stress that none of the considered methods makes assumptions
about the motion model, therefore in this sense the ParCRLB
is an overly optimistic bound.
Furthermore, the numerical simulation results related to the
model in (29)-(30) are shown in Fig. 9, in comparison with
the derived PCRLB in (31). The method denoted as “EKF
CV” is implemented using the motion model in (29), and the
same model is used to generate the trajectories. Specifically,
a standard EKF is implemented [25], by linearizing about
the current estimate of the mean and covariance, where the
Jacobian matrices are computed by numerical differentiation
at each iteration. In this case, we notice that EKF CV out-
performs the proposed sensor fusion method based on Pareto
optimization and approaches the PCRLB. This confirms that
the knowledge of the specific underlying motion model is
beneficial for estimation. However, the proposed sensor fusion
method still provides good performance, improving over the
LC-KF. Conversely to EKF CV, the proposed method does not
assume any specific motion model, thus it can be applied to a
wider class of trajectories, proving to be robust with respect to
model mismatch. Nevertheless, the same Pareto optimization
approach presented in this paper could be applied to derive an
estimator which exploits the knowledge of the motion model,
although such derivation is outside of the scope of this paper.
VI. CONCLUSION
In this paper, the fundamental limitations of position esti-
mation for a mobile node by the fusion of information from
ranging, velocity and angular measurements were investigated.
Upper and lower bounds of the posterior Cramér-Rao lower
bound were derived. A sensor fusion method based on Pareto
optimization was proposed. The analysis was validated by
Monte Carlo simulations, which showed that the proposed
method provides better performance than Kalman filtering-
based methods when no knowledge of the dynamic model of
the underlying system is assumed or when the model is time
varying.
Future work includes the extension of the proposed esti-
mator to cases when the mobility model is known. The good
performance achieved by the proposed sensor fusion method
based on Pareto optimization in the case of no motion model
assumption lets us think that the use of a motion model would
substantially farther improve the performance of the proposed
method.
APPENDIX A
PROOF OF LEMMA 3.2
By taking the mean of Eq. (16), we obtain:
E
n
ˆx(r)o
=
AT WA

−1AT W
×

1 E {˜r2
M} −E
n
˜r2
1, . . . , ˜r2
M−1
T o
+ a

=x +
AT WA

−1 AT W
×

1σ2
wrM −
h
σ2
wr1 , . . . , σ2
wrM−1
iT 
=x + E {wr} ,
whereby the lemma follows.
APPENDIX B
PROOF OF LEMMA 3.3
By taking the square of Eq. (1) we get ˜r2
i = r2
i + w2
ri +
2riwri . Then, using this expression in (14) and substituting in
(16), we have the following expression for the ranging-based
estimate:
ˆx(r) =
AT WA

−1AT W
n r2
M −r2
1, . . . , r2
M −r2
M−1
T
+Bk + a} ,
where
Bk =


w2
rM + 2rMwrM −w2
r1 −2r1wr1
...
w2
rM + 2rMwrM −w2
rM−1 −2rM−1wrM−1

.
Therefore
the
error
can
be
expressed
as
wr
=
AT WA

−1AT WBk , and its correlation is E

w2
r
	
=
AT WA

−1AT W E

BkBT
k
	
WT A
AT WA

−1. We let
C ≜E

BkBT
k
	
, whereas the lemma follows after some
algebraic manipulation.
APPENDIX C
PROOF OF LEMMA 3.4
Since ˜Vk and ˜φk are independent,
E
n
˜Vk cos

˜φk
o
= E
n
˜Vk
o
E
n
cos

˜φk
o
=Vk E {cos(φk + wφk)}
=Vk (cos(φk) E {cos(wφ,k)}
−sin(φk) E {sin(wφ,k)}) .
(52)
14
By using the cosine series,
E {cos(wφ,k)}= E
(
1 −
w2
φ,k
2
+
w4
φ,k
4!
+ · · · + (−1)n w2n
φ,k
(2n)!
)
= 1 −
σ2
φ
2 +
σ4
φ
22 · 2! + · · · + (−1)n σ2n
φ
2n · n! = e−
σ2
φ
2 ,
(53)
where we have used that ˜φk is Gaussian. Similarly, by using
the sine series, it results
E {sin(wφ,k)} = E
(
wφ −
w3
φ,k
3!
+
w5
φ,k
5!
−
w7
φ,k
7!
+ · · · + (−1)n w2n+1
φk
(2n + 1)!
)
= 0 .
(54)
Then Eq. (18) follows by substituting Eqs. (53) and (54) in
Eq. (52). Furthermore, since ˜Vk and ˜φk are independent, we
have
E
n
˜V 2
k cos2(˜φk)
o
= E
n
˜V 2
k
o
E
n
cos2(˜φk)
o
.
(55)
For the second factor of the right-end side of this expression
we have
E
n
cos2(˜φk)
o
= E
1
2 + 1
2 cos(2˜φk)

=1
2 + 1
2 E {cos (2 (φk + wφ,k))}
=1
2 + 1
2 E {cos (2φk) cos (2wφ,k)
−sin (2φk) sin (2wφ,k)}
=1
2 + 1
2 cos (2φk) E {cos (2wφ,k)} .
(56)
Using the cosine series, we write
E {cos (2wφ,k)} = E
(
1 −22 w2
φk
2! + 24 3w4
φk
4!
−26 5w6
φk
6!
+ · · · + (−1)n (2wφ)2n
(2n)!
)
=1 −2
σ2
φ
1! + 22 σ4
φ
2! −23 σ6
φ
3!
+ · · · + (−1)n 2nσ2n
φ
n!
= e−2σ2
φ .
(57)
By substituting (57) in (56), we obtain:
E
n
cos2(˜φk)
o
= 1
2 + 1
2 cos (2φk) e−2σ2
φ .
(58)
Eq. (19) follows by substituting (58) in (55), therefore con-
cluding the proof.
APPENDIX D
PROOF OF LEMMA 3.5
The expectation of (8) is
E {ˆx1,k+1|k+1} = (1 −β1,k) xk+1 + (1 −β1,k) E
n
w(r)
x1,k+1
o
+ β1,k

xk + E {wx1,k} + T Vk cos (φk) e−
σ2
φ
2

= (1 −β1,k) xk+1 + β1,kxk+1 + (1 −β1,k) E
n
w(r)
x1,k+1
o
+ β1,k E {wx1,k} + β1,kT Vk cos (φk)

e−
σ2
φ
2 −1

= xk+1 + E {wx1,k+1}.
whereby the lemma follows.
APPENDIX E
PROOF OF PROPOSITION 3.7
The cost function in (11) is always positive for any choice
of β1,k ∈R and of ρ1,k, and it is a parabola in β1,k, which
is due to that ak > 0. It follows that the cost function is
convex. Therefore, the optimal solution is given by computing
the derivative, and observing that the optimal solution must lie
in the interval B = [−1, 1].
APPENDIX F
PROOF OF LEMMA 4.1
We note that the speed at time k can be written as Vk = V0+
v3,1 + . . . + v3,k−1 . Since the process noise in the speed, v3,
is white, it follows that Vk ∼N
V0, (k −1) σ2
3


. Similarly,
for the orientation φk, we have that φk ∼N
φ0, (k −1) σ2
4


.
Then the Lemma follows by using Lemma 3.4, after simple
algebraic and trigonometric calculations.
APPENDIX G
PROOF OF PROPOSITION 4.2
First, we start by proving the following initial lemma:
Lemma G.1: Let a, b be two independent random variables
with non-zero average. Suppose that
lim
k→∞
E {(b −E {b})k}
( E {b})k
= 0 .
(59)
Then,
E
na
b
o
= E {a}
E {b}

1 +
∞
X
k=1
(−1)k E
n
(b −E {b})ko
( E {b})k


Proof: Since a and b are independent, we have that
E
na
b
o
= E {a} E
1
b

.
(60)
Furthermore,
E
1
b

= E

b−1	
= E
(
1 + b −E {b}
E {b}
−1)
1
E {b} .
(61)
15
By using Taylor series expansion as in [48, Section 5-4], we
can write
E
(
1 + b −E {b}
E {b}
−1)
=1+
∞
X
k=1
(−1)k E {(b −E {b})k}
( E {b})k
(62)
The sum in the previous equation converges as a consequence
of assumption (59). Finally, by substituting (62) in (61) and
(60), the lemma follows.
Now, using Lemma G.1 above, we write the expectation as
E

q2
q2 + z2

= E







1
1 + z2
q2







=
1
1 + E

z2
q2

×


1 +
∞
X
k=1
(−1)k
E
(z2
q2 −E
z2
q2
k)

1 + E

z2
q2
k


.
(63)
We are left with the derivation of E

z2/q2	
. We start noting
that, from Lemma G.1, we have
E
z2
q2

= E

z2	
E
 1
q2

= E {z2}
E {q2}
×

1 +
∞
X
k=1
(−1)k E
nq2 −E

q2	
ko
( E {q2})k

.
(64)
The variable q2 has a non-central chi-squared distribution of
1 degree of freedom, thus E

q2	
= σ2
q + µ2
q, and
E
nq2 −E

q2	
ko
=σ2k
q E
( q2
σ2q
−E
 q2
σ2q
k)
=σ2k
q µk .
(65)
This implies that
lim
k→∞
E
nq2 −E

q2	
ko
( E {q2})k
= lim
k→∞
σ2k
q µk
(µ2q + σ2q)k
= lim
k→∞
µk
 µ2q
σ2q + 1
k = 0 ,
since µk is bounded and the denominator grows unbounded,
and condition (59) holds, which is needed for the convergence
of the sum in Eq. (64). Therefore, Eq. (64) becomes
E
z2
q2

= µ2
z + σ2
z
µ2
q + σ2
q
"
1 +
∞
X
k=1
(−1)k
σ2k
q µk
µ2q + σ2q

k
#
.
(66)
We are left with the computation of the expectation in the
numerator of the summation in (63). This reduces to the k-th
central moment of a doubly noncentral F-distribution. With
this goal in mind, we first define χz ≜z2/σ2
z and χq ≜
q2/σ2
q. We note that χz and χq are distributed according to a
noncentral chi-squared distribution with 1 degree of freedom.
Therefore z2/q2 = χz/χq, since σ2
z = σ2
q, where the random
variable χz/χq is distributed according to a doubly noncentral
F-distribution. Then, we have that
E
(χz
χq
−E
χz
χq
k)
= µF,k .
(67)
The proposition follows by substituting (66) and (67) in
(63), and noting that as required by (59),
lim
k→∞
E
(z2
q2 −E
z2
q2
k)

1 + E

z2
q2
k
= 0 ,
since the numerator is given by (67) and is bounded, whereas
the denominator diverges because is the power of a number
greater than 1 due to that the expectation in the denominator
is positive.
APPENDIX H
PROOF OF PROPOSITION 4.3
The upper bound follows immediately from the fact that
q2/
q2 + z2

≤1, ∀q ∈R, ∀z ∈R. Regarding the lower
bound, by noting that ln (·) is a concave function and applying
Jensen’s inequality, we obtain that
e E{ln(q2)}−ln( E{q2}+ E{z2}) ≤E

q2
q2 + z2

.
(68)
Moreover,
E

ln
q2
	
=
−γe −ln (2) −2 ln (σq) −
M
0, 1
2, −µ2
q/2σ2
q


, whereas the proposition follows by sub-
stituting the previous equation in (68).
APPENDIX I
PROOF OF PROPOSITION 4.5
We begin by observing that ∀i, j the two matrices J(ub) and
J(lb) satisfy J(lb)
ij
≤Jij ≤J(ub)
ij
. Consider the element-wise
upper bound J(ub)
ij
. It follows that there exists a matrix G
such that Jij + Gij ≥J(ub)
ij
∀i, j and J + G = J(ub,G), where
G is symmetric and diagonally dominant with real and non-
negative elements. From the Gershgorin circle theorem, [49], it
follows that all its eigenvalues are non-negative, i.e., λ(G) ≥
0. Thus, λ
J(ub,G)

= λ (J) + λ (G) because both J and G
are Hermitian. Therefore, since J is positive semidefinite and
λ
J(ub,G)

> λ (J), J ⪯J(ub,G) and 0 ⪯J(ub,G). The proof
is concluded by noting that the same steps may be employed
for the lower bound J(lb,G).
REFERENCES
[1] N. Patwari, J. N. Ash, S. Kyperountas, A. O. Hero III, R. L. Moses, and
N. S. Correal, “Locating the nodes: cooperative localization in wireless
sensor networks,” IEEE Signal Processing Magazine, vol. 22, no. 4, pp.
54–69, July 2005.
[2] J. Rantakokko, J. Rydell, P. Strömbäck, P. Händel, J. Callmer, D. Törn-
qvist, F. Gustafsson, M. Jobs, and M. Gruden, “Accurate and reliable
soldier and first responder indoor positioning: multisensor systems and
cooperative localization,” IEEE Wireless Communications, vol. 18, no. 2,
pp. 10–18, Apr. 2011.
[3] H. Liu, H. Darabi, P. Banerjee, and J. Liu, “Survey of wireless indoor
positioning techniques and systems,” IEEE Transactions on Systems,
Man, and Cybernetics, Part C: Applications and Reviews, vol. 37, no. 6,
pp. 1067–1080, Nov. 2007.
16
[4] N. Patwari, I. Hero, A.O., M. Perkins, N. Correal, and R. O’Dea, “Rela-
tive location estimation in wireless sensor networks,” IEEE Transactions
on Signal Processing, vol. 51, no. 8, pp. 2137–2148, Aug. 2003.
[5] U. A. Khan, S. Kar, and J. M. F. Moura, “DILAND: An algorithm for
distributed sensor localization with noisy distance measurements,” IEEE
Transactions on Signal Processing, vol. 58, no. 3, pp. 1940–1947, Mar.
2010.
[6] X. Chen, A. Edelstein, Y. Li, M. Coates, M. Rabbat, and A. Men,
“Sequential Monte Carlo for simultaneous passive device-free tracking
and sensor localization using received signal strength measurements,”
in 10th International Conference on Information Processing in Sensor
Networks (IPSN), Apr. 2011, pp. 342–353.
[7] H. Wymeersch, J. Lien, and M. Win, “Cooperative localization in
wireless networks,” Proceedings of the IEEE, vol. 97, no. 2, pp. 427–
450, Feb. 2009.
[8] A. Conti, M. Guerra, D. Dardari, N. Decarli, and M. Win, “Network
experimentation for cooperative localization,” IEEE Journal on Selected
Areas in Communications, vol. 30, no. 2, pp. 467–475, February 2012.
[9] A. Bertrand and M. Moonen, “Seeing the bigger picture: How nodes
can learn their place within a complex ad hoc network topology,” IEEE
Signal Processing Magazine, vol. 30, no. 3, pp. 71–82, 2013.
[10] M. Rabbat, R. Nowak, and J. Bucklew, “Robust decentralized source
localization via averaging,” in IEEE International Conference on Acous-
tics, Speech, and Signal Processing, vol. 5, 2005, pp. 1057–1060.
[11] U. Khan, S. Kar, and J. Moura, “Distributed sensor localization in
random environments using minimal number of anchor nodes,” IEEE
Transactions on Signal Processing, vol. 57, no. 5, pp. 2000–2016, 2009.
[12] U. Hammes, E. Wolsztynski, and A. Zoubir, “Robust tracking and ge-
olocation for wireless networks in NLOS environments,” IEEE Journal
of Selected Topics in Signal Processing, vol. 3, no. 5, pp. 889–901, 2009.
[13] J. Wang, Q. Gao, Y. Yu, H. Wang, and M. Jin, “Toward robust
indoor localization based on bayesian filter using chirp-spread-spectrum
ranging,” IEEE Transactions on Industrial Electronics, vol. 59, no. 3,
pp. 1622–1629, 2012.
[14] D. Zachariah, A. De Angelis, S. Dwivedi, and P. Handel, “Schedule-
based sequential localization in asynchronous wireless networks,”
EURASIP Journal on Advances in Signal Processing, no. 1, p. 16, 2014.
[Online]. Available: http://asp.eurasipjournals.com/content/2014/1/16
[15] S. Gezici, Z. Tian, G. Giannakis, H. Kobayashi, A. Molisch, H. Poor,
and Z. Sahinoglu, “Localization via ultra-wideband radios: a look at
positioning aspects for future sensor networks,” IEEE Signal Processing
Magazine, vol. 22, no. 4, pp. 70–84, July 2005.
[16] Y. Shen and M. Win, “Fundamental limits of wideband localization part
I: A general framework,” IEEE Transactions on Information Theory,
vol. 56, no. 10, pp. 4956–4980, Oct. 2010.
[17] M. Spirito, “On the accuracy of cellular mobile station location estima-
tion,” IEEE Transactions on Vehicular Technology, vol. 50, no. 3, pp.
674–685, May 2001.
[18] B. Ristic and M. S. Arulampalam, “Tracking a manoeuvring target
using angle-only measurements: algorithms and performance,” Signal
Processing, vol. 83, no. 6, pp. 1223–1238, June 2003.
[19] N. Bergman, “Recursive Bayesian estimation navigation and tracking
applications,” Ph.D. dissertation, Linköping Studies in Science and
Technology, Sweden, 1999, dissertation n.579.
[20] L. Arienzo and M. Longo, “Posterior Cramér-Rao bound for range-
based target tracking in sensor networks,” in IEEE/SP 15th Workshop
on Statistical Signal Processing (SSP), Sept. 2009, pp. 541–544.
[21] K. Cheung, H. So, W.-K. Ma, and Y. Chan, “Received signal strength
based mobile positioning via constrained weighted least squares,” in
IEEE International Conference on Acoustics, Speech, and Signal Pro-
cessing (ICASSP), vol. 5, Apr. 2003, pp. V – 137–140.
[22] A. Speranzon, C. Fischione, K. Johansson, and A. Sangiovanni-
Vincentelli, “A distributed minimum variance estimator for sensor net-
works,” IEEE Journal on Selected Areas in Communications, vol. 26,
no. 4, pp. 609–621, May 2008.
[23] F. Gustafsson, Statistical sensor fusion.
Studentliteratur, 2010.
[24] I. Skog and P. Händel, “In-car positioning and navigation technologies
- a survey,” IEEE Transactions on Intelligent Transportation Systems,
vol. 10, no. 1, pp. 4–21, Mar. 2009.
[25] T. Kailath, A. H. Sayed, and B. Hassibi, Linear Estimation, 1st ed.
Prentice Hall, Apr. 2000.
[26] D. B. Jourdan, J. Deyst, J. J., M. Z. Win, and N. Roy, “Monte Carlo
localization in dense multipath environments using UWB ranging,” in
IEEE International Conference on Ultra-Wideband (ICUWB), 5–8 Sept.
2005, pp. 314–319.
[27] J. Hol, F. Dijkstra, H. Luinge, and T. Schön, “Tightly coupled
UWB/IMU pose estimation,” in IEEE International Conference on
Ultra-Wideband (ICUWB), Sept. 2009, pp. 688 –692.
[28] A. De Angelis, J. O. Nilsson, I. Skog, P. Händel, and P. Carbone, “Indoor
positioning by ultrawide band radio aided inertial navigation,” Metrology
and Measurement Systems, vol. XVII, pp. 447–460, 2010.
[29] R. Niu, R. Blum, P. Varshney, and A. Drozd, “Target localization and
tracking in noncoherent multiple-input multiple-output radar systems,”
IEEE Transactions on Aerospace and Electronic Systems, vol. 48, no. 2,
pp. 1466 –1489, april 2012.
[30] A. De Angelis, C. Fischione, and P. Händel, “A sensor fusion algorithm
for mobile node localization,” in 18th World Congress of the Interna-
tional Federation of Automatic Control (IFAC), Milan, Italy, Aug. 28 –
Sept. 2 2011.
[31] A. De Angelis and C. Fischione, “A distributed information fusion
method for localization based on Pareto optimization,” in 7th IEEE
International Conference on Distributed Computing in Sensor Systems
(DCOSS), Barcelona, Spain, June 27–29 2011.
[32] P. Tichavský, C. Muravchik, and A. Nehorai, “Posterior Cramér-Rao
bounds for discrete-time nonlinear filtering,” IEEE Transactions on
Signal Processing, vol. 46, no. 5, pp. 1386–1396, May 1998.
[33] A. De Angelis, M. Dionigi, R. Moschitta, A. Giglietti, and P. Carbone,
“Characterization and modeling of an experimental UWB pulse-based
distance measurement system,” IEEE Transactions on Instrumentation
and Measurement, vol. 58, no. 5, pp. 1479–1486, May 2009.
[34] M. Gholami, S. Gezici, and E. Strom, “Tdoa based positioning in the
presence of unknown clock skew,” Communications, IEEE Transactions
on, vol. 61, no. 6, pp. 2522–2534, June 2013.
[35] A. Mourikis and S. I. Roumeliotis, “Performance analysis of multirobot
cooperative localization,” IEEE Transactions on Robotics, vol. 22, no. 4,
pp. 666–681, Aug. 2006.
[36] A. Georgiev and P. K. Allen, “Localization methods for a mobile robot
in urban environments,” IEEE Transactions on Robotics, vol. 20, no. 5,
pp. 851–864, Oct. 2004.
[37] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge
University Press, 2004.
[38] I. Guvenc, S. Gezici, and Z. Sahinoglu, “Fundamental limits and im-
proved algorithms for linear least-squares wireless position estimation,”
Wireless Communications and Mobile Computing, 2010.
[39] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation
Theory.
Englewood Cliffs, NJ: Prentice Hall, 1993.
[40] J. Taylor, “The Cramér-Rao estimation error lower bound computation
for deterministic nonlinear systems,” IEEE Transactions on Automatic
Control, vol. 24, no. 2, pp. 343–344, Apr. 1979.
[41] B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman Filter:
Particle Filters for Tracking Applications.
Artech House, 2004.
[42] T. Bréhard and J.-P. Le Cadre, “Closed-form posterior Cramér-Rao
bound for bearings-only tracking,” IEEE Transactions on Aerospace and
Electronic Systems, vol. 42, no. 4, pp. 1198–1223, 2006.
[43] A. Bishop, B. Fidan, B. Anderson, K. Dogancay, and P. Pathirana, “Op-
timality analysis of sensor-target localization geometries,” Automatica,
vol. 46, no. 3, pp. 479–492, Mar. 2010.
[44] I. Shames, B. Fidan, B. Anderson, and H. Hmam, “Cooperative self-
localization of mobile agents,” IEEE Transactions on Aerospace and
Electronic Systems, vol. 47, no. 3, pp. 1926–1947, July 2011.
[45] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions
with Formulas, Graphs, and Mathematical Tables.
Dover Publications,
1972.
[46] E. Abbott and D. Powell, “Land-vehicle navigation using GPS,” Pro-
ceedings of the IEEE, vol. 87, no. 1, pp. 145–162, Jan. 1999.
[47] S. Julier and J. Uhlmann, “Unscented filtering and nonlinear estimation,”
Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, Mar. 2004.
[48] A. Papoulis, Probability, Random Variables and Stochastic Processes.
McGraw-Hill Companies, February 1991.
[49] R. Horn and C. Johnson, Matrix analysis.
Cambridge University Press,
1990.