Keywords: Data fusion, State Estimation, Dempster-Shafer theory, Evidence theory, Interval analysis, Bounded error estimation, Localization

A new approach to non-linear state estimation based on belief function theory and interval analysis is presented. This method uses belief structures composed of a finite number of axis-aligned boxes with associated masses. Such belief structures can represent partial information on model and measurement uncertainties, more accurately than can the bounded error approach alone. Focal sets are propagated in system equations using interval arithmetics and constraint satisfaction techniques, thus generalizing pure interval analysis. This model was used to locate a land vehicle using a dynamic fusion of GPS measurements with dead reckoning sensors. The method has been shown to provide more accurate estimates of vehicle position than does the bounded error method while retaining what is essential: providing guaranteed computations. The performances of our method were also slightly better than those of a particle filter, with comparable running time. These results suggest that our method is a viable alternative to both bounded error and probabilistic Monte-Carlo approaches for vehicle localization applications.

Introduction

Estimating the state of a dynamical system based on noisy sensor measurements is a common problem in various areas such as mobile robot navigation or target tracking [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] [START_REF] Gustafsson | Particle filters for positioning, navigation, and tracking[END_REF]. If several sensors or other information sources are available, then data fusion algorithms have to be used. A critical issue when designing such systems is the representation of uncertainties that pervade both sensor measurements and the state space model itself. Two broad categories of state estimation methods can be distinguished in that respect.

Methods in the first category are based on a probabilistic description of uncertainty and assume measurement noise and state perturbations to be realizations of random variables with known statistical properties. This is the mainstream approach. In a linear context, the most popular method is by far the application of the Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], which assumes model and observation noise to be normal. In non-linear cases, the two main methods are based on the Extended Kalman filter (EKF) and the Particle filter (PF). The EKF linearizes the state and measurement equations and then applies the Kalman filter to obtain state estimates [START_REF] Bar-Shalom | Multitarget-Multisensor Tracking: Applications and Advances[END_REF][START_REF] Gelb | Applied Optimal Estimation[END_REF], assuming process and observation noises to be normal. The state posterior probability distribution is approximated using a Gaussian distribution that is propagated analytically through linearized system equations. However, linearization is inherently local and may fail to produce reliable estimates, especially when the state model is highly non-linear. The statistical interpretation of covariance matrices is also unclear in this approach, as the statistical properties of perturbations though non-linear equations are not well-known. Recently, sequential Monte Carlo methods for recursive Bayesian filtering, or particle filter (PF) methods have emerged as useful tools for problems requiring non-linear dynamic state estimation [START_REF] Chen | Bayesian Estimation by Sequential Monte-Carlo Sampling[END_REF][15] [START_REF] Gustafsson | Particle filters for positioning, navigation, and tracking[END_REF]. The PF estimator usually provides more accurate information about the state posterior probability distribution than does the EKF, especially if it takes a multimodal shape or if noise distributions are non-Gaussian. The efficiency and accuracy of the PF method depends mostly on the number of particles used in the estimation and on the propagation used to re-allocate weights to these particles at each iteration. To cope with high uncertainty of measurement, a large number of particles have to be used, especially when the number of dimensions of the state vector is high; this may be an issue for real-time implementation of the PF method. Several authors have tried to overcome these shortcomings by combining approaches (see [START_REF] Gustafsson | Particle filters for positioning, navigation, and tracking[END_REF] and references therein) or by adapting the size of sample sets during the estimation process [START_REF] Fox | Adapting the sample size in particle filters through KLD-sampling[END_REF].

The second group of methods corresponds to state bounding. Assuming that all variables belong to known compact sets, these methods attempt to build simple sets, such as ellipsoids or boxes, guaranteed to contain all state vectors consistent with given constraints [START_REF] Jaulin | Applied Interval Analysis[END_REF]. In a linear context, this approach has been investigated by several authors since the 1960's (see, e.g., [START_REF] Bertsekas | Recursive state estimation for a set-membership description of uncertainty[END_REF][START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Maksarov | State bounding with ellipsoidal set description of the uncertainty[END_REF][START_REF] Schweppe | Recursive state estimation: unknown but bounded errors and system inputs[END_REF][START_REF] Servi | Recursive estimation in the presence of uniformly distributed measurement noise[END_REF][START_REF] Witsenhausen | Sets of possible states of linear systems given perturbed observations[END_REF] and references therein). In the non-linear context, the methodology is less developed. Non-linear state-bounding methods have been presented in [START_REF] Ashokaraj | Sensor based robot localisation and navigation: Using interval analysis and extended Kalman filter[END_REF] [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][26] [START_REF] Kieffer | Guaranteed recursive nonlinear state estimation using interval analysis[END_REF]. Recently, a relatively simple and fast method based on constraint propagation and interval analysis was introduced and applied to vehicle localization by [START_REF] Gning | Dynamic vehicle localization using constraints propagation techniques on intervals. a comparison with Kalman filtering[END_REF] [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF]. The main interest in the bounded-error approach arises from the fact that it allows so-called validated computations: computations with guaranteed accuracy taking into account all possible sources of error, from data imprecision to rounding errors due to computer calculations. The major drawback of this approach is the difficulty of determining noise bounds. Indeed, if the bounds are too tight, the data may become inconsistent with the system equations, in which case, the method fails to provide a solution. On the other hand, if the bounds are overestimated, the estimated state becomes very imprecise, and the method becomes overly pessimistic. A related problem with this approach is its lack of robustness against outlying observations, which may be a severe limitation in cases where sensor information is unreliable.

In this paper, we propose to replace the set-based representation of uncertainty used in the state bounding approach by a more general formalism based on belief functions. The theory of belief functions, also known as Dempster-Shafer (DS) or Evidence theory, was introduced by Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] in statistical inference, and was further developed by Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF] and Smets [START_REF] Ph | The Transferable Belief Model[END_REF] as a theory of epistemic uncertainty. The theory of belief functions makes it possible to model various states of knowledge, ranging from complete ignorance to probabilistic uncertainty [START_REF] Yager | Classic Works of the Dempster-Shafer Theory of Belief Functions[END_REF]. A belief function may be viewed both as a generalized set [START_REF] Dubois | Possibility Theory: An approach to computerized processing of uncertainty[END_REF] and as a generalized probability measure: DS theory thus encompasses set-based as well as probabilistic formalisms.

Until now, the application of belief functions to state estimation has been limited. Dempster [START_REF] Dempster | Normal belief functions and the Kalman filter[END_REF] showed that computations in the Kalman filter can be interpreted in a graphical model in terms of the propagation of Gaussian belief functions (a special kind of belief function characterized by a normal density over linear subspaces). Ristic and Smets [START_REF] Ph | Kalman filter and joint tracking and classification based on belief functions in the TBM framework[END_REF] have provided another interpretation of the Kalman filter under a belief function paradigm and have presented an approach to joint tracking and classification based on continuous belief functions. Mahler [START_REF] Mahler | Can the Bayesian and Dempster-Shafer approaches be reconciled? Yes[END_REF] has also discussed various extensions of the Kalman filter that make it possible to handle uncertainties expressed by belief functions. The main motivations of this study have been to cope with situations where the available information about the system under study or about sensor measurements is too weak to be easily handled using a purely probabilistic approach. However, these situations have been confined mostly to linear models; this may be explained by the fact that the manipulation of continuous belief functions is delicate and quickly becomes intractable when applied to more complex models.

In this paper, we propose a different approach -one that generalizes the bounding approach to state estimation. Taking as a starting point the bounded error state estimation method proposed in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], we present a data filtering and fusion method that uses Dempster-Shafer mass functions to assign belief masses to a finite number of focal sets, taking the form of axis-aligned boxes. Such mass functions can be seen as "generalized boxes" composed of a collection of boxes with associated weights. Focal sets are propagated in the system equations using interval arithmetics and constraint satisfaction techniques [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Applied to state estimation in dynamical systems, this approach allows us to compute a belief function for the system state at each time step. This approach extends the pure interval approach, making it more robust and accurate, while retaining the validity of computations.

The article is organized as follows. Section 2 first presents the background on interval analysis and bounded error state estimation. The necessary concepts in belief function theory are then recalled in Section 3. Our approach is introduced in Section 4, and in Section 5 we apply our method to dynamic land vehicle localization using GPS, gyrometer and odometer measurements. Finally, we conclude and discuss the main contributions of the paper in Section 6.

Interval Analysis

In this subsection, we briefly recall the main concepts pertaining to interval analysis and constraint satisfaction techniques that will be used in the rest of this work.

Basic Concepts

Usually, interval analysis is used to model quantities that vary around a central value within certain bounds [START_REF] Jaulin | Applied Interval Analysis[END_REF]. A (closed) real interval, denoted [x], is defined as a closed and connected subset of R:

[x] = [x, x] = {x ∈ R | x ≤ x ≤, x}, (1)
where x and x are the minimal and maximal bounds of [x]. Set theory operations, such as intersection or union, can be applied to intervals. The intersection of two intervals is always an interval, whereas the union is not. The interval union of two intervals [x] and [y] can be defined as the interval hull of [x] ∪ [y], i.e., the smallest interval containing [x] ∪ [y], which is denoted as follows [START_REF] Jaulin | Applied Interval Analysis[END_REF]:

[x] ⊔ [y] = [[x] ∪ [y]] . (2)
The four classical arithmetic operations can be extended to intervals. For any such binary operator ⋄, the interval [x] ⋄ [y] is defined for any intervals [x] and [y] as [START_REF] Jaulin | Applied Interval Analysis[END_REF] [x]

⋄ [y] = [{x ⋄ y ∈ R|x ∈ [x], y ∈ [y]}] . (3)
When ⋄ is continuous, as is the case for the arithmetic operations, the set {x ⋄ y ∈ R|x ∈ [x], y ∈ [y]} is an interval; we thus have the following equivalent definition for

[x] ⋄ [y]: [x] ⋄ [y] = {x ⋄ y ∈ R|x ∈ [x], y ∈ [y]}. ( 4
)
When ⋄ is one of the four arithmetic operations, then the bounds of [x] ⋄ [y] can be computed from the bounds of [x] and [y] using simple formulas. The following equations hold [START_REF] Jaulin | Applied Interval Analysis[END_REF][33]:

[x] + [y] = [x + y, x + y] (5) [x] -[y] = [x -y, x -y] (6)
If the bounds of [x] and [y] are positive, we also have

[x]/[y] = [x/y, x/y]. (8)
Elementary functions such as exp, ln, sin, etc., can also be extended to intervals. If f is a function from R to R, its interval counterpart [f ] is defined by:

[f ]([x]) = [{f (x)|x ∈ [x]}]. (9)
For non-monotonic functions, however, the computation of [f ]([x]) is usually not as simple.

An interval real vector [x] of R n , also called a box, is defined as a Cartesian product of n intervals:

[x] = [x 1 ] × . . . × [x n ]. (11)
Geometrically, it may be figured as an n dimensional, axis-aligned parallelepiped. The set of n-dimensional interval real vectors will be denoted IR n . The notions recalled above may be easily extended to boxes. For instance, the interval hull of a subset A of R n is the smallest box of R n containing A. It is noted [A]. The interval union of two boxes [x] and [y] may be computed as:

[x] ⊔ [y] = ([x 1 ] ⊔ [y 1 ]) × . . . × ([x n ] ⊔ [y n ]). (12)
Classical operations for real vectors can be extended to interval vectors. For instance,

[x] + [y] = ([x 1 ] + [y 1 ]) × . . . × ([x n ] + [y n ]). ( 13
)
Let f be a function from R n to R m . It is often useful to compute the image f ([x]) of a box [x] by f . This image is usually not a box (see Figure 1), and its expression may be difficult to obtain. However, it may sometimes be approximated by an inclusion functi on. An interval function [f ] from IR n to IR m is said to be an inclusion function

for f if f ([x]) ⊆ [f ]([x]), ∀[x] ∈ IR n . ( 14
)
Inclusion functions may be very pessimistic, as shown by Figure 1. The inclusion function

[f ] is minimal if, for any x, [f ]([x]) is the interval hull of f ([x]
). The minimal inclusion function for f is unique and may be noted as [f ] * . Finding inclusion functions that can be evaluated reasonably quickly and such that [f ]([x]) is not too large for most x is one of the main purposes of interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Different algorithms have been proposed to reduce the size of boxes enclosing f ([x]). For the fusion problem considered here, we have chosen to use constraint propagation techniques [START_REF] Jaulin | Applied Interval Analysis[END_REF], because of the great redundancy of data and equations.

Contraint Satisfaction Problems

Consider n variables x i , i ∈ {1, . . . , n} linked by m relations (or constraints) of the form

f j (x 1 , . . . , x n ) = 0, j = 1, . . . , m. (15)
If we denote by x the vector (x 1 , x 2 , • • • , x n ) T and by f the function whose coordinate functions are the f i s:

f = (f 1 , f 2 , • • • , f m ) T
, we can write these m constraints in vector notation as

f (x) = 0. (16)
Let us assume that vector x is known to belong to some prior domain [x], and we want to compute the set of all x in the prior domain verifying the constraints. This defines a constraint satisfaction problem (CSP), which can be denoted as

H : (f (x) = 0, x ∈ [x]). ( 17
)
The solution set of H is defined as: [START_REF] Jaulin | Applied Interval Analysis[END_REF]Chapter 4], including Gauss elimination, the Gauss-Seidel algorithm and linear programming, among others. Each of these methods may be more suitable for some types of CSP. Although the approach introduced in this paper is not bound to any particular contractor, we will focus on the forward-backward propagation (FBP) technique [START_REF] Waltz | Generating semantic descriptions from drawings of scenes with shadows[END_REF][9][10] [26][25]. This method was chosen because it is known to be especially efficient in cases of high redundancy of data and equations -a situation that will be encountered in the application described in Section 5. Furthermore, this method can handle any kind of non-linearity and its implementation is quite simple. However, other contractors could be used as well.

S = {x ∈ [x] | f (x) = 0}. (18)
The FBP contractor is based on the propagation of primitive constraints for real variables. A constraint is said to be primitive if it involves a single operator (such as +, -, * or \) or a single function (such as cos, sin or sinh). The method proceeds by contracting H with respect to each of the primitive constraints until convergence on a minimal domain. The complete description of the FBP algorithm goes beyond the scope of this paper. It can be found in [26, page 77]. Here, it will simply be described though an example.

Example 1 Let us consider the constraint z = x exp(y) and the initial domain

[z] = [0, 3], [x] = [1, 7] and [y] = [0, 1].
The constraint is not primitive, as it involves one arithmetic operation and one elementary function. However, it can be decomposed into two primitive constraints by introducing an auxiliary variable a:

a = exp(y) z = a • x. ( 19
)
The domain of a may be initialized by [a] = [0, +∞). The FBP algorithm alternates between two phases, called forward and backward propagation.

In the forward propagation phase, the domains of a and z are contracted using [START_REF] Fox | Adapting the sample size in particle filters through KLD-sampling[END_REF]. This is achieved in two steps F 1 and F 2 :

• F 1 : [a] ← [a] ∩ [exp]([y]) = [0, +∞) ∩ [1, e] = [1, e] • F 2 : [z] ← [z] ∩ [x] • [a] = [0, 3] ∩ [1, 7] • [1, e] = [1, 3].
Then, backward propagation is performed, updating the domains associated with all variables. This is achieved in three steps B 3 , B 4 and B 5 :

• B 3 : [x] ← [x] ∩ ([z]/[a]) = [1, 7] ∩ [1, 3]/[1, e] = [1, 3] • B 4 : [a] ← [a] ∩ ([z]/[x]) = [1, e] ∩ [1, 3]/[1, 3] = [1, e] • B 5 : [y] ← [y] ∩ [ln]([a]) = [0, 1] ∩ [0, 1] = [0, 1].

After one forward-backward propagation cycle, the domains of the variables have thus been reduced to

[z] = [1, 3], [x] = [1, 3] and [y] = [0, 1].
It can be verified that these domains are no longer changed after another iteration of the procedure. Figure 3 shows a graphical representation of the problem, where each node corresponds to a variable; there is a link between two nodes if the corresponding variables are linked by a primitive constraint. It may be shown that the FPB contractor yields the smallest possible domain if this graph contains no cycle [START_REF] Jaulin | Guaranteed nonlinear estimation using constraint propagation on sets[END_REF].

Bounded Error Estimation

In this section, we briefly present the bounded error state estimation method (BEE) introduced in [21] [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF]. Consider a non-linear dynamical system represented in state space by the following discrete-time equations:

x k+1 = f (x k , u k ) ( 20
)
z k+1 = g(x k+1 ), (21)
where f : R nx × R nu × R nv -→ R nx is a possibly non-linear function that relates the state x k+1 at time k +1 to the previous state x k at time k with the input u k . Function g : R nx × R nw -→ R nz expresses the relationship between the observed output z k+1 and the state x k+1 . Note that we could easily consider a more complex model including, for instance, additional state and measurement perturbation variables. However, the simple model above will be sufficient to illustrate this approach, as well as for the application we develop in Section 5.

When working with bounded errors, the initial state x 0 is assumed to belong to some prior box [x 0 ] ⊂ IR nx and, at each time step k, the input u k and the measurement vector z k are assumed to belong to boxes [u k ] and [z k ], respectively. The problem considered in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF] is to provide at each time step k + 1 a guaranteed estimate of the state x k+1 , in the form of a box [x k+1 ] that certainly contains the true state vector x k+1 , assuming the system equations (20-21) to be correct and taking into account all the information available at time step k. This problem can be seen as a CSP

H 0:k+1 : (F(ξ 0:k+1 ) = 0, ξ 0:k+1 ∈ [ξ 0:k+1 ]), (22)
with ξ 0:k+1 = (x k+1 , . . . , x 0 , u k , . . . , u 0 , z k+1 , . . . , z 1 ) and

F : ξ 0:k+1 → (x k+1 -f (x k , u k ), . . . , x 1 -f (x 0 , u 0 ), z k+1 -g(x k+1 ), . . . , z 1 -g(x 1 )) . ( 23
)
The FBP contractor described in the previous section can be applied to this problem after all the constraints in H 0:k+1 have been decomposed into primitive ones. The forward and backward propagations are then iterated until convergence.

In practice, however, considering all equations from step 0 to k + 1 is often not feasible in real-time, as noted in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF]. This problem can be solved by considering only variables and constraints within a limited time window. In the simplest (and fastest) variation considered here, we start from the previous state estimate [x k ] at step k, and we determine [

x k+1 ] from [x k ], [u k ], [z k+1
] and the constraints ( 20)-( 21) between time steps k and k + 1. Formally, this corresponds to the following CSP:

H k+1 : (F(ξ k+1 ) = 0, ξ k+1 ∈ [ξ k+1 ]), (24)
with ξ k+1 = (x k+1 , x k , u k , z k+1 ) and

F : ξ k+1 → (x k+1 -f (x k , u k ), z k+1 -g(x k+1 )) . ( 25
)
The method is summarized in Algorithm 1.

Algorithm 1 Bounded Error Estimation (BEE) algorithm

Input: % Forward propagation % 6:

[x k ],[u k ], [z k+1 ] Output: [x k+1 ] 1: Read the input
[x k+1 ] ← [x k+1 ] ∩ [f ]([x k ], [u k ]) 7: [z k+1 ] ← [z k+1 ] ∩ [g]([x k+1 ]) 8:
% Backward propagation %

9:

Back-propagate the constraint

x k+1 -f (x k , u k ) = 0 10:
Back-propagate the constraint z k+1g(x k+1 ) = 0 11: end while

In steps 1 and 2, intervals are determined around the input u k and the output z k+1 using error bounds. These bounds depend on the sensors used and are assumed to be known. As noted in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], they may sometimes be determined from constructor data sheets, or they may be derived from estimated imprecision provided by the sensors in real time.

Initially, nothing is known about x k+1 : [x k+1 ] is thus initialized as the largest possible box in step 3. Boxes [x k+1 ], [x k ], [u k ] and [z k+1 ] are then contracted using the FBP algorithm described in the previous section. Step 6 corresponds to a prediction step, as it computes a new set estimate of x k+1 from the current set estimates of x k and u k . We can remark that an inclusion function [f ] for f has to be used, as the image of [x k ] and [u k ] by f is usually not a box. Similarly, step 7 computes the intersection between the "predicted output" [g]([x k+1 ]) computed from [x k+1 ] using an inclusion function [g] for g and the measured output. These two steps are illustrated in Figure 4.

Steps 9 and 10 correspond to the backward propagation phase of the FBP algorithm. They consist in contracting all boxes, including [x k ] and [u k ], by expressing the corresponding variables as functions of the other variables. As already noted in Section 2.2, the FBP algorithm requires that all constraints be decomposed into primitive ones.

The algorithm stops when no interval is contracted during the most recent iteration. The function computed by this algorithm will be denoted as follows:

[x k+1 ] = C (f , g, [x k ], [u k ], [z k+1 ]) . ( 26
)
As noted in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], the main difficulty in using this algorithm, and the bounded error approach in general [START_REF] Jaulin | Applied Interval Analysis[END_REF]Chapter 6], is the determination of error bounds. If they are overestimated, then the estimation may be overly pessimistic, resulting in highly imprecise state estimates. On the other hand, if they are underestimated, then inconsistencies can occur in step 7 of Algorithm 1, as the intersection between the predicted and measured output can become empty, in which case the procedure has to be stopped.

In this paper, we propose to cope with this problem by using a more general formalism to describe the uncertainty of the variables of interest. This formalism, referred to as the theory of belief functions or the Dempster-Shafer theory, can be seen as generalizing both probability theory and set theory. It will allow us to describe errors no longer by single boxes, but by families of boxes with associated belief masses. This formalism is introduced in the next section.

Theory of Belief Functions

The Dempster-Shafer (DS) theory of belief functions [START_REF] Shafer | A mathematical theory of evidence[END_REF] is a formal tool suitable for representing partial knowledge as well as uncertain and imprecise information. It is based on the formalism of belief functions, which generalize both set-based representations (including interval analysis) and probabilistic representations. In this section, we introduce the main concepts of this theory and some necessary notions that will be used in the proposed approach.

Basic Definitions

Let ω denote a variable taking values in a domain Ω called the frame of discernment. Let A 1 , . . . , A p be p subsets of Ω. A discrete mass function m with focal sets A 1 , . . . , A p is a function from 2 Ω to [0, 1] verifying m(A i ) = m i > 0 for all i ∈ {1, . . . , p}, m(A) = 0 for all A ∈ {A 1 , . . . , A p }, and:

p i=1 m i = 1. ( 27
)
A mass function is said to be normal if m(∅) = 0, categorical if p = 1 and consonant if the focal sets are nested. In the following, all mass functions will be assumed to be normal, unless otherwise specified.

In most presentations of DS theory, Ω is assumed to be finite. However, the theory remains basically unchanged if Ω is infinite (even uncountable), as long as the number of focal sets remains finite. If Ω = R, the focal sets are usually assumed to be intervals [START_REF] Ferson | Constructing probability boxes and Dempster-Shafer structures[END_REF][45] [START_REF] Petit-Renaud | Nonparametric regression analysis of uncertain and imprecise data using belief functions[END_REF]. In the multidimensional case where Ω = R n , this approach can be extended by assuming focal sets to be n-dimensional boxes. The case of a discrete probability distribution is recovered when each A i is reduced to a point. We also note that categorical mass functions are in one-to-one correspondence with subsets of Ω: to each subset A ⊆ Ω corresponds a mass function m A such that m A (A) = 1. A mass function can thus be seen as a generalized set [START_REF] Dubois | A set-theoretic view of belief functions: logical operations and approximations by fuzzy sets[END_REF]. This observation is the basis of the extension of interval analysis into the realm of belief functions, which will be presented in the next section.

Given a normal mass function m with focal sets A 1 , . . . , A p , the corresponding belief and plausibility functions are defined, respectively, as:

bel(A) = {i|A i ⊆A} m i (28) pl(A) = {i|A i ∩A =∅} m i , (29)
for all A ⊆ Ω. In particular, if m is consonant, it may be shown [START_REF] Shafer | A mathematical theory of evidence[END_REF] that

pl(A ∪ B) = max(pl(A), pl(B)), ∀A, B ⊆ Ω. ( 30
)
In this case, the plausibility function pl associated with m is thus a possibility measure, in the sense of possibility theory as introduced by Zadeh [47] [START_REF] Dubois | Possibility Theory: An approach to computerized processing of uncertainty[END_REF]. The corresponding possibility distribution π is then the function from Ω to [0, 1] defined by π(x) = pl({x}), for all x ∈ Ω. Mass functions may be used to model different situations of uncertainty about the values of variables describing a given problem. In the so-called multivalued mapping framework initially introduced by Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], we have an underlying space Ψ = {ψ 1 , . . . , ψ p } and a multivalued mapping Γ : Ψ → 2 Ω such that Γ(ψ i ) = A i , for i = 1, . . . , p, where, for simplicity, the A i s are assumed to be all distinct. We assume our beliefs on Ψ to be modeled by a probability measure P , with P (ψ i ) = m i . For instance, P might represent some stochastic mechanism that selects an element from Ψ. If ψ i is selected, all we know about ω is that it is in A i . Our beliefs concerning Ω may then be represented by the mass function induced by P through the multivalued mapping Γ, so that m(Γ(ψ i )) = P (ψ i ) = m i for all i ∈ {1, . . . , p}. We then have bel

(A) = P ({ψ i ∈ Ψ | Γ(ψ i ) ⊆ A}) (31)
and

pl(A) = P ({ψ i ∈ Ψ | Γ(ψ i ) ∩ A = ∅}), (32)
for all A ⊆ Ω. The elements of Ψ may be viewed as different ways of explaining or interpreting the evidence that we have for a given problem, or as different mechanisms that may have produced this evidence. In some cases, they may be given a precise physical interpretation. For instance, assume that ω is a measurement error and Ω ⊆ R. The sensor that delivers the measurement can be in the nominal state (ψ 1 ) or in a defective state (ψ 2 ). In the first case, all we know is that |ω| ≤ b, whereas in the second case ω can take any value in Ω. The probability that the sensor is in the nominal state is m 1 . This information about ω may be represented as a mass function m such that m([-b, b]) = m 1 and m(Ω) = 1m 1 .

Extending Interval Analysis to Mass Functions

In [START_REF] Yager | Arithmetic and other operations on Dempster-Shafer structures[END_REF], Yager proposed a simple scheme for extending arithmetic and other operations to mass functions (see also [START_REF] Regan | Equivalence of methods for uncertainty propagation of real-valued random variables[END_REF]). This approach will be summarized in this section. Because our goal here is to generalize interval analysis, all mass functions will be assumed to have interval focal sets.

Let us assume that we have two real variables x and y, and our knowledge about these variables is modeled by mass functions denoted by m x and m y , respectively. We are interested in the value of z = x⋄y, where ⋄ is a binary operator from R 2 to R. If m x and m y are categorical and such that m x ([x]) = 1 and m y ([y]) = 1 for some intervals [x] and [y] in IR, then interval analysis can be applied and our knowledge about z can be represented by a categorical mass function m z with focal interval [x] ⋄ [y], defined by (3). In the more general case, where m x and m y are mass functions with arbitrary numbers of focal intervals, then some extension of interval analysis has to be proposed, based on additional assumptions.

Let us assume that m x and m y are induced by multivalued mappings Γ x and Γ y from underlying spaces Ψ x and Ψ y , with [x i ] = Γ x (ψ x i ) for i = 1, . . . , p x and [y j ] = Γ y (ψ y j ) for j = 1, . . . , p y . Let us further assume that our beliefs about Ψ x and Ψ y are modeled by probability measures P x and P y , and that the joint probability measure P x,y for Ψ x × Ψ y is simply the product of P x and P y :

P x,y (ψ x i , ψ y j ) = P x (ψ x i )P y (ψ y j ) = m x i m y j , (33)
for all i and j. If P x and P y correspond to some stochastic mechanisms selecting elements from Ψ x and Ψ y , then [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF] simply means that these two mechanisms are independent. Under these assumptions, what can we say about z? If we know that ψ x i and ψ y j have been selected, then all we know about x, y and z

is that x ∈ [x i ], y ∈ [y j ], and z = x ⋄ y ∈ [x i ] ⋄ [y j ].
As this happens with probability m x i m y j , we can see that our beliefs about z can be represented by the mass function m z , defined as follows:

m z ([z]) = {i,j|[x i ]⋄[y j ]=[z]} m x i m y j , (34)
for all [z] ∈ IR, which generalizes (3). This approach can easily be extended to the case of a function of n variables. Let f be a function from R n to R, z = f (x 1 , . . . , x n ), and [f ] an inclusion function for f . Assume that our knowledge of each variable x i is represented by a mass function m x i . If each m x i ([x i ]) = 1 for some interval [x i ], then our knowledge of z can be represented by:

[z] = [f ]([x 1 ], . . . , [x n ]). (35)
In the general case, using the same line of reasoning as above, we can derive a mass function m z on z under assumption of independence, using the following formula:

m z ([z]) = {i 1 ,...,in|[z]=[f ]([x i 1 ],...,[x in ])} m x 1 i 1 . . . m xn in . ( 36
)
Example 2 Assume that m x and m y have the following expressions:

m x ([1, 2]) = 0.7, m x ([0, 3]) = 0.3, (37)
m y ([0, 1]) = 0.6, m y ([0, 2]) = 0.4. ( 38
)
Let z be a variable related to x and y by the equation z = x exp(y). What are our beliefs regarding z? Using [START_REF] Regan | Equivalence of methods for uncertainty propagation of real-valued random variables[END_REF], a mass function m z on z is then derived as

m z ([1, 2] • [exp]([0, 1])) = m z ([1, 2e]) = 0.7 × 0.6 = 0.42 (39) m z ([1, 2] • [exp]([0, 2])) = m z ([1, 2e 2 ]) = 0.7 × 0.4 = 0.28 (40) m z ([0, 3] • [exp]([0, 1])) = m z ([0, 3e]) = 0.3 × 0.6 = 0.18 (41) m z ([0, 3] • [exp]([0, 2])) = m z ([0, 3e 2 ]) = 0.3 × 0.4 = 0.12. ( 42
)

Summarization

Despite its success as a model of human reasoning under uncertainty, one of the arguments often raised against the theory of belief functions is its relatively high computational complexity. Indeed, the number of focal sets obtained when computing a function of n mass functions using (36) grows exponentially with n, which may lead to serious problems when both the number of steps in the computations and the size of focal sets become large [START_REF] Denoeux | Inner and outer approximation of belief structures using a hierarchical clustering approach[END_REF].

A natural way to decrease this complexity is therefore to decrease the number of focal sets by grouping similar or unimportant ones. This strategy is at the basis of the summarization method [START_REF] Lowrance | A framework for evidentialreasoning systems[END_REF], in which the focal sets with the smallest mass are aggregated and the sum of their masses transferred to their union.

More precisely, let m be a mass function on Ω with p focal sets A 1 , . . . , A p , p > 1. We note, as before, m i = m(A i ), i = 1, . . . , p. Assume that the focal sets are ordered such that m 1 ≥ m 2 ≥ . . . ≥ m p , and let 1 ≤ q < p. A mass function m ′ with q focal sets can be constructed from m as follows:

m ′ (A i ) = m i , i = 1, . . . , q -1, (43)
m ′ (A 0 ) = p i=q m i , (44)
where A 0 is defined as p i=q A i , or p i=q A i when the A i are intervals or boxes. The summarization method is a simple approximation method making it possible to control the complexity of computations in DS theory by bounding the number of focal sets. More sophisticated approximation schemes have been proposed in, for example, [START_REF] Denoeux | Inner and outer approximation of belief structures using a hierarchical clustering approach[END_REF] and [START_REF] Denoeux | Approximating the combination of belief functions using the fast moebius transform in a coarsened frame[END_REF]. Such refined schemes allow a better approximation at the cost of higher computing cost. The simpler and faster summarization method will be adopted in this paper.

Expectations

Let us assume in this subsection that Ω = R, and our state of knowledge regarding some variable x is represented by a mass function m. How can the concept of the expected value of x be defined in this context? An answer to this question is more easily found by interpreting m as defining a set of probability measures P such that bel(A) ≤ P (A) ≤ pl(A), [START_REF] Yager | Cumulative distribution functions from Dempster-Shafer belief structures[END_REF] for any measurable subset A of R. A probability measure P verifying the above inequalities is said to be compatible with m; the set of such probability measures will be noted as P(m).

The lower and upper expectations [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] of x with respect to m may then be defined, respectively, as

E * (m) = inf P ∈P(m) E(P ), (46)
E * (m) = sup P ∈P(m) E(P ). ( 47
)
When the focal sets are real intervals [x i ], then the lower-upper expectation interval [E](m) = [E * (m), E * (m)] (hereafter referred to as the interval expectation of m) may be shown [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] to be equal to

[E](m) = p i=1 m i • [x i ]. (48)
A particular value in this interval may be selected by averaging the center c i of each focal interval [x i ], which will be denoted:

E(m) = p i=1 m i c i . (49)
The quantity E(m) happens to be the expectation with respect to the pignistic probability measure associated with m [START_REF] Ph | The Transferable Belief Model[END_REF][START_REF] Petit-Renaud | Nonparametric regression analysis of uncertain and imprecise data using belief functions[END_REF]; it will be referred to as the pignistic expectation of m. Equations ( 48) and (49) can easily be extended to the multidimensional case where Ω = R n . The interval expectation of m is then an n-dimensional box, and its pignistic expectation is a vector of R n .

Belief State Estimation

As mentioned in Section 2.3, this interval framework guarantees the validity of the calculations under the assumptions that the system equations are correct and measurement errors remain within given bounds. However, this approach may be very conservative, which may result in very imprecise estimates of state variables. As noted in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], "the major implementation problem of such an approach is to determine correctly the bounds of the noises. Indeed, if these bounds are underestimated, the contractor may lead to no solution. On the contrary, if the bounds are overestimated, the estimated boxes can be very large (the estimates are then very pessimistic)".

In this section, we propose to solve this problem by replacing intervals by general mass functions with focal intervals. The construction of mass functions using prior knowledge of measurement error will first be addressed in Section 4.1. The bounded error estimation method recalled in Section 2.3 will then be extended to the Dempster-Shafer framework in Subsection 4.2. A simple illustrative example is given in Appendix A.

Construction of Mass Functions

As recalled in Section 1, classical state estimation techniques such as the Kalman filter and its extensions assume process and observation noises to be normal. In contrast, the BEE method (Section 2.3) makes the much weaker assumption that perturbations are bounded and that nothing is known except the bounds. While the former approach is very restrictive, the latter may be argued to be too conservative -as noted at the beginning of this section -leading to overly imprecise estimates of the system state. In many cases, however, additional realistic assumptions about error distributions can be made: for instance, we can assume some measure of central tendency such as the mode, the mean or the median to be equal to zero. As will be shown below, the belief function formalism allows to us to refine the BEE approach by making use of such additional assumptions.

Let us consider a continuous real random variable X characterized by a probability distribution P X with known support [a, b] and mode c ∈ (a, b). Let π denote the triangular possibility distribution

π(x) =      x-a c-a if a ≤ x < c, b-x b-c if c ≤ x < b, 0 otherwise. ( 50
)
As shown in [START_REF] Baudrit | Practical representations of incomplete probabilistic knowledge[END_REF], the following inequalities hold for every measurable set A:

Π(A) ≤ P X (A) ≤ 1 -Π(A), (51)
where A denotes the complement of A, and Π is the possibility measure associated with π:

Π(A) = sup x∈A π(A). (52)
The possibility distribution π can thus be seen as an approximation of the set of probability measures with support [a, b] and mode c. Such a continuous possibility measure can be itself approximated by a discrete mass function m with p focal intervals defined as follows.

For any α ∈ (0, 1], the α-cut of π is the set of values x such that π(x) ≥ α. This is the interval [π α ] with bounds:

π α = a + α(c -a) (53)
π α = c + (1 -α)(b -c). (54
m([π 0 ]) = α 1 , (55) m([π α 1 ]) = α 2 -α 1 , (56) . . . ( 57
) m([π α p-1 ]) = 1 -α p-1 . ( 58
)
The principle of this construction is illustrated in Figure 5. In this figure, there are two α-cuts, [π α 1 ] and [π α 2 ], associated with α 1 and α 2 respectively. By considering these α-cuts and [π 0 ] = [a, b], we can construct a mass function with p = 3 focal

elements: m([π 0 ]) = α 1 , m([π α 1 ]) = α 2 -α 1 and m([π α 2 ]) = 1 -α 2 .
It is clear that m corresponds to a possibility distribution π ′ that approximates π, and such that π ′ ≥ π. Consequently, if bel and pl denote, respectively, the belief and plausibility functions associated with m, the following inequalities hold for every event A: bel

(A) ≤ P X (A) ≤ pl(A), (59)
which means that m can be considered an approximation of the set of probability measures with known support [a, b] and mode c. We may remark that, in the above formalism, all values outside the interval [a, b] are considered to be absolutely impossible. In many cases, such an assumption may be argued to be too strong. For instance, if X denotes a measurement error, the bounds a and b may be derived from knowledge of the measurement process, which may be valid most of the time except in rare situations (for instance, in case of sensor failure). To account for the fact that values outside [a, b] are not impossible but have very small plausibilities, we may discount m by a small discount rate ǫ [39, page 251]. Discounting m means transferring to the frame Ω a fraction ǫ of the mass m(A) given to all proper subsets A of Ω. The resulting mass function ǫ is thus defined as

ǫ m(A) = (1 -ǫ)m(A) if A ⊂ Ω, (1 -ǫ)m(Ω) + ǫ if A = Ω. ( 60
)
We note that ǫ bel(A) = (1ǫ)bel(A) ≤ bel(A) and ǫ pl(A) = (1ǫ)pl(A) + ǫ ≥ pl(A) for all strict subsets A of Ω. Consequently, discounting has the effect of enlarging the intervals [bel(A), pl(A)], and (59) remains valid.

Example 3 Let U be a real random variable modeling the observation noise in some measurement process. Assume that expectation of U is zero, and its standard deviation σ is given. In bounded error approaches, it is common to take the interval [-3σ, +3σ] as the support of U . If we assume that the mode of U is zero, then we can better approximate its probability distribution by a symmetrical triangular possibility distribution, which can itself be approximated by the following discrete mass function with p = 3 focal intervals:

m([-3σ, 3σ]) = 1/3, (61) m
([-2σ, 2σ]) = 1/3, (62) m
([-σ, σ]) = 1/3. ( 63
)
Obviously, a better approximation of the continuous triangular possibility distribution can be obtained by increasing the number p of focal intervals, at the expense of higher computational complexity. Discounting m at rate ǫ and assuming Ω = R yields the following discounted mass function:

ǫ m([-3σ, 3σ]) = (1 -ǫ)/3, ( 64
) ǫ m([-2σ, 2σ]) = (1 -ǫ)/3, ( 65
) ǫ m([-σ, σ]) = (1 -ǫ)/3 (66) ǫ m(R) = ǫ. (67)

Extending the BEE Approach

Let us consider again the dynamical system described by Equations ( 20)- [START_REF] Gning | Dynamic vehicle localization using constraints propagation techniques on intervals. a comparison with Kalman filtering[END_REF]. However, in contrast with the interval analysis performed in Section 2.3, we now assume the initial state x 0 as well as the input vector u k and the observation vector z k at each time step k to be described by mass functions with interval focal sets. To denote a mass function, the variable will be indicated as a superscript, and the time step as a subscript. For instance, m x k will denote the mass function regarding the state x k at time step k.

The BEE algorithm can be generalized using the extension of interval analysis to belief functions, as recalled from Section 3.2. We recall that Algorithm 1 computes a box estimate [x k+1 ] of the state at each time step k + 1 from the imprecise input [u k ] and the output [z k+1 ], using [START_REF] Jaulin | Applied Interval Analysis[END_REF]. The theory recalled from Section 3.2 allows us to extend function C from intervals to mass functions.

Let p x k , p u and p z denote, respectively, the number of focal sets of m x k , m u k , m z (without loss of generality, we assume the mass functions on u and z to be computed using the procedure described in Section 4.1, with a constant number of focal intervals). The i-th focal element of m x k will be noted as [x i k ], and similarly for the other mass functions. Using [START_REF] Regan | Equivalence of methods for uncertainty propagation of real-valued random variables[END_REF], Equation ( 26) can be generalized to

m x k+1 ([x]) = {i,j,ℓ|[x]=C(f ,g,[x i k ],[u j k ],[z ℓ k+1 ])} m x k ([x i k ]) • m u k ([u j k ]) • m z k+1 ([z ℓ k+1 ]). (68)
The detailed procedure is described in Algorithm 2. It can be summarized as follows. The algorithm depends on three parameters: the number of focal elements

Algorithm 2 Belief State Estimation (BSE) algorithm. Input: m x k , m u k , m z k+1 , p u , p z , q Output: m x k+1 , [ x k+1 ],
x k+1 1: Read the input u k and its error. Deduce m u k with p u focal elements 2: Read the output z k+1 and its error. Deduce m z k+1 with p z focal elements 3: r ← 0 4: for i = 1 to p x k do 5:

for j = 1 to p u do 6:

for ℓ = 1 to p z do 7:

r ← r + 1

8:

[x r k+1 ] ← C f , g, [x i k ], [u j k ], [z ℓ k+1 ] .
9:

m x k+1 ([x r k+1 ]) ← m x k ([x i k ]) • m u k ([u j k ]) • m z k+1 ([z ℓ k+1 ]) 10:
end for 11:

end for 12: end for 13: Normalize m x k+1 14: Summarize m x k+1 to keep at most q focal elements. 15: % Computation of interval and point estimates %

16: [ x k+1 ] ← [E](m x k+1 ) 17: x k+1 ← E(m x k+1 )
p u and p z for the input and the output, respectively, and the maximum number q of focal elements for m x x+1 used by the summarization algorithm recalled from Section 3.3. In steps 1 and 2, the mass functions on the input at time k and on the output at time k + 1 are constructed using the method described in Section 4. has been computed, it is normalized (step 13) and summarized (step 14) using the summarization algorithm described in Section 3.3 with a maximum number q of focal elements. Normalization consists in dividing each mass assigned to a non empty focal set by 1m x k+1 (∅) and setting m x k+1 (∅) to zero. Finally, interval and pignistic expectations are computed using (48) and (49).

A simple illustrative example is given in Appendix A. This algorithm, hereafter referred to as the Belief State Estimation (BSE) algorithm, thus amounts to running the BEE algorithm several times (Algorithm 1). The maximum number of runs of BEE at each time step is p x p z q: the choice of these three parameters is thus crucial for controlling the execution time of the algorithm. The sensitivity to these three parameters of the results will be studied in the experiments presented in Section 5.

As explained in Section 3.2, Equation (68) results from the assumption of independence between mass functions m x k ([

x i k ]), m u k ([u j k ]
) and m z k+1 , in the precise sense outlined in Section 3.2. This hypothesis is the counterpart of the statistical independence assumptions that are commonly accepted when models are designed with random perturbations.

As a final remark, we note that the above algorithm may fail if each call to the BEE procedure in step 8 returns the empty set, in which case the normalization operation performed in step 13 is undefined. Such a situation may arise if the input u k or the output z k+1 are outliers and are inconsistent with the current state estimate encoded in mass function m x k . To avoid such a situation, we may discount m u k and/or m x k+1 by a small rate ǫ using (60), making the estimation procedure robust against outliers.

Application

In this section, we apply the belief state estimation (BSE) algorithm introduced above to dynamic localization of a land vehicle. The application and the data are identical to those used in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF] and [START_REF] Abdallah | Box particle filtering for non linear state estimation using interval analysis[END_REF] 1 . We will first introduce the problem and describe the experimental settings and then present the simulation results.

Problem Description

As noted in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], dynamic localization of a land vehicle in outdoor environments is a key issue in many applications. It consists of computing, at each time step k, an estimation of the vehicle position based on sensor measurements. Usually, two kinds of sensors are used: dead reckoning sensors (such as odometers, gyrometers, etc.) and absolute sensors, such as telemeters, goniometers, vision or Global Positioning System (GPS), among others. A major issue is to exploit the complementarity and redundancy between these sensors in order to achieve higher precision, availability and integrity. For instance, the quality of GPS positioning depends on the configuration of visible satellites, and GPS information can even become unavailable because of masking effects occurring in forests, tunnels, cities, etc. Fusing GPS information with dead reckoning sensor measurements makes it possible to filter the GPS estimates, thus increasing the performance of the localizer. The conventions and notations used in this section are illustrated in Figure 6. The mobile frame origin M is chosen as the middle of the rear axle. Let (x k , y k ) be the position of the vehicle and θ k its heading angle at time step k. Let x k = (x k , y k , θ k ) T be the state of the vehicle. As shown in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF], its evolution can be described by the following state equations:

     x k+1 = x k + δ S,k cos(θ k + δ θ,k
2 )

y k+1 = y k + δ S,k sin(θ k + δ θ,k 2 ) θ k+1 = θ k + δ θ,k , (69)
where δ S,k and δ θ,k denote, respectively, the elementary displacement and rotation of the vehicle at time step k. These quantities can be obtained with good precision with two rear-wheel ABS sensors and a fiber optic gyrometer, respectively. They are considered as the input to the system. With previous notations, we can thus denote u k = (δ S,k , δ θ,k ) T . The measure of the position z k = (x GP S , y GP S ) is obtained by a GPS receiver, after converting each estimated point (longitude, latitude) in a Cartesian local frame [START_REF] Gning | Dynamic vehicle localization using constraints propagation techniques on intervals. a comparison with Kalman filtering[END_REF]. We note that the heading angle θ k is not observed.

Experimental Settings

The results reported here were obtained using real sensor measurements collected using an experimental car (see [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF] for details about the sensors and data processing). The vehicle was equipped with two GPS receivers: a Trimble AgGPS 132 to estimate its position in real time, and a high-precision Thales Navigation GPS receiver to compute estimation errors. The Thales receiver was used in post-processed kinematic (PPK) mode with a local base (a Trimble 7400). This system was able to provide reference positions at a sampling rate of 1 Hz. Because the constellation of the satellites was satisfactory during all the trials, kinematic ambiguities were resolved, and an accuracy of a few centimeter was achieved. To avoid confusion, this reference data will be referred to as "PPK", while the term "GPS" will be reserved for data from the Trimble AgGPS 132 GPS receiver.

Experiments were carried out on a test track in Versailles (France). We report hereafter the analysis of a 4.7 km trajectory with a mean speed of 50 km/h (see the left-hand side of Figure 7). The zoomed part of the trajectory (right-hand side of Figure 7) shows the GPS and estimated positions.

As in [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF] and [START_REF] Abdallah | Box particle filtering for non linear state estimation using interval analysis[END_REF], the GPS error bounds were taken to be plus or minus three times the estimated standard deviation computed in real time by the receiver. Assuming the mode on the GPS error to be zero, mass functions m x 0 and m z k were constructed to approximate triangular possibility distributions, as explained in Section 4.1. The heading angle θ was initialized as [θ 0 ] = [0, 2π]. The error bounds of the input vector u k were computed from characteristics of the ABS sensor and the gyrometer (±3 degrees for δ θ,k ).

The three parameters of the BSE algorithm are p z , p u and q. In the results presented in the next section, the default values of these parameters are p z = 4, p u = 1 and q = 20. A sensitivity analysis of the influence of these parameters will also be presented.

Results

Figure 8 shows the focal sets of m x k at time steps k = 10 and 100. In this figure, the vehicle trajectory estimated by the GPS alone is plotted as a solid line and the estimated positions given by the BSE method are plotted as black ( * ) points. The interval expectation [E](m x k ) is represented in bold. For the sake of clarity, only the focal sets with mass greater than 0.05 are shown. The summarization process is illustrated by Figure 9, which shows the focal sets of m x k before and after summarization. Figures 10 and11 show the interval errors on x and y for the GPS as well as the BEE and BSE methods. We can see that the BSE method provides narrower intervals that still contain the true positions along both coordinates. The BSE method is also much less affected than is the BEE method by degradation of the GPS signal due to masking effects, as occurred around time steps 150, 210 and 400.

Figures 12 and13 show the distance (in meters) between estimated and actual positions of the vehicle for the following estimation methods: GPS, BEE and BSE (Figure 12); GPS, particle filter (PF) and BSE (Figure 13). Implementation and parameterization of the PF method were the same as reported in [START_REF] Abdallah | Box particle filtering for non linear state estimation using interval analysis[END_REF]. The PF was run with 1000 particles. We can see that the BSE method is clearly much more accurate than BEE. The difference is smaller between PF and BSE, with a small advantage exhibited by BSE most of the time.

These results are confirmed by Table 1, which reports the mean squared errors (MSE) on both coordinates x and y as well as the mean one-step running time (on a PC with Matlab) for the GPS, the PF (with 1000, 2000 and 3000 particles), bounded error estimator (BEE) and the belief state estimator (BSE) with the default parameter values (p z = 4, p u = 1 and q = 20). We can see that the BSE method significantly outperforms the BEE method, and is also slightly more accurate than the PF method, with comparable running time. Increasing the number of particles in the PF improves its performance, at a higher computational cost. Table 2 displays the average one-step-ahead running times as well as the mean squared errors on both spatial coordinates for the BSE algorithm with different choices for the three parameters p z , p u , q. The setting p z = p u = q = 1 corresponds to the BEE method. Parameter q seems to have a critical impact on both computational cost and performance. However, increasing it beyond 20 does not significantly improve estimation accuracy.

We can remark that simulations were done in Matlab on a PC computer and that the code was not optimized. Consequently, the reported running times in Tables 1 and 2 can only be used to compare the methods. We can conclude from these first results that our approach is no more costly, but slightly more accurate than the PF approach for the kind of application considered here. Although more costly than the BEE, our algorithm remains compatible with a real-time implementation.

Conclusion

This paper has proposed a new approach to non-linear state estimation based on belief function theory and interval analysis. This method uses mass functions composed of a finite number of axis-aligned boxes with associated masses. Such mass functions can encode partial information for model and measurement uncertainties more accurately than the bounded error approach alone. Focal sets are propagated in the system equations using interval arithmetics and constraint satisfaction techniques, thus generalizing pure interval analysis. Applied to state estimation in dynamical systems, this approach makes it possible to compute a belief function for the system state each time step.

This method was applied to the localization of a land vehicle based on the dynamic fusion of GPS measurements with ABS sensors and a gyrometer. The method was shown to provide more accurate estimates of vehicle position than did the bounded error method, while retaining the essential property of interval analysis of providing guaranteed computations (as the box provided by the BEE method is always one of the focal sets computed by the BSE method). The performance of the BSE approach is also slightly better than that of the particle filter method, with comparable running time. This suggests that our method might be a viable alternative to both bounded error and probabilistic Monte-Carlo approaches to non-linear state estimation, at least for the kind of applications considered here. The preliminary results reported in this paper will obviously have to be confirmed by further experiments related to more complex estimation problems.

A further advantage of the belief function approach is that it makes possible combining computed belief structures with additional information, such as digital road network data, which can also be represented in the belief function framework. A belief-function based map matching method that constructs a belief structure from a set of candidate roads extracted from an existing local map has already been presented in [START_REF] Nassreddine | Map matching algorithm using belief function theory[END_REF]. Further results in this direction will be reported in future publications.

  [x] • [y] = [min(xy, xy, xy, xy), max(xy, xy, xy, xy)].

Figure 1 :

 1 Figure 1: Images of a box [x] by a function f , a pessimistic inclusion function [f ] and the minimal inclusion function [f ] * .

Figure 2 :

 2 Figure 2: A CSP in R 2 .

Figure 3 :

 3 Figure 3: Graph corresponding to the CSP of Example 1.

Figure 4 :

 4 Figure 4: Forward propagation in the BEE algorithm. A set estimate of x k+1 is computed using the current set estimates of x k and u k . The "predicted output" [g]([x k+1 ]) is then intersected with the measured output, which allows us to contract [z k+1 ]. The constraints are then back-propagated to contract [x k+1 ],[x k ] and [u k ].

)

  By convention, let [π0 ] = [a, b]. Let us consider p -1 distinct values for α: 0 = α 0 < α 1 < α 2 < . . . < α p-1 < 1,and let m be the consonant mass function with focal intervals [π α k ], k = 0, . . . , p -1, defined by:

Figure 5 :

 5 Figure 5: Approximation of a triangular possibility distribution π by a discrete consonant mass function. The focal set of m are defined as the support of [π 0 ] of π, and α-cuts [π α k ] for k = 1, . . . , p -1 (p = 3 in this figure). The resulting possibility distribution π ′ dominates π.

  1. The BEE algorithm (Algorithm 1) is then applied for each choice of focal element of [x i k ] of m x k , focal interval [u j k ] of m u k and focal element [z ℓ k+1 ] of m z k+1 (steps 4-12). The belief mass assigned to the r-th focal element [x r k+1 ] of m x k+1 is the product of the masses assigned to [x i k ], [u j k ] and [z ℓ k+1 ] (line 9). Once the whole mass function m x k+1

Figure 6 :

 6 Figure 6: Vehicle representation at time steps k and k + 1.

Figure 7 :

 7 Figure 7: Test trajectory (left), with a zoom on part of the trajectory (right) showing the GPS positions (solid line) and the estimated positions (*).

Figure 8 :

 8 Figure 8: Focal sets of the mass function for x k , with mass greater than 0.05, at time steps k = 10 (up) and k = 100 (down). The GPS and estimated positions are shown as a solid lines and * symbols, respectively. The interval expectation [E](m x k ) is represented in bold. The vehicle is moving downwards.

Figure 9 :

 9 Figure 9: Focal sets of the mass function on x k before (left) and after (right) summarization. The line width is proportional to the mass.

Figure 10 :

 10 Figure 10: Interval errors on x for GPS (dashed line), BSE (bold line) and BEE (thin line).

Figure 11 :

 11 Figure 11: Interval errors on y for GPS (dashed line), BSE (bold line) and BEE (thin line).

Figure 12 :

 12 Figure 12: Distance in meters between the estimated and actual vehicle positions for the GPS (bold line), BEE (dotted line) and BSE (thin line).

Figure 13 :

 13 Figure 13: Distance in meters between the estimated and actual vehicle positions for the GPS (bold line), PF (dotted line) and BSE (thin line).

  u k and its error. Deduce [u k ]. 2: Read the output z k+1 and its error. Deduce [z k+1 ].

3: [x k+1 ] ← (-∞, +∞) nx 4: while at least one variable is contracted do 5:

Table 1 :

 1 Mean squared errors (MSE) on x and y, and one-step running times for the GPS without filtering, and the following methods: particle filter (PF) with 3000, 2000, 1000 particles, bounded error estimation (BEE) and belief state estimation (BSE) with q = 15 and 20.

		GPS PF(3000) PF(2000) PF(1000) BEE BSE (20)
	MSE on x (m 2 )	0.134	0.119	0.121	0.125	0.123	0.118
	MSE on y (m 2 )	0.374	0.215	0.232	0.243	0.249	0.199
	Running time (ms)	-	639	526	401	136	409

Table 2 :

 2 Mean squared errors (MSE) on x and y (in m 2 ), and average one-step running times (in ms) for the belief state estimation (BSE).

	p z p u q	running time MSE on x MSE on y
	1	1	1	136	0.123	0.249
	2	1 10	158	0.121	0.240
	3	1 10	198	0.120	0.236
	4	1 10	241	0.120	0.229
	2	2 10	199	0.121	0.236
	2	4 10	236	0.120	0.229
	2	6 10	316	0.120	0.221
	4	1 15	386	0.120	0.225
	4	1 20	409	0.118	0.199
	4	1 22	426	0.118	0.197
	4	1 25	452	0.116	0.197

We thank Philippe Bonnifait and his group for providing and allowing us to use these data.

A Illustrative Example

To illustrate the computations performed in the BSE algorithm, let us consider as a very simple illustrative example the linear motion of a vehicle on an axis x. The vehicle position x k at each time step k is assumed to be governed by the following equations:

where u k is the vehicle speed at time step k measured by an odometer sensor with an uncertainty of ±3 m/s, z k is the measured position of the vehicle given by a GPS receiver and T is the time sampling taken, equal to 1 s. At each time step k, we receive mass functions m u k on the vehicle speed and m z k on its position. Furthermore, an initial estimate of the position x 0 at time k = 0 is given in the form of a mass function m x 0 . Assume that m x 0 , m u 0 and m z 1 have the following expressions:

We will detail the calculations needed to compute the mass function m x 1 regarding the position x 1 of the vehicle at time step k = 1. Because m x 0 and m z 1 have two focal intervals each and m u has only one, the FBP algorithm has to be run four times, corresponding to the four different ways of choosing one focal interval for each of the three BSs. Calculations are detailed below. [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Moore | Methods and Applications of Interval Analysis[END_REF]. The forward steps of the FBP algorithm are as follows:

[

We can verify that this interval is no longer contracted by further steps of the FBP algorithm. The corresponding mass is:

Again, this interval is not contracted by further steps of the FBP algorithm. The corresponding mass is: [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Moore | Methods and Applications of Interval Analysis[END_REF]. The forward steps of the FBP algorithm are:

We thus have

and [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF]. The predicted state at k = 1 is the same as above: [x 1 ] = [START_REF] Denoeux | Approximating the combination of belief functions using the fast moebius transform in a coarsened frame[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF], and the intersection with [z 2 1 ] yields:

[z 1 ] ← [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF] ∩ [START_REF] Denoeux | Approximating the combination of belief functions using the fast moebius transform in a coarsened frame[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF] = [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF]. (87)

We thus have [x 4 1 ] = [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF]. Again, this interval is not contracted by the FBP algorithm. The corresponding mass is:

We thus obtain a non normal mass function with 4 focal intervals. After normalization, we get: