- Email: [email protected]
Continuous–discrete confidence interval observer – Application to vehicle positioning
Information Fusion 14 (2013) 541–550
Contents lists available at SciVerse ScienceDirect
Information Fusion journal homepage: www.elsevier.com/locate/inffus
Continuous–discrete confidence interval observer – Application to vehicle positioning G. Goffaux, M. Remy, A. Vande Wouwer ⇑ Service d’Automatique, Université de Mons (UMONS), Mons, Belgium
a r t i c l e
i n f o
Article history: Received 27 November 2012 Received in revised form 12 January 2013 Accepted 13 February 2013 Available online 28 February 2013 Keywords: Intervals State observer Positioning systems Railways Vehicles
a b s t r a c t In vehicle positioning applications, the confidence level in the position and velocity estimates can be even more significant than accuracy. In this study, a probabilistic interval method is proposed, which combines, through union and intersection operations, the information from a possibly uncertain predictor (the vehicle model) and measurement sensors. The proposed method is compared to Kalman filtering and to guaranteed interval estimation in the context of railway vehicles where security is the key objective. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Position and velocity estimation is generally performed using data fusion algorithms which combine information from a mathematical model and sensor measurements. Popular estimation methods include Kalman filtering based on a kinematic model of the vehicle and measurements from GPS (Global Positioning System) and INS (Inertial Navigation System) [1], particle filtering for positioning, navigation and tracking [2] or nonlinear filtering in marine applications [3]. Accuracy and security are two conflicting objectives. However, security requirements are particularly important in vehicle positioning as collisions could cause human and economical loss. Hence, the estimated variables have to be defined by confidence intervals in applications such as safe railway positioning [4]. Recently, state estimation methods based on interval algebra [5] have been proposed to compute guaranteed intervals upper bounding the solution sets in the case of uncertain discrete-time models. Subsequent important contributions in this area of research include the development of methods for nonlinear continuous-time systems, see for instance [6–10]. These developments have been of particular interest in the study of biological systems where the models are nonlinear and uncertain [11–14]. However, these latter results do not give a realistic measure of the confidence level as they consider guaranteed intervals. Other studies have preferred a statistical approach for tackling the prob⇑ Corresponding author. Tel.: +32 (0)65374141; fax: +32 (0)65374136. E-mail address: [email protected] (A. Vande Wouwer). 1566-2535/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.inffus.2013.02.006
lem of continuous–discrete estimation with faulty sensors [15], either by considering the manipulation of confidence intervals in a static configuration [16,17], or uncertainties in the probabilities (for example, a variance interval) [18,19]. In this study, a continuous–discrete confidence interval estimator is proposed, which proceeds in two steps. The first step uses a continuous-time predictor to compute the evolution of lower and upper bounds for each state variable between two discrete measurement times. The second step occurs at a measurement instant, where intervals from the predictor and from hardware sensors are combined using union and intersection operations so as to achieve specified confidence levels. In [20], a preliminary version of the algorithm has been designed in the context of a marine vehicle described by a nonlinear dynamic model. The effect of state transformations on the nonlinear predictor has been studied and a security measure has been defined which corresponds to a given level of error tolerance within a given time span. In this study, the method is further developed in the context of a railway application where the security demands are more severe. The prediction step uses a dynamic model based on the equilibrium of forces. In addition, several sensors are considered, providing position, velocity or acceleration information. Moreover, a comparison with two well-established techniques, Kalman filtering and guaranteed interval estimation, is achieved so as to highlight the advantages of the proposed probabilistic interval approach. In the sequel, the paper is organized in four main sections. In Section 2, some preliminary model assumptions are introduced
542
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
along with the concept of confidence intervals. Section 3 is devoted to the general structure of the continuous–discrete confidence interval observer. In Section 4, the observer is applied to a test example, the positioning of a railway vehicle. After a description of the vehicle model, a dedicated interval observer is proposed. In Section 5, the probabilistic interval observer is compared with an extended Kalman filter (EKF) and a guaranteed interval observer.
Consider the vehicle-sensors system described by the following equations for t P t0:
ð RÞ :
To compute the resulting confidence levels, the intervals are associated with Boolean random variables denoted B which are defined as follows [21]:
Bx ¼ true if x 2 ½x ; xþ x; x ; xþ 2 R;
x < xþ
Bx ¼ false otherwise
ð7Þ
Consequently, the probability of enclosing a variable x by an interval is related to the corresponding Boolean random variable:
2. Preliminaries
(
2.1. Intersection and union operations
_ xðtÞ ¼ fðxðtÞ; uðtÞ; hÞ xðt0 Þ ¼ x0
ð1Þ
yk ¼ yðt k Þ ¼ Cðt k Þxðt k Þ þ ðt k Þ
The first part of (R) models the time evolution of positioning variables and consists of (possibly nonlinear) differential equations. xðtÞ 2 X Rnx is the vector of state variables (including position, velocity and acceleration) with X a compact set of Rnx . x0 is the initial state vector. In this study, prior confidence intervals are defined by the probability (denoted P) that the corresponding bounds include the unknown values, at the initial time with confidence level b0.
P x0 2 x0 ; xþ0 P b0 () h i h i P x0;1 2 x0;1 ; xþ0;1 ; . . . ; x0;nx 2 x0;nx ; xþ0;nx P b0
ð8Þ
From this information, union and intersection operations can be defined in terms of Boolean random variables [21,22]. We consider two measuring a variable x with the intervals independent þsensors x1 ; xþ 1 and x2 ; x2 . Bx1 and Bx2 are the corresponding Boolean random variables and bx1 and bx2 are the lower bounds on the integrity levels.
PðBx1 ¼ trueÞ ¼ Pðx 2 x1 ; xþ1 Þ P bx1 þ PðBx2 ¼ trueÞ ¼ Pðx 2 x2 ; x2 Þ P bx2 PðBx1 ¼ falseÞ ¼ Pðx R x1 ; xþ1 Þ 6 1 bx1 PðBx2 ¼ falseÞ ¼ Pðx R x2 ; xþ2 Þ 6 1 bx2
ð2Þ Intersection:
nu
uðtÞ 2 U R is the vector of system inputs including, for instance, forces controlling the vehicle motion. U is the set of admissible input signals, a subset of the space of measurable bounded piecewise constant functions. Consequently, utk !t is a piecewise constant function between tk and t which can be determined by the knowledge of a finite number of input values u(tk), . . . , u(t). In this study, they are known with some uncertainties such that confidence intervals are defined by
h i P utk !t 2 utk !t ; uþtk !t ()
PðBx ¼ trueÞ ¼ Pðx 2 ½x ; xþ Þ P bx
P x 2 x1 ; xþ1 \ x2 ;xþ2 ¼ PðBx1 ¼ true and Bx2 ¼ trueÞ P bx1 bx2
ð9Þ
Union: P x 2 x1 ;xþ1 [ x2 ;xþ2 ¼ PðBx1 ¼ true or Bx2 ¼ trueÞ ¼ 1 PðBx1 ¼ false and Bx2 ¼ falseÞ P 1 ð1 bx1 Þð1 bx2 Þ ¼ bx1 þ bx2 bx1 bx2
ð10Þ
Pðuðt k Þ 2 ½u ðt k Þ; u ðt k Þ; . . . ; uðtÞ 2 ½u ðtÞ; u ðtÞÞ P buk ðtÞ
Remark 2.1. If x is a nx-dimension vector, its confidence level is defined by the joint probability that the intervals related to its components include the unknown values. x ¼ ½x1 ; . . . ; xnx T
Moreover, h 2 H 2 Rnh is a model parameter vector and again confidence intervals are given (e.g. by a parameter identification procedure).
Pðx 2 ½x ; xþ Þ ¼ P x1 2 x1 ; xþ1 ; . . . ; xnx 2 xnx ; xþnx
Pðh 2 ½h ; hþ Þ P bh
Remark 2.2. If there is a dependence between two intervals, (9) and (10) become
þ
ð3Þ
þ
ð4Þ
The second part of (R) describes the measurement process. In vehicle positioning, sensors include accelerometers, gyroscopes, radars, GNSS (Global Navigation Satellite System) receivers providing partial information on vehicle position, velocity and acceleration. yðtk Þ 2 Rny is the vector of sampled measurements. Usually, measurements are modelled using additive random errors (tk) with Gaussian distribution, but other situations can be considered as well. Confidence intervals are given by
P yk 2 yk ; yþk P byk
P x 2 x1 ; xþ1 \ x2 ;xþ2 ¼ PðBx1 ¼ true and Bx2 ¼ trueÞ ¼ PðBx1 ¼ trueÞ þ PðBx2 ¼ trueÞ PðBx1 ¼ true or Bx2 ¼ trueÞ PPðBx1 ¼ trueÞ PðBx2 ¼ falseÞ Pbx1 þ bx2 1
ð5Þ
ð11Þ
Union:
Typically, each sensor has its own sampling period and measurements are asynchronous. Thus, tk is a measurement instant where C(tk) is the observation matrix corresponding to the active sensors. From these definitions, a confidence interval observer is designed in order to derive bounds, x(t) and x+(t), with a specified confidence level bobj (t P t0):
PðxðtÞ 2 ½x ðtÞ; xþ ðtÞÞ P bobj
Intersection:
ð6Þ
P x 2 x1 ; xþ1 [ x2 ;xþ2 ¼ PðBx1 ¼ true or Bx2 ¼ trueÞ P maxðbx1 ; bx2 Þ
ð12Þ
Remark 2.3. If the intervals are disconnected, the union operation gives the convex hull. For instance, [2,3] and [4,5] gives [2,5]. In the sequel, these operations are used to develop a confidence interval observer.
543
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
3. Continuous–discrete confidence interval observers
x_ i ðt k Þ 6
The algorithm introduced in [20] proceeds in two steps. Between two measurements times, the state intervals are propagated (or predicted) using the model equations. Then, the predicted intervals and the measurement intervals are combined to produce corrected intervals.
x_ þi ðt k Þ P
(
Ti ðtk Þ :
x_ P ðtÞ ¼ f xP ðtÞ; xþP ðtÞ; u ðtÞ; uþ ðtÞ; h ; hþ þ x_ þP ðtÞ ¼ f xP ðtÞ; xþP ðtÞ; u ðtÞ; uþ ðtÞ; h ; hþ
ð13Þ
xþP ðt k Þ ¼ xþC ðtk Þ
where f and f+ are derived from the model equations, and depend on bounds on the state variables, inputs and parameters. Before computing the resulting confidence level, the predictor equations have to satisfy the following property. Property 3.1 (Enclosure of the prediction equations). If + þ þ xðtk Þ 2 x P ðt k Þ; xP ðt k Þ ; uðtÞ 2 ½u ðtÞ; u ðtÞ and h 2 [h ,h ], then the solution of the prediction equations f and f+ (given by (13)) encloses þ the state variables if it satisfies xðtÞ 2 x P ðtÞ; xP ðtÞ 8t P t k . This property is related to the definition of an error vector eðtÞ 2 R2nx which has to remain nonnegative.
eðtÞ ¼
þ
xP ðtÞ xðtÞ eþ ðtÞ ¼ P 0 8t P t k xðtÞ xP ðtÞ e ðtÞ
ð14Þ
Proof. It is sufficient to check the following condition, assuming e(tk) > 0
If 9t P tk jei ðt Þ ¼ 0;
ej ðt Þ P 0 8j – i 2 f1; . . . ; 2nx g; then e_ i ðt Þ P 0
ð15Þ
Each time a component ei of the error vector vanishes, its derivative must be nonnegative in order to prevent a negative error [20]. f and f+ are derived from the evolution equation f using Müller’s theorem [23,7]. h T
Theorem 3.1. Assume that the vector of functions f ¼ ½f1 ; . . . ; fnx in (1) is continuous on a domain
8 > < t0 6 t 6 tf T : x ðtÞ 6 xðtÞ 6 xþ ðtÞ; > : h 2 H; uðtÞ 2 U
x ðtÞ; xðtÞ; xþ ðtÞ 2 X
T þ where x ðtÞ ¼ and x ðtÞ ¼ xþ are 1 ðtÞ; . . . ; xnx ðtÞ continuous on [t0, tf] and such that
T x 1 ðtÞ; . . . ; xnx ðtÞ
Tþi ðtk Þ :
þ
þ þ 1. x ðt0 Þ ¼ x 0 and x ðt 0 Þ ¼ x0 2. the derivatives of x(tk) and x+(tk) at tk (t0 6 tk 6 tf) satisfy for i = 1, . . . , nx,
max fi ðxðt k Þ; uðt k Þ; hÞ
ð17Þ
xðtk Þ2Tþ ðt k Þ i
8 t ¼ tk > > > > > > > < xi ðtÞ ¼ xi ðtÞ > > xj ðtÞ 6 xj ðtÞ 6 xþj ðtÞ; > > > > > : h 2 H; uðtÞ 2 U
j–i
8 t ¼ tk > > > > > < xi ðtÞ ¼ xþ ðtÞ i > xj ðtÞ 6 xj ðtÞ 6 xþj ðtÞ; > > > > : h 2 H; uðtÞ 2 U
j–i
þ Then, for any x0 2 x 0 ; x0 , a solution to (1) exists, which remains in
E:
xP ðt k Þ ¼ xC ðtk Þ
ð16Þ
þ where T i ðt k Þ and Ti ðt k Þ are subsets of T respectively defined by:
3.1. Prediction step The prediction step aims at propagating the interval bounds for each state variable from a start time tk using the information from the model equations only. For the construction of the confidence interval predictor, it is first assumed that guaranteed intervals exist for the initial state vector, the inputs and the parameters. Advanced methods for the design of guaranteed interval predictors can be employed [8–10]. Subsequently, this assumption is relaxed considering confidence levels for the variables and parameters. The typical form of an interval predictor consists of two coupled dynamic systems (t P tk) which are initialized with intervals ob þ tained in the previous correction step: x C ðt k Þ; xC ðt k Þ .
min fi ðxðt k Þ; uðt k Þ; hÞ
xðt k Þ2T ðtk Þ i
t0 6 t 6 tf
ð18Þ
x ðtÞ 6 xðtÞ 6 xþ ðtÞ
Proof. To prove relation (18), condition (15) is checked. By definiþ tion of T i ðt k Þ and Ti ðt k Þ, condition (15) is satisfied. For instance, at ðt Þ are the true state vector and the true ðt Þ and u time t = t⁄, if x input respectively along with x+(t⁄) the estimated upper bounds, one has:
if eþi ðt Þ ¼ xþi ðt Þ xi ðt Þ ¼ 0 and xj ðt Þ 6 xj ðt Þ 6 xþj ðt Þ ðj – iÞ
e_ þi ðt Þ P
ðt Þ; hÞ P 0 ðt Þ; u max fi ðxðt Þ; uðt Þ; hÞ fi ðx
xðt Þ2Tþ ðt Þ i
with
8 t ¼ t > > > < xi ðtÞ ¼ xþi ðtÞ ¼ xi ðtÞ Tþi ðt Þ : þ > > xj ðtÞ 6 xj ðtÞ 6 xj ðtÞ; > : h 2 H; uðtÞ 2 U
j–i
Theorem 3.1 can be extended to uncertain parameters and inputs. Moreover, a specific version of Theorem 3.1 can be obtained when f satisfies monotonicity properties [24]. Finally, advanced guaranteed interval predictors could be easily applied [13,14,10], if it is required by the complexity of the system model. Remark 3.1. Interval predictors can be designed from additional models such as kinematic and dynamic models of a vehicle. The resulting intervals can be combined if each of the interval predictor satisfies Property 3.1. Extending the previous theorem to confidence intervals, the following corollary gives a lower bound of the confidence level related to the intervals computed by an interval predictor [20]. Corollary 3.1. If confidence intervals are known at time tk (related to h i þ þ and the the state vector x C ðt k Þ; xC ðt k Þ , the inputs ut k !t ; ut k !t parameters [h,h+]) and if the prediction equations (f and f+) satisfy Property 3.1 of enclosure, the intervals on the predicted state þ variables, computed by the interval predictor x P ðtÞ; xP ðtÞ , satisfy the following lower bound of the integrity level [25] (tk < t 6 tk+1)
544
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
P xðtÞ 2 xP ðtÞ; xþP ðtÞ P P xðt k Þ 2 xC ðtk Þ; xþC ðtk Þ ; utk !t h i 2 utk !t ; uþtk !t ; h 2 ½h ; hþ P bxk buk ðtÞbh
ð19Þ
with bxk ; bh and buk the integrity related, respectively, to the state vector at tk, the parameters and the inputs (independence is assumed). Moreover, the inequality is also satisfied for each component of the state vector. Indeed,
h i P xi ðtÞ 2 xi;P ðtÞ; xþi;P ðtÞ P 8i 2 1; . . . ; nx h i P xðtk Þ 2 xC ðt k Þ; xþC ðt k Þ ; utk !t 2 utk !t ; uþtk !t ; h 2 ½h ; hþ P
3 intervals (i = 3): there is only one selection (j = 1) of intervals (3 elements out of a set of 3 elements) but there are four possible sequences of operations (from k = 1 to k = 4): {\, \}, {[, [}, {[, \}, {\, [} and also three possible orderings of the intervals when applying this sequence of operations: {1, 2, 3}, {1, 3, 2}, {2, 3, 1}. Consequently, a fourth dimension is added to the notation: c(i, j, k, l) which represents the possible orderings of the selected intervals. Hence, one has the following operations (the operations are performed from left to right):
cð3; 1; 1Þ ¼ ½x1 \ ½x2 \ ½x3 ¼ ½5:3; 13:6 b P 1 3 103
ð20Þ
cð3; 1; 2Þ ¼ ½x1 [ ½x2 [ ½x3 ¼ ½3:9; 15:4 b ¼ 1 109
bxk buk ðtÞbh
cð3; 1; 3; 1Þ ¼ ½x1 \ ½x2 [ ½x3 ¼ ½3:9; 15:4 b P 1 2 106
Proof. It is a consequence of Property 3.1. Indeed, f and f+ are constructed in such a way that the computed bounds x P ðtÞ and xþ P ðtÞ enclose the unknown state trajectory provided that the inputs, the parameter and the initial state vector are within their respective interval. h
cð3; 1; 3; 2Þ ¼ ½x1 \ ½x3 [ ½x2 ¼ ½4:4; 14:1 b P 1 2 106 cð3; 1; 3; 3Þ ¼ ½x2 \ ½x3 [ ½x1 ¼ ½4:4; 14:1 b P 1 2 106 cð3;1; 4; 1Þ ¼ ½x1 [ ½x2 \ ½x3 ¼ ½4:4; 14:1 b P 1 103 106 cð3;1; 4; 2Þ ¼ ½x1 [ ½x3 \ ½x2 ¼ ½4:4; 13:6 b P 1 103 106
3.2. Correction step
cð3;1; 4; 3Þ ¼ ½x2 [ ½x3 \ ½x1 ¼ ½5:3; 14:1 b P 1 103 106
The correction step proceeds by union and intersection operations between intervals from hardware sensors and intervals given by the predictor at the considered measurement time. A lower bound of the resulting confidence levels is computed by (9) and (10) in case of independence. For each state variable, intervals are built considering each possible combination (union and intersection) among the available confidence intervals representing the state variable. Then, the best resulting interval is derived selecting the smallest nonzero interval satisfying the integrity objective. To illustrate the methodology, an example is considered with 3 intervals, [x1] = [5.3, 14.1], [x2] = [4.4, 13.6] and [x3] = [3.9, 15.4]. Each of them has a confidence level of 1–103. At first, there are 3 possible choices involving 1, 2 or 3 intervals. 1 interval: the set of combined intervals, denoted c, only contains the original intervals with their corresponding confidence level. The notation is as follows: c(i, j) has two dimensions, i corresponding to the number of considered intervals (here i = 1) and j to the selected interval (here, j = 1, 2, or 3). 3
cð1; 1Þ ¼ ½x1 ¼ ½5:3; 14:1 b ¼ 1 10 cð1; 2Þ ¼ ½x2 ¼ ½4:4; 13:6 b ¼ 1 103 cð1; 3Þ ¼ ½x3 ¼ ½3:9; 15:4 b ¼ 1 103 2 intervals: with two intervals (i = 2), there are two possible operations, either an intersection or a union. The number of paired intervals is given by the number of combinations of 2 elements in a set of three elements. Now, c(i, j, k) has three dimensions. j indicates the paired interval and k the selected operations: intersection (k = 1) or union (k = 2). Consequently, one has:
cð2; 1; 1Þ ¼ ½x1 \ ½x2 ¼ ½5:3; 13:6 cð2; 1; 2Þ ¼ ½x1 [ ½x2 ¼ ½4:4; 14:1 cð2; 2; 1Þ ¼ ½x1 \ ½x3 ¼ ½5:3; 14:1 cð2; 2; 2Þ ¼ ½x1 [ ½x3 ¼ ½3:9; 15:4 cð2; 3; 1Þ ¼ ½x2 \ ½x3 ¼ ½4:4; 13:6 cð2; 3; 2Þ ¼ ½x2 [ ½x3 ¼ ½3:9; 15:4
b P 1 2 103 b ¼ 1 106 b P 1 2 103 b ¼ 1 106 b P 1 2 103 b ¼ 1 106
The resulting confidence level of an intersection is obtained from (9): b = (1–103)2 P 1–2 103. The resulting confidence level of a union is obtained from (10): b = 2(1–103)– (1–103)2 = 1–106. Also, with more than 2 intervals, the order of the operations has an importance, as it will be seen in the sequel.
The resulting confidence level of two consecutive intersections is obtained by: b = (1–103)3 P 1–3 103 whereas for two consecutive unions: b = (1–103) + (1–106)–(1– 103)(1–106) = 1–109. Also, one can note that consecutive unions or intersections are independent of the order of the intervals. The resulting confidence level of the sequence {\, [} is obtained by: b = (1–103)2 + (1–103)–(1–103)2(1–103) = 1– 2 106 + 109 P 1–2 106 whereas for the sequence {[, \}: b = (1–106)(1–103) P 1–103–106. To conclude the example, if the integrity objective is 1–106, the smallest interval satisfying the objective is [x1] [ [x2] = [4.4, 14.1]. A similar procedure can be performed for a larger number of intervals. Usually, the number of sensors is limited and combinations of 3 or 4 intervals are sufficient. 4. Vehicle positioning: a railway vehicle example In this section, the model of the illustrative example is introduced. Then, the confidence interval observer is developed and simulation tests illustrate the method potentials. 4.1. Railway vehicle: a dynamic model Considering the equilibrium of forces acting on the railway vehicle (Fig. 1) along the railroad and applying the fundamental principles of dynamics, one has [26,27]
M
dv ¼ F j R Mg sin a dt
ð21Þ
In this equation, v is the vehicle velocity and Fj is the force transmitted by the motor to the vehicle wheels. In this simplified model, no loss of adhesion is considered and consequently, the torque delivered by the motor is totally converted into force. R contains all the resisting forces opposing the vehicle movement. In general, these forces depend on the vehicle velocity and can be represented by a formula including an independent term, a term proportional to the velocity and a quadratic term:
R ¼ A þ Bv þ C v 2
ð22Þ
545
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
In (25), r(p) is the derivative of the slope angle with respect to the vehicle position. In the sequel, the state vector x is used to represent the vehicle position p, the vehicle velocity v and the slope angle a. Fig. 2 simulates the evolution equations during 600 s with the parameters given by Table 1. The slope changes with the vehicle position and the force Fj is computed so as to achieve a constant velocity. 4.2. Confidence interval observer Fig. 1. Forces applied to the vehicle.
Table 1 Model parameters (Eurostar train with 20 vehicles) [27]. Description
Parameter
Value
Unit
Vehicle mass Earth’s gravity
M g
816 10
t m/s2
Coefficients of resisting forces R
A B C
4.82 0.23508 0.013608
kN t/s t/m
The term A + Bv is related to mechanical and rolling frictions. A depends on the load applied to the wheels and B gives information on the quality of the railway track and on the vehicle stability. Finally, Cv2 represents the aerodynamic forces in the direction of the vehicle movement. C depends on the drag coefficient. The latter varies with the design of the vehicle profile. Moreover, the gravity force projected into the direction of motion is given by Mg sin a where a is the track slope varying according to the train position. Numerical values used in this simulation example are given in Table 1. The evolution equations can be summarized as follows:
p_ ¼ v 1 v_ ¼ ðF j A Bv C v 2 Þ g sin a M a_ ¼ rðpÞv
In the following, a confidence interval observer is designed based on the evolution Eqs. (23–25) and sensors measuring position, velocity and acceleration. At first, the interval predictor is proposed based on Müller’s theorem. Then, the correction step and the combination of confidence intervals are introduced. 4.2.1. Interval predictor Lower and upper bounds are computed by integrating a couple of differential equations related to the vehicle position, velocity and acceleration. The equations related to the vehicle position are obtained maximizing the vehicle velocity as a consequence of Müller’s theorem:
p_ P ðtÞ ¼ p_ þP ðtÞ ¼
ð25Þ
v ðtÞ ¼ v P ðtÞ
max
v ðtÞ ¼ v þP ðtÞ
v ðtÞ2½v P ðtÞ;v þP ðtÞ
ð26Þ
Moreover, assuming no uncertainties on the parameter values and the input Fj, applying Müller’s theorem gives the following derivative for the lower bound of the vehicle velocity.
v_ P ðtÞ ¼
ð23Þ ð24Þ
min
v ðtÞ2½v P ðtÞ;v þP ðtÞ
v_ P ðtÞ ¼
Fj A B C 2 min v ðtÞ M v ðtÞg sin aðtÞ pðtÞ 2 pP ðtÞ;pþP ðtÞ M M M v ðtÞ ¼ v ðtÞ P þ aðtÞ 2 aP ðtÞ; aP ðtÞ Fj A B C v ðtÞ v P 2 ðtÞ g sin aþP ðtÞ M M M M P
Similarly, one has
Fig. 2. Simulation of the dynamic model.
ð27Þ
546
v_ þP ðtÞ ¼
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
Fj A B C 2 v þ ðt Þ v þP ðtÞ g sin aP ðtÞ M M M M P
ð28Þ
In this formulation, a(t) lies between p2 and p2 such that sina is monotonic. In the same way, the derivatives of the bounds of the slope angle are computed:
a_ P ðtÞ ¼
min
v ðtÞ 2 v P ðtÞ; v þP ðtÞ
rðtÞ v ðtÞ
rðtÞ 2 ½rP ðtÞ; rþP ðtÞ a_ þP ðtÞ ¼
max rðtÞv ðtÞ v ðtÞ 2 v P ðtÞ; v þP ðtÞ rðtÞ 2 rP ðtÞ; rþP ðtÞ
ð29Þ
In (29), the optimum values are obtained from the product of two intervals. The first one comes from the vehicle velocity v(t) and the second one from the derivative r(t) of the slope angle with respect to the position p(t). The latter is obtained from the interval of the vehicle position considering the worst slope angle derivatives. For instance, r P ðtÞ ¼ minpðtÞ2½p ðtÞ;pþ ðtÞ rðpðtÞÞ. P
P
Remark 4.1. The model equations can be more complicated depending on the application and the desired precision. However, Müller’s theorem can still be applied in such cases. As an example, consider a part of the 2D evolution equations of a vessel [3]
p_ x ðtÞ ¼ v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
v1(t) and v2(t) are the velocity components given in a frame related to the vessel movement. /(t) is the vehicle yaw angle whereas px is the vessel position in the x-axis of the reference frame (see [3] for more details). In these equations, /(t) appears in a non-monotonic way. The application of Müller’s theorem gives the following result:
p_ x ðtÞ ¼
p_ þx ðtÞ ¼
min v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ
v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
max
v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ
The optima can be computed by an optimization algorithm or whenever possible using analytic formulas. Otherwise, suboptimal approaches can be considered to simplify the computation. For instance, in [28], this approach is based on the decomposition of a non-monotonic function into the difference of two increasing functions provided some Lipschitz conditions are satisfied. Another choice could be given by the decoupling into two simpler terms. One has thus:
p_ þx ðtÞ ¼
max
fv 1 ðtÞ cos /ðtÞg
v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ
/ðtÞ 2 / ðtÞ; /þ1 ðtÞ min fv 2 ðtÞ sin /ðtÞg v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ Each term can be more easily computed and the process can be re_ _þ peated for each derivative (e.g. p_ x ; py ; py ; . . .). The results are suboptimal and the intervals may grow in the worst case. In this context, on-going investigations are performed [29,30] and meanwhile, hardware sensors are usually added to improve performance. Finally, a lower bound of the confidence level is computed using þ (19) which is to P xðtÞ 2 x P P ðtÞ; xP ðtÞ reduced þ P xðtk Þ 2 xC ðtk Þ; xC ðt k Þ as the parameters and the input are
known. Consequently, a lower bound of the confidence level for each state at a time t only depends on the joint confidence level at time tk < t of the vehicle position, the vehicle velocity and the slope angle. For instance, in the case of the estimation of the vehi þ þ cle position, 2 p P P xðtÞ 2 x P ðtÞ; pP ðtÞ P ðtÞ;xP ðtÞ P pðtÞ þ þ P Pðpðtk Þ 2 pC ðtk Þ; pC ðt k Þ; v ðtk Þ 2 v C ðt k Þ; v C ðtk Þ; aðtk Þ 2 aC ðt k Þ; aþ C ðt k ÞÞ. Hence, if the intervals at the correction step have a confidence level of 1–104 and are independent, then PðpðtÞ 2 4 3 þ pP ðtÞ; pþ P ðtÞ Þ P P xðt k Þ 2 xC ðt k Þ; xC ðt k ÞÞ ¼ ð1 10 Þ P 1 3 104 . Remark 4.2. If the intervals are dependent, one has the following þ relation: P xðt k Þ 2 x P 1 3 104 . Moreover, as C ðt k Þ; xC ðt k Þ the probability of an error is small, the resulting confidence level is close to the case of independent intervals: (1–104)3 1–3 104.
4.2.2. Correction step In the considered example, a sensor measures the vehicle position (GPS for instance) with some uncertainty and h another (airadar þ for instance) measures the vehicle velocity: y and p ðt k Þ; yp ðt k Þ yv ðtk Þ; yþ ðt Þ . k v At each measurement instant tk, the sensor interval is combined with the predictor interval by union and intersection operations to find the best interval satisfying an integrity objective. As there are two intervals, the set of combinations are: the predicted interval, the sensor interval, the intersection of intervals or the union.
h
i yp ðtk Þ; yþp ðtk Þ and pP ðt k Þ; pþP ðt k Þ ! pC ðtk Þ; pþC ðt k Þ : yv ðt k Þ; yþv ðtk Þ and v P ðt k Þ; v þP ðt k Þ ! v C ðtk Þ; v þC ðtk Þ : If the confidence level of the measurement interval is 1–104 and the one of the predicted interval is 1–3 104, the best interval combination satisfying the objective level 1–104 is the measurement interval only. If the objective level is higher, 1–105 for example, the only possibility is the union of the intervals giving 1–3 108 as a minimum value of the confidence level. Furthermore, the vehicle acceleration along the railway is measured by accelerometers. These measurements allow to determine measurements of the slope angle isolating a in (21):
sin a ¼
1 F j A Bv C v 2 M v_ gM
ð30Þ
Consequently, na intervals measuring acceleration are combined according to the rules explained in Section 3.2 in order to obtain the corrected acceleration interval:
h i h i ya;1 ðtk Þ; yþa;1 ðt k Þ . . . ya;na ðtk Þ; yþa;na ðtk Þ
! aC ðt k Þ; aþC ðtk Þ :
Then, the interval of the slope angle is computed from the accel þ eration interval a ðt Þ; a ðt Þ and the velocity interval k k C C v C ðtk Þ; v þC ðtk Þ obtained in the correction step. The resulting confidence level is the joint probability of the corrected intervals, i.e., the product of the individual integrities as the intervals related to the velocity and the acceleration are independent.
2 1 F j A Bv þC ðt k Þ C v þC ðt k Þ MaþC ðt k Þ gM 2 1 F j A Bv C ðt k Þ C v C ðt k Þ MaC ðt k Þ sin yþa ðt k Þ ¼ gM
sin ya ðt k Þ ¼
ð31Þ ð32Þ
Finally, the interval of the slope angle is obtained by combining the þ predicted interval a computed from P ðt k Þ; aP ðt k Þ and the interval þ the acceleration measurements y a ðt k Þ; ya ðt k Þ .
ya ðt k Þ; yþa ðtk Þ and aP ðtk Þ; aþP ðt k Þ ! aC ðt k Þ; aþC ðtk Þ :
Prediction and correction steps are repeated as new measurements become available.
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
4.2.3. Numerical results The algorithm is illustrated with the simulated trajectory (Fig. 2) with the following characteristics: The initial intervals have the following values: [p(t0), p+(t0)] = [10 m, 10 m], [v(t0), v+(t0)] = [0 m/s, 0.1 m/s] and [a(t0), a+(t0)] = [0.05 rad, 0.05 rad]. Moreover, their confidence level is chosen as 1–104. The vehicle position, velocity and acceleration are measured with a uniform distribution ðUða; bÞÞ of errors and 1 s as sampling period. The size of the uniform window (b–a) is 5 m for the vehicle position, 1 m/s for the vehicle velocity and 0.1 m/ s2 for the acceleration. One sensor is available for position and velocity measurements whereas 3 accelerometers are considered. It is a reasonable configuration as the price of accelerometers (such as MEMS-based accelerometers) is competitive. Each measurement interval has 1–104 as confidence level. In the correction step, the integrity objective is 1–104 for the position and velocity, 1–106 for the acceleration and 1–104 for the slope angle. Fig. 3 shows the results of the estimation errors with the confidence interval observer. In these graphs, the error is represented with respect to the real trajectory. The computed intervals have the following properties: At each measurement instant, the chosen intervals for the vehicle position and velocity come from the sensors and satisfy the
547
integrity objective as the predicted interval has a lower confidence level than the measurement one. Consequently, the corrected position and velocity intervals have 1–104 as confidence level. The interval based on the combination of 3 accelerometers with 1–106 as integrity objective results from the consecutive application of an intersection and a union (c(3, 1, 3, 1), c(3, 1, 3, 2) and c(3, 1, 3, 3) according to the notation in Section 3.2) and a lower bound of the resulting confidence level is 1–2 108. The integrity of the slope angle interval at the correction step is obtained by the product of the corrected velocity interval (1– 104) and the corrected acceleration 1–2 108: (1–104)(1– 2 108) P 1–104–2 108 1–104. In the prediction step, the integrity of the position, the velocity and the slope angle are obtained from the integrity of the state vector at the initial instant, i.e., the integrity of the corrected position (1–104), the corrected velocity (1–104) and the corrected slope angle (1–104). As the intervals are not independent (the corrected slope angle is computed using the corrected velocity), formula (11) is used and one has the following lower bound: 1–3 104. The computational cost of the confidence interval observer is not prohibitive in the example as it uses a limited amount of intervals (4 intervals at most). Indeed, 2 intervals (from the prediction and from the sensor) are compared for the position and the velocity whereas 3 intervals are combined for the acceleration. Then, the intervals bounds are propagated by 2 sets of differential equations.
Fig. 3. Estimation with a confidence interval observer: error intervals. Red: lower and upper bounds. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
548
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
A large number of intervals would be unrealistic as it would imply a large number of sensors and an increase of the investment and operational costs. The following section compares the proposed confidence interval observer with two well-established techniques, Kalman filtering and guaranteed interval estimation. 5. Comparison with a guaranteed interval estimator and an extended Kalman filter In this section, the use of a guaranteed interval estimator and an extended Kalman filter is discussed as compared with the probabilistic observer. The guaranteed and probabilistic interval estimators both propagate intervals in the prediction step, assuming that the true values of state lie in the initial intervals. However, they differ in the correction step as the guaranteed interval observer only performs intersections. This procedure can cause failure of the algorithm in case of ‘‘unlucky’’ outliers implying an empty intersection. On the other hand, Kalman filtering [31] is a popular state estimation method for linear stochastic systems with a Gaussian distribution of the state and measurement noises. EKF, its extension to nonlinear systems, is widely used and generally gives satisfactory results. Based on a linearization along the estimated trajectory and the local application of Kalman filter equations, EKF can be considered as an exponential observer under some conditions [32] but it can fail because of large deviations in the initial conditions, model uncertainties or non-Gaussian noise [33] (uniform, multimodal, unsymmetrical noise distributions for instance).
In the sequel, two situations are considered: The measurement noise is described by a uniform distribution. An outlier is simulated at a given time instant. 5.1. Test 1: uniform measurement noise The first example considers a situation where the measurements are corrupted by noises described by uniform distributions (on given intervals). A typical example is the propagation errors in the ionosphere and troposphere modifying the measurement of a GNSS receiver. An extended Kalman filter is tested using the same set of measurements (given by intervals) as in Section 4.2.3 (Fig. 3). To describe the covariance matrix of the measurement noise (used as a tuning parameter in the EKF), a standard deviation is chosen for each measured variable in such a way that the 99.99%-confidence intervals computed with a Gaussian distribution are in the range of the measurement intervals given by the sensors according to a uniform distribution: 0.65 m for the position, 0.13 m/s for the velocity and 0.013 m/s2 for the acceleration. As in the previous example, three accelerometers are considered. Fig. 4 shows the estimation results with 99.99%-confidence intervals for the estimation errors. As compared to the confidence interval observer (Fig. 3), Kalman filtering generally provides narrower intervals, which is interesting at first sight. However these narrower intervals do not always contain the true values, which is detrimental in security applications. For instance, 8 position
Fig. 4. Estimation with an extended Kalman filter: error intervals. Red: lower and upper bounds. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
549
Fig. 5. Effect of an outlier on a position measurement. Left: EKF, right: confidence interval observer.
Fig. 6. Effect of an outlier on a velocity measurement. Left: EKF, right: confidence interval observer.
Fig. 7. Effect of an outlier on an acceleration measurement. Left: EKF, right: confidence interval observer.
intervals do not enclose the true position at a measurement instant (601 measurements in the simulation). Moreover, 10 velocity intervals and 32 acceleration intervals do not enclose the true values. To alleviate this problem the measurement noise covariance could be increased resulting in larger confidence intervals, still with no guarantee that they contain the true values.
5.2. Test 2: outlier The second example deals with the presence of an outlier. An aberrant measurement can happen according to environmental conditions or due to a sensor fault. For instance, in a GNSS receiver, the multipath effect can provide abnormal measurements because of the local geometrical configuration. Three cases are shown:
550
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
An outlier is simulated for a position measurement at t = 25 s (Fig. 5). An outlier is simulated for a velocity measurement at t = 25 s (Fig. 6). An outlier is simulated for an acceleration measurement at t = 25 s (Fig. 7). The EKF is tested with measurements following a uniform distribution as in the previous test. One can notice that with the considered configuration of the confidence interval observer, recovery after the outlier is faster than with an EKF. It is explained by the fact that the measurement is selected at each measurement instant. Moreover, in the third case, the confidence interval observer is not sensitive to the outlier provided by an accelerometer thanks to the combination of 3 accelerometers. Furthermore, the presence of outliers can cause the failure of a guaranteed interval observer in the correction step. Indeed, the correction step proceeds with intersections and an outlier can lead to a void intersection. 6. Conclusion This paper deals with the design of confidence interval observers applied to vehicle positioning. The main objective is to satisfy security criteria which are critical in the railway or aeronautics area. The method uses a model-based interval predictor and a combination of intersection and union operations at each measurement instant between probabilistic intervals from hardware sensors and the interval predictor. The algorithm is compared to the extended Kalman filter. Even though Kalman filtering is a reference method in estimation, it is not directly applicable for mainly two reasons: (a) Kalman filtering assumes zero-mean white noise on the measurements and the state and does not allow uncertainties on the model parameters and system inputs to be accounted for, (b) Kalman filtering provides confidence intervals on the state estimate trajectories based on the assumption of a Gaussian distribution of the noise, which is not always verified in practice (particularly if GPS signals are used). Moreover, outliers can occur and a prior test should be performed before applying KF equations. Furthermore, the method is also compared with guaranteed interval observers, which are less robust to outliers and do not allow the consideration of security objectives specified in terms of probability. Acknowledgements This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its authors. References [1] M. Grewal, L. Weill, A. Andrews, Global Positioning Systems, Inertial Navigation and Integration, Second ed., John Wiley & Sons, Inc., New Jersey, 2007. [2] B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter – Particle Filters for Tracking Applications, Artech House, Boston, 2004. [3] T.I. Fossen, Marine Control Systems: Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles, Marine Cybernetics AS, Trondheim, 2002.
[4] A. Simsky, F. Wilms, J.-P. Franckart, GNSS-based fail-safe train positioning system for low-density traffic lines based on one-dimensional positioning algorithm, in: Proceedings of NAVITEC 2004 Workshop, Noordwijk, The Netherlands, 2004. [5] L. Jaulin, M. Kieffer, O. Didrit, E. Walter, Applied Interval Analysis, Springer, London, England, 2001. [6] M. Kieffer, E. Walter, Guaranteed nonlinear state estimator for cooperative systems, Numerical Algorithms 37 (1–4) (2004) 187–198. [7] M. Kieffer, E. Walter, Guaranteed nonlinear state estimation for continuoustime dynamical models from discrete-time measurements, in: Proceedings of the 5th IFAC Symposium on Robust Control Design, 2006. [8] L. Jaulin, Robust set-membership state estimation; application to underwater robotics, Automatica 45 (1) (2009) 202–206. [9] N. Ramdani, N. Meslem, Y. Candau, A hybrid bounding method for computing an over-approximation for the reachable set of uncertain nonlinear systems, IEEE Transactions on Automatic Control 54 (10) (2009) 2352–2364. [10] T. Raı¨ssi, G. Videau, A. Zolghadri, Interval observer design for consistency checks of nonlinear continuous-time systems, Automatica 46 (3) (2010) 518– 527. [11] O. Bernard, J.-L. Gouzé, Closed loop observers bundle for uncertain biotechnological models, Journal of Process Control 14 (7) (2004) 765–774. [12] M. Moisan, O. Bernard, J.-L. Gouzé, Near optimal interval observers bundle for uncertain bioreactors, Automatica 45 (1) (2009) 291–295. [13] G. Goffaux, A. Vande Wouwer, O. Bernard, Improving continuous–discrete interval observers with application to microalgae-based bioprocesses, Journal of Process Control 19 (7) (2009) 1182–1190. [14] N. Meslem, N. Ramdani, Y. Candau, Using hybrid automata for set-membership state estimation with uncertain nonlinear continuous-time systems, Journal of Process Control 20 (4) (2010) 481–489. [15] O. Bilenne, Integrity-directed sequential state estimation: assessing high reliability requirements via safe confidence intervals, Information Fusion 8 (1) (2007) 40–55. [16] E. Druet, O. Bilenne, M. Massar, F. Meers, F. Lefèvre, Interval-based state estimation for safe train positioning, in: 7th World Congress on Railway Research, Montreal, Canada, 2006. [17] Y. Zhu, B. Li, Optimal interval estimation fusion based on sensor interval estimates with confidence degrees, Automatica 42 (1) (2006) 101–108. [18] D.J. Berleant, M. Ceberio, G. Xiang, V. Kreinovich, Towards adding probabilities and correlations to interval computations, International Journal of Approximate Reasoning 46 (3) (2007) 499–510. [19] D.J. Berleant, O. Kosheleva, V. Kreinovich, H.T. Nguyen, Unimodality, independence lead to NP-hardness of interval probability problems, Reliable Computing 13 (3) (2007) 261–282. [20] G. Goffaux, A. Vande Wouwer, M. Remy, Vehicle positioning by a confidence interval observer – application to an autonomous underwater vehicle, in: Proceedings of the IEEE Intelligent Vehicles Symposium, Istanbul, Turkey, 2007. [21] O. Bilenne, Integrity-directed sensor fusion by robust risk assessment of faulttolerant interval functions, in: Proceedings of the 9th International Conference on Information Fusion, Florence, Italy, 2006. [22] M. Fréchet, Généralisation du théorème des probabilités totales, Fundamenta Mathematicae 25 (1935) 379–387. [23] M. Müller, Über das fundamentaltheorem in der theorie der gewöhnlichen differentialgleichungen, Mathematische Zeitschrift 26 (1926) 619–645. [24] H. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, American Mathematical Society, 1995. [25] M. Fruchard, O. Bernard, J.-L. Gouzé, Interval observers with guaranteed confidence levels. Application to the activated sludge process, in: Proceedings of the IFAC World Congress, Barcelona, Spain, 2002. [26] R. Kaller, J.-M. Allenbach, Traction Electrique, vol. 1, Presses Polytechniques et Universitaires Romandes, 1995. [27] C. Courtois, J. Coumel, Traction électrique ferroviaire: dynamique ferroviaire et sous-stations, Techniques de l’ingénieur (2009) 1–16. [28] M. Moisan, O. Bernard, Interval observers for non monotone systems. application to bioprocess models, in: 16th IFAC World Congress, Prague, Czech Republic, 2005. [29] F. Mazenc, O. Bernard, Interval observers for linear time-invariant systems with disturbances, Automatica 47 (1) (2011) 140–147. [30] N. Meslem, N. Ramdani, Interval observer design based on nonlinear hybridization and practical stability analysis, International Journal of Adaptive Control and Signal Processing 25 (3) (2011) 228–248. [31] R.E. Kalman, R.S. Bucy, New Results in Linear Filtering and Prediction Theory, Transactions of the ASME, Journal of Basic Engineering 83 (1) (1961) 95–108. [32] K. Reif, R. Unbehauen, The extended Kalman filter as an exponential observer for nonlinear systems, IEEE Transactions in Signal Processing 47 (8) (1999) 2324–2328. [33] E.L. Haseltine, J.B. Rawlings, Critical evaluation of extended Kalman filtering and moving-horizon estimation, Industrial & Engineering Chemistry Research 44 (8) (2005) 2451–2460.
Contents lists available at SciVerse ScienceDirect
Information Fusion journal homepage: www.elsevier.com/locate/inffus
Continuous–discrete confidence interval observer – Application to vehicle positioning G. Goffaux, M. Remy, A. Vande Wouwer ⇑ Service d’Automatique, Université de Mons (UMONS), Mons, Belgium
a r t i c l e
i n f o
Article history: Received 27 November 2012 Received in revised form 12 January 2013 Accepted 13 February 2013 Available online 28 February 2013 Keywords: Intervals State observer Positioning systems Railways Vehicles
a b s t r a c t In vehicle positioning applications, the confidence level in the position and velocity estimates can be even more significant than accuracy. In this study, a probabilistic interval method is proposed, which combines, through union and intersection operations, the information from a possibly uncertain predictor (the vehicle model) and measurement sensors. The proposed method is compared to Kalman filtering and to guaranteed interval estimation in the context of railway vehicles where security is the key objective. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Position and velocity estimation is generally performed using data fusion algorithms which combine information from a mathematical model and sensor measurements. Popular estimation methods include Kalman filtering based on a kinematic model of the vehicle and measurements from GPS (Global Positioning System) and INS (Inertial Navigation System) [1], particle filtering for positioning, navigation and tracking [2] or nonlinear filtering in marine applications [3]. Accuracy and security are two conflicting objectives. However, security requirements are particularly important in vehicle positioning as collisions could cause human and economical loss. Hence, the estimated variables have to be defined by confidence intervals in applications such as safe railway positioning [4]. Recently, state estimation methods based on interval algebra [5] have been proposed to compute guaranteed intervals upper bounding the solution sets in the case of uncertain discrete-time models. Subsequent important contributions in this area of research include the development of methods for nonlinear continuous-time systems, see for instance [6–10]. These developments have been of particular interest in the study of biological systems where the models are nonlinear and uncertain [11–14]. However, these latter results do not give a realistic measure of the confidence level as they consider guaranteed intervals. Other studies have preferred a statistical approach for tackling the prob⇑ Corresponding author. Tel.: +32 (0)65374141; fax: +32 (0)65374136. E-mail address: [email protected] (A. Vande Wouwer). 1566-2535/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.inffus.2013.02.006
lem of continuous–discrete estimation with faulty sensors [15], either by considering the manipulation of confidence intervals in a static configuration [16,17], or uncertainties in the probabilities (for example, a variance interval) [18,19]. In this study, a continuous–discrete confidence interval estimator is proposed, which proceeds in two steps. The first step uses a continuous-time predictor to compute the evolution of lower and upper bounds for each state variable between two discrete measurement times. The second step occurs at a measurement instant, where intervals from the predictor and from hardware sensors are combined using union and intersection operations so as to achieve specified confidence levels. In [20], a preliminary version of the algorithm has been designed in the context of a marine vehicle described by a nonlinear dynamic model. The effect of state transformations on the nonlinear predictor has been studied and a security measure has been defined which corresponds to a given level of error tolerance within a given time span. In this study, the method is further developed in the context of a railway application where the security demands are more severe. The prediction step uses a dynamic model based on the equilibrium of forces. In addition, several sensors are considered, providing position, velocity or acceleration information. Moreover, a comparison with two well-established techniques, Kalman filtering and guaranteed interval estimation, is achieved so as to highlight the advantages of the proposed probabilistic interval approach. In the sequel, the paper is organized in four main sections. In Section 2, some preliminary model assumptions are introduced
542
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
along with the concept of confidence intervals. Section 3 is devoted to the general structure of the continuous–discrete confidence interval observer. In Section 4, the observer is applied to a test example, the positioning of a railway vehicle. After a description of the vehicle model, a dedicated interval observer is proposed. In Section 5, the probabilistic interval observer is compared with an extended Kalman filter (EKF) and a guaranteed interval observer.
Consider the vehicle-sensors system described by the following equations for t P t0:
ð RÞ :
To compute the resulting confidence levels, the intervals are associated with Boolean random variables denoted B which are defined as follows [21]:
Bx ¼ true if x 2 ½x ; xþ x; x ; xþ 2 R;
x < xþ
Bx ¼ false otherwise
ð7Þ
Consequently, the probability of enclosing a variable x by an interval is related to the corresponding Boolean random variable:
2. Preliminaries
(
2.1. Intersection and union operations
_ xðtÞ ¼ fðxðtÞ; uðtÞ; hÞ xðt0 Þ ¼ x0
ð1Þ
yk ¼ yðt k Þ ¼ Cðt k Þxðt k Þ þ ðt k Þ
The first part of (R) models the time evolution of positioning variables and consists of (possibly nonlinear) differential equations. xðtÞ 2 X Rnx is the vector of state variables (including position, velocity and acceleration) with X a compact set of Rnx . x0 is the initial state vector. In this study, prior confidence intervals are defined by the probability (denoted P) that the corresponding bounds include the unknown values, at the initial time with confidence level b0.
P x0 2 x0 ; xþ0 P b0 () h i h i P x0;1 2 x0;1 ; xþ0;1 ; . . . ; x0;nx 2 x0;nx ; xþ0;nx P b0
ð8Þ
From this information, union and intersection operations can be defined in terms of Boolean random variables [21,22]. We consider two measuring a variable x with the intervals independent þsensors x1 ; xþ 1 and x2 ; x2 . Bx1 and Bx2 are the corresponding Boolean random variables and bx1 and bx2 are the lower bounds on the integrity levels.
PðBx1 ¼ trueÞ ¼ Pðx 2 x1 ; xþ1 Þ P bx1 þ PðBx2 ¼ trueÞ ¼ Pðx 2 x2 ; x2 Þ P bx2 PðBx1 ¼ falseÞ ¼ Pðx R x1 ; xþ1 Þ 6 1 bx1 PðBx2 ¼ falseÞ ¼ Pðx R x2 ; xþ2 Þ 6 1 bx2
ð2Þ Intersection:
nu
uðtÞ 2 U R is the vector of system inputs including, for instance, forces controlling the vehicle motion. U is the set of admissible input signals, a subset of the space of measurable bounded piecewise constant functions. Consequently, utk !t is a piecewise constant function between tk and t which can be determined by the knowledge of a finite number of input values u(tk), . . . , u(t). In this study, they are known with some uncertainties such that confidence intervals are defined by
h i P utk !t 2 utk !t ; uþtk !t ()
PðBx ¼ trueÞ ¼ Pðx 2 ½x ; xþ Þ P bx
P x 2 x1 ; xþ1 \ x2 ;xþ2 ¼ PðBx1 ¼ true and Bx2 ¼ trueÞ P bx1 bx2
ð9Þ
Union: P x 2 x1 ;xþ1 [ x2 ;xþ2 ¼ PðBx1 ¼ true or Bx2 ¼ trueÞ ¼ 1 PðBx1 ¼ false and Bx2 ¼ falseÞ P 1 ð1 bx1 Þð1 bx2 Þ ¼ bx1 þ bx2 bx1 bx2
ð10Þ
Pðuðt k Þ 2 ½u ðt k Þ; u ðt k Þ; . . . ; uðtÞ 2 ½u ðtÞ; u ðtÞÞ P buk ðtÞ
Remark 2.1. If x is a nx-dimension vector, its confidence level is defined by the joint probability that the intervals related to its components include the unknown values. x ¼ ½x1 ; . . . ; xnx T
Moreover, h 2 H 2 Rnh is a model parameter vector and again confidence intervals are given (e.g. by a parameter identification procedure).
Pðx 2 ½x ; xþ Þ ¼ P x1 2 x1 ; xþ1 ; . . . ; xnx 2 xnx ; xþnx
Pðh 2 ½h ; hþ Þ P bh
Remark 2.2. If there is a dependence between two intervals, (9) and (10) become
þ
ð3Þ
þ
ð4Þ
The second part of (R) describes the measurement process. In vehicle positioning, sensors include accelerometers, gyroscopes, radars, GNSS (Global Navigation Satellite System) receivers providing partial information on vehicle position, velocity and acceleration. yðtk Þ 2 Rny is the vector of sampled measurements. Usually, measurements are modelled using additive random errors (tk) with Gaussian distribution, but other situations can be considered as well. Confidence intervals are given by
P yk 2 yk ; yþk P byk
P x 2 x1 ; xþ1 \ x2 ;xþ2 ¼ PðBx1 ¼ true and Bx2 ¼ trueÞ ¼ PðBx1 ¼ trueÞ þ PðBx2 ¼ trueÞ PðBx1 ¼ true or Bx2 ¼ trueÞ PPðBx1 ¼ trueÞ PðBx2 ¼ falseÞ Pbx1 þ bx2 1
ð5Þ
ð11Þ
Union:
Typically, each sensor has its own sampling period and measurements are asynchronous. Thus, tk is a measurement instant where C(tk) is the observation matrix corresponding to the active sensors. From these definitions, a confidence interval observer is designed in order to derive bounds, x(t) and x+(t), with a specified confidence level bobj (t P t0):
PðxðtÞ 2 ½x ðtÞ; xþ ðtÞÞ P bobj
Intersection:
ð6Þ
P x 2 x1 ; xþ1 [ x2 ;xþ2 ¼ PðBx1 ¼ true or Bx2 ¼ trueÞ P maxðbx1 ; bx2 Þ
ð12Þ
Remark 2.3. If the intervals are disconnected, the union operation gives the convex hull. For instance, [2,3] and [4,5] gives [2,5]. In the sequel, these operations are used to develop a confidence interval observer.
543
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
3. Continuous–discrete confidence interval observers
x_ i ðt k Þ 6
The algorithm introduced in [20] proceeds in two steps. Between two measurements times, the state intervals are propagated (or predicted) using the model equations. Then, the predicted intervals and the measurement intervals are combined to produce corrected intervals.
x_ þi ðt k Þ P
(
Ti ðtk Þ :
x_ P ðtÞ ¼ f xP ðtÞ; xþP ðtÞ; u ðtÞ; uþ ðtÞ; h ; hþ þ x_ þP ðtÞ ¼ f xP ðtÞ; xþP ðtÞ; u ðtÞ; uþ ðtÞ; h ; hþ
ð13Þ
xþP ðt k Þ ¼ xþC ðtk Þ
where f and f+ are derived from the model equations, and depend on bounds on the state variables, inputs and parameters. Before computing the resulting confidence level, the predictor equations have to satisfy the following property. Property 3.1 (Enclosure of the prediction equations). If + þ þ xðtk Þ 2 x P ðt k Þ; xP ðt k Þ ; uðtÞ 2 ½u ðtÞ; u ðtÞ and h 2 [h ,h ], then the solution of the prediction equations f and f+ (given by (13)) encloses þ the state variables if it satisfies xðtÞ 2 x P ðtÞ; xP ðtÞ 8t P t k . This property is related to the definition of an error vector eðtÞ 2 R2nx which has to remain nonnegative.
eðtÞ ¼
þ
xP ðtÞ xðtÞ eþ ðtÞ ¼ P 0 8t P t k xðtÞ xP ðtÞ e ðtÞ
ð14Þ
Proof. It is sufficient to check the following condition, assuming e(tk) > 0
If 9t P tk jei ðt Þ ¼ 0;
ej ðt Þ P 0 8j – i 2 f1; . . . ; 2nx g; then e_ i ðt Þ P 0
ð15Þ
Each time a component ei of the error vector vanishes, its derivative must be nonnegative in order to prevent a negative error [20]. f and f+ are derived from the evolution equation f using Müller’s theorem [23,7]. h T
Theorem 3.1. Assume that the vector of functions f ¼ ½f1 ; . . . ; fnx in (1) is continuous on a domain
8 > < t0 6 t 6 tf T : x ðtÞ 6 xðtÞ 6 xþ ðtÞ; > : h 2 H; uðtÞ 2 U
x ðtÞ; xðtÞ; xþ ðtÞ 2 X
T þ where x ðtÞ ¼ and x ðtÞ ¼ xþ are 1 ðtÞ; . . . ; xnx ðtÞ continuous on [t0, tf] and such that
T x 1 ðtÞ; . . . ; xnx ðtÞ
Tþi ðtk Þ :
þ
þ þ 1. x ðt0 Þ ¼ x 0 and x ðt 0 Þ ¼ x0 2. the derivatives of x(tk) and x+(tk) at tk (t0 6 tk 6 tf) satisfy for i = 1, . . . , nx,
max fi ðxðt k Þ; uðt k Þ; hÞ
ð17Þ
xðtk Þ2Tþ ðt k Þ i
8 t ¼ tk > > > > > > > < xi ðtÞ ¼ xi ðtÞ > > xj ðtÞ 6 xj ðtÞ 6 xþj ðtÞ; > > > > > : h 2 H; uðtÞ 2 U
j–i
8 t ¼ tk > > > > > < xi ðtÞ ¼ xþ ðtÞ i > xj ðtÞ 6 xj ðtÞ 6 xþj ðtÞ; > > > > : h 2 H; uðtÞ 2 U
j–i
þ Then, for any x0 2 x 0 ; x0 , a solution to (1) exists, which remains in
E:
xP ðt k Þ ¼ xC ðtk Þ
ð16Þ
þ where T i ðt k Þ and Ti ðt k Þ are subsets of T respectively defined by:
3.1. Prediction step The prediction step aims at propagating the interval bounds for each state variable from a start time tk using the information from the model equations only. For the construction of the confidence interval predictor, it is first assumed that guaranteed intervals exist for the initial state vector, the inputs and the parameters. Advanced methods for the design of guaranteed interval predictors can be employed [8–10]. Subsequently, this assumption is relaxed considering confidence levels for the variables and parameters. The typical form of an interval predictor consists of two coupled dynamic systems (t P tk) which are initialized with intervals ob þ tained in the previous correction step: x C ðt k Þ; xC ðt k Þ .
min fi ðxðt k Þ; uðt k Þ; hÞ
xðt k Þ2T ðtk Þ i
t0 6 t 6 tf
ð18Þ
x ðtÞ 6 xðtÞ 6 xþ ðtÞ
Proof. To prove relation (18), condition (15) is checked. By definiþ tion of T i ðt k Þ and Ti ðt k Þ, condition (15) is satisfied. For instance, at ðt Þ are the true state vector and the true ðt Þ and u time t = t⁄, if x input respectively along with x+(t⁄) the estimated upper bounds, one has:
if eþi ðt Þ ¼ xþi ðt Þ xi ðt Þ ¼ 0 and xj ðt Þ 6 xj ðt Þ 6 xþj ðt Þ ðj – iÞ
e_ þi ðt Þ P
ðt Þ; hÞ P 0 ðt Þ; u max fi ðxðt Þ; uðt Þ; hÞ fi ðx
xðt Þ2Tþ ðt Þ i
with
8 t ¼ t > > > < xi ðtÞ ¼ xþi ðtÞ ¼ xi ðtÞ Tþi ðt Þ : þ > > xj ðtÞ 6 xj ðtÞ 6 xj ðtÞ; > : h 2 H; uðtÞ 2 U
j–i
Theorem 3.1 can be extended to uncertain parameters and inputs. Moreover, a specific version of Theorem 3.1 can be obtained when f satisfies monotonicity properties [24]. Finally, advanced guaranteed interval predictors could be easily applied [13,14,10], if it is required by the complexity of the system model. Remark 3.1. Interval predictors can be designed from additional models such as kinematic and dynamic models of a vehicle. The resulting intervals can be combined if each of the interval predictor satisfies Property 3.1. Extending the previous theorem to confidence intervals, the following corollary gives a lower bound of the confidence level related to the intervals computed by an interval predictor [20]. Corollary 3.1. If confidence intervals are known at time tk (related to h i þ þ and the the state vector x C ðt k Þ; xC ðt k Þ , the inputs ut k !t ; ut k !t parameters [h,h+]) and if the prediction equations (f and f+) satisfy Property 3.1 of enclosure, the intervals on the predicted state þ variables, computed by the interval predictor x P ðtÞ; xP ðtÞ , satisfy the following lower bound of the integrity level [25] (tk < t 6 tk+1)
544
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
P xðtÞ 2 xP ðtÞ; xþP ðtÞ P P xðt k Þ 2 xC ðtk Þ; xþC ðtk Þ ; utk !t h i 2 utk !t ; uþtk !t ; h 2 ½h ; hþ P bxk buk ðtÞbh
ð19Þ
with bxk ; bh and buk the integrity related, respectively, to the state vector at tk, the parameters and the inputs (independence is assumed). Moreover, the inequality is also satisfied for each component of the state vector. Indeed,
h i P xi ðtÞ 2 xi;P ðtÞ; xþi;P ðtÞ P 8i 2 1; . . . ; nx h i P xðtk Þ 2 xC ðt k Þ; xþC ðt k Þ ; utk !t 2 utk !t ; uþtk !t ; h 2 ½h ; hþ P
3 intervals (i = 3): there is only one selection (j = 1) of intervals (3 elements out of a set of 3 elements) but there are four possible sequences of operations (from k = 1 to k = 4): {\, \}, {[, [}, {[, \}, {\, [} and also three possible orderings of the intervals when applying this sequence of operations: {1, 2, 3}, {1, 3, 2}, {2, 3, 1}. Consequently, a fourth dimension is added to the notation: c(i, j, k, l) which represents the possible orderings of the selected intervals. Hence, one has the following operations (the operations are performed from left to right):
cð3; 1; 1Þ ¼ ½x1 \ ½x2 \ ½x3 ¼ ½5:3; 13:6 b P 1 3 103
ð20Þ
cð3; 1; 2Þ ¼ ½x1 [ ½x2 [ ½x3 ¼ ½3:9; 15:4 b ¼ 1 109
bxk buk ðtÞbh
cð3; 1; 3; 1Þ ¼ ½x1 \ ½x2 [ ½x3 ¼ ½3:9; 15:4 b P 1 2 106
Proof. It is a consequence of Property 3.1. Indeed, f and f+ are constructed in such a way that the computed bounds x P ðtÞ and xþ P ðtÞ enclose the unknown state trajectory provided that the inputs, the parameter and the initial state vector are within their respective interval. h
cð3; 1; 3; 2Þ ¼ ½x1 \ ½x3 [ ½x2 ¼ ½4:4; 14:1 b P 1 2 106 cð3; 1; 3; 3Þ ¼ ½x2 \ ½x3 [ ½x1 ¼ ½4:4; 14:1 b P 1 2 106 cð3;1; 4; 1Þ ¼ ½x1 [ ½x2 \ ½x3 ¼ ½4:4; 14:1 b P 1 103 106 cð3;1; 4; 2Þ ¼ ½x1 [ ½x3 \ ½x2 ¼ ½4:4; 13:6 b P 1 103 106
3.2. Correction step
cð3;1; 4; 3Þ ¼ ½x2 [ ½x3 \ ½x1 ¼ ½5:3; 14:1 b P 1 103 106
The correction step proceeds by union and intersection operations between intervals from hardware sensors and intervals given by the predictor at the considered measurement time. A lower bound of the resulting confidence levels is computed by (9) and (10) in case of independence. For each state variable, intervals are built considering each possible combination (union and intersection) among the available confidence intervals representing the state variable. Then, the best resulting interval is derived selecting the smallest nonzero interval satisfying the integrity objective. To illustrate the methodology, an example is considered with 3 intervals, [x1] = [5.3, 14.1], [x2] = [4.4, 13.6] and [x3] = [3.9, 15.4]. Each of them has a confidence level of 1–103. At first, there are 3 possible choices involving 1, 2 or 3 intervals. 1 interval: the set of combined intervals, denoted c, only contains the original intervals with their corresponding confidence level. The notation is as follows: c(i, j) has two dimensions, i corresponding to the number of considered intervals (here i = 1) and j to the selected interval (here, j = 1, 2, or 3). 3
cð1; 1Þ ¼ ½x1 ¼ ½5:3; 14:1 b ¼ 1 10 cð1; 2Þ ¼ ½x2 ¼ ½4:4; 13:6 b ¼ 1 103 cð1; 3Þ ¼ ½x3 ¼ ½3:9; 15:4 b ¼ 1 103 2 intervals: with two intervals (i = 2), there are two possible operations, either an intersection or a union. The number of paired intervals is given by the number of combinations of 2 elements in a set of three elements. Now, c(i, j, k) has three dimensions. j indicates the paired interval and k the selected operations: intersection (k = 1) or union (k = 2). Consequently, one has:
cð2; 1; 1Þ ¼ ½x1 \ ½x2 ¼ ½5:3; 13:6 cð2; 1; 2Þ ¼ ½x1 [ ½x2 ¼ ½4:4; 14:1 cð2; 2; 1Þ ¼ ½x1 \ ½x3 ¼ ½5:3; 14:1 cð2; 2; 2Þ ¼ ½x1 [ ½x3 ¼ ½3:9; 15:4 cð2; 3; 1Þ ¼ ½x2 \ ½x3 ¼ ½4:4; 13:6 cð2; 3; 2Þ ¼ ½x2 [ ½x3 ¼ ½3:9; 15:4
b P 1 2 103 b ¼ 1 106 b P 1 2 103 b ¼ 1 106 b P 1 2 103 b ¼ 1 106
The resulting confidence level of an intersection is obtained from (9): b = (1–103)2 P 1–2 103. The resulting confidence level of a union is obtained from (10): b = 2(1–103)– (1–103)2 = 1–106. Also, with more than 2 intervals, the order of the operations has an importance, as it will be seen in the sequel.
The resulting confidence level of two consecutive intersections is obtained by: b = (1–103)3 P 1–3 103 whereas for two consecutive unions: b = (1–103) + (1–106)–(1– 103)(1–106) = 1–109. Also, one can note that consecutive unions or intersections are independent of the order of the intervals. The resulting confidence level of the sequence {\, [} is obtained by: b = (1–103)2 + (1–103)–(1–103)2(1–103) = 1– 2 106 + 109 P 1–2 106 whereas for the sequence {[, \}: b = (1–106)(1–103) P 1–103–106. To conclude the example, if the integrity objective is 1–106, the smallest interval satisfying the objective is [x1] [ [x2] = [4.4, 14.1]. A similar procedure can be performed for a larger number of intervals. Usually, the number of sensors is limited and combinations of 3 or 4 intervals are sufficient. 4. Vehicle positioning: a railway vehicle example In this section, the model of the illustrative example is introduced. Then, the confidence interval observer is developed and simulation tests illustrate the method potentials. 4.1. Railway vehicle: a dynamic model Considering the equilibrium of forces acting on the railway vehicle (Fig. 1) along the railroad and applying the fundamental principles of dynamics, one has [26,27]
M
dv ¼ F j R Mg sin a dt
ð21Þ
In this equation, v is the vehicle velocity and Fj is the force transmitted by the motor to the vehicle wheels. In this simplified model, no loss of adhesion is considered and consequently, the torque delivered by the motor is totally converted into force. R contains all the resisting forces opposing the vehicle movement. In general, these forces depend on the vehicle velocity and can be represented by a formula including an independent term, a term proportional to the velocity and a quadratic term:
R ¼ A þ Bv þ C v 2
ð22Þ
545
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
In (25), r(p) is the derivative of the slope angle with respect to the vehicle position. In the sequel, the state vector x is used to represent the vehicle position p, the vehicle velocity v and the slope angle a. Fig. 2 simulates the evolution equations during 600 s with the parameters given by Table 1. The slope changes with the vehicle position and the force Fj is computed so as to achieve a constant velocity. 4.2. Confidence interval observer Fig. 1. Forces applied to the vehicle.
Table 1 Model parameters (Eurostar train with 20 vehicles) [27]. Description
Parameter
Value
Unit
Vehicle mass Earth’s gravity
M g
816 10
t m/s2
Coefficients of resisting forces R
A B C
4.82 0.23508 0.013608
kN t/s t/m
The term A + Bv is related to mechanical and rolling frictions. A depends on the load applied to the wheels and B gives information on the quality of the railway track and on the vehicle stability. Finally, Cv2 represents the aerodynamic forces in the direction of the vehicle movement. C depends on the drag coefficient. The latter varies with the design of the vehicle profile. Moreover, the gravity force projected into the direction of motion is given by Mg sin a where a is the track slope varying according to the train position. Numerical values used in this simulation example are given in Table 1. The evolution equations can be summarized as follows:
p_ ¼ v 1 v_ ¼ ðF j A Bv C v 2 Þ g sin a M a_ ¼ rðpÞv
In the following, a confidence interval observer is designed based on the evolution Eqs. (23–25) and sensors measuring position, velocity and acceleration. At first, the interval predictor is proposed based on Müller’s theorem. Then, the correction step and the combination of confidence intervals are introduced. 4.2.1. Interval predictor Lower and upper bounds are computed by integrating a couple of differential equations related to the vehicle position, velocity and acceleration. The equations related to the vehicle position are obtained maximizing the vehicle velocity as a consequence of Müller’s theorem:
p_ P ðtÞ ¼ p_ þP ðtÞ ¼
ð25Þ
v ðtÞ ¼ v P ðtÞ
max
v ðtÞ ¼ v þP ðtÞ
v ðtÞ2½v P ðtÞ;v þP ðtÞ
ð26Þ
Moreover, assuming no uncertainties on the parameter values and the input Fj, applying Müller’s theorem gives the following derivative for the lower bound of the vehicle velocity.
v_ P ðtÞ ¼
ð23Þ ð24Þ
min
v ðtÞ2½v P ðtÞ;v þP ðtÞ
v_ P ðtÞ ¼
Fj A B C 2 min v ðtÞ M v ðtÞg sin aðtÞ pðtÞ 2 pP ðtÞ;pþP ðtÞ M M M v ðtÞ ¼ v ðtÞ P þ aðtÞ 2 aP ðtÞ; aP ðtÞ Fj A B C v ðtÞ v P 2 ðtÞ g sin aþP ðtÞ M M M M P
Similarly, one has
Fig. 2. Simulation of the dynamic model.
ð27Þ
546
v_ þP ðtÞ ¼
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
Fj A B C 2 v þ ðt Þ v þP ðtÞ g sin aP ðtÞ M M M M P
ð28Þ
In this formulation, a(t) lies between p2 and p2 such that sina is monotonic. In the same way, the derivatives of the bounds of the slope angle are computed:
a_ P ðtÞ ¼
min
v ðtÞ 2 v P ðtÞ; v þP ðtÞ
rðtÞ v ðtÞ
rðtÞ 2 ½rP ðtÞ; rþP ðtÞ a_ þP ðtÞ ¼
max rðtÞv ðtÞ v ðtÞ 2 v P ðtÞ; v þP ðtÞ rðtÞ 2 rP ðtÞ; rþP ðtÞ
ð29Þ
In (29), the optimum values are obtained from the product of two intervals. The first one comes from the vehicle velocity v(t) and the second one from the derivative r(t) of the slope angle with respect to the position p(t). The latter is obtained from the interval of the vehicle position considering the worst slope angle derivatives. For instance, r P ðtÞ ¼ minpðtÞ2½p ðtÞ;pþ ðtÞ rðpðtÞÞ. P
P
Remark 4.1. The model equations can be more complicated depending on the application and the desired precision. However, Müller’s theorem can still be applied in such cases. As an example, consider a part of the 2D evolution equations of a vessel [3]
p_ x ðtÞ ¼ v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
v1(t) and v2(t) are the velocity components given in a frame related to the vessel movement. /(t) is the vehicle yaw angle whereas px is the vessel position in the x-axis of the reference frame (see [3] for more details). In these equations, /(t) appears in a non-monotonic way. The application of Müller’s theorem gives the following result:
p_ x ðtÞ ¼
p_ þx ðtÞ ¼
min v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ
v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
max
v 1 ðtÞ cos /ðtÞ v 2 ðtÞ sin /ðtÞ
v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ
The optima can be computed by an optimization algorithm or whenever possible using analytic formulas. Otherwise, suboptimal approaches can be considered to simplify the computation. For instance, in [28], this approach is based on the decomposition of a non-monotonic function into the difference of two increasing functions provided some Lipschitz conditions are satisfied. Another choice could be given by the decoupling into two simpler terms. One has thus:
p_ þx ðtÞ ¼
max
fv 1 ðtÞ cos /ðtÞg
v 1 ðtÞ 2 v 1 ðtÞ; v þ1 ðtÞ
/ðtÞ 2 / ðtÞ; /þ1 ðtÞ min fv 2 ðtÞ sin /ðtÞg v 2 ðtÞ 2 v 2 ðtÞ; v þ2 ðtÞ /ðtÞ 2 / ðtÞ; /þ1 ðtÞ Each term can be more easily computed and the process can be re_ _þ peated for each derivative (e.g. p_ x ; py ; py ; . . .). The results are suboptimal and the intervals may grow in the worst case. In this context, on-going investigations are performed [29,30] and meanwhile, hardware sensors are usually added to improve performance. Finally, a lower bound of the confidence level is computed using þ (19) which is to P xðtÞ 2 x P P ðtÞ; xP ðtÞ reduced þ P xðtk Þ 2 xC ðtk Þ; xC ðt k Þ as the parameters and the input are
known. Consequently, a lower bound of the confidence level for each state at a time t only depends on the joint confidence level at time tk < t of the vehicle position, the vehicle velocity and the slope angle. For instance, in the case of the estimation of the vehi þ þ cle position, 2 p P P xðtÞ 2 x P ðtÞ; pP ðtÞ P ðtÞ;xP ðtÞ P pðtÞ þ þ P Pðpðtk Þ 2 pC ðtk Þ; pC ðt k Þ; v ðtk Þ 2 v C ðt k Þ; v C ðtk Þ; aðtk Þ 2 aC ðt k Þ; aþ C ðt k ÞÞ. Hence, if the intervals at the correction step have a confidence level of 1–104 and are independent, then PðpðtÞ 2 4 3 þ pP ðtÞ; pþ P ðtÞ Þ P P xðt k Þ 2 xC ðt k Þ; xC ðt k ÞÞ ¼ ð1 10 Þ P 1 3 104 . Remark 4.2. If the intervals are dependent, one has the following þ relation: P xðt k Þ 2 x P 1 3 104 . Moreover, as C ðt k Þ; xC ðt k Þ the probability of an error is small, the resulting confidence level is close to the case of independent intervals: (1–104)3 1–3 104.
4.2.2. Correction step In the considered example, a sensor measures the vehicle position (GPS for instance) with some uncertainty and h another (airadar þ for instance) measures the vehicle velocity: y and p ðt k Þ; yp ðt k Þ yv ðtk Þ; yþ ðt Þ . k v At each measurement instant tk, the sensor interval is combined with the predictor interval by union and intersection operations to find the best interval satisfying an integrity objective. As there are two intervals, the set of combinations are: the predicted interval, the sensor interval, the intersection of intervals or the union.
h
i yp ðtk Þ; yþp ðtk Þ and pP ðt k Þ; pþP ðt k Þ ! pC ðtk Þ; pþC ðt k Þ : yv ðt k Þ; yþv ðtk Þ and v P ðt k Þ; v þP ðt k Þ ! v C ðtk Þ; v þC ðtk Þ : If the confidence level of the measurement interval is 1–104 and the one of the predicted interval is 1–3 104, the best interval combination satisfying the objective level 1–104 is the measurement interval only. If the objective level is higher, 1–105 for example, the only possibility is the union of the intervals giving 1–3 108 as a minimum value of the confidence level. Furthermore, the vehicle acceleration along the railway is measured by accelerometers. These measurements allow to determine measurements of the slope angle isolating a in (21):
sin a ¼
1 F j A Bv C v 2 M v_ gM
ð30Þ
Consequently, na intervals measuring acceleration are combined according to the rules explained in Section 3.2 in order to obtain the corrected acceleration interval:
h i h i ya;1 ðtk Þ; yþa;1 ðt k Þ . . . ya;na ðtk Þ; yþa;na ðtk Þ
! aC ðt k Þ; aþC ðtk Þ :
Then, the interval of the slope angle is computed from the accel þ eration interval a ðt Þ; a ðt Þ and the velocity interval k k C C v C ðtk Þ; v þC ðtk Þ obtained in the correction step. The resulting confidence level is the joint probability of the corrected intervals, i.e., the product of the individual integrities as the intervals related to the velocity and the acceleration are independent.
2 1 F j A Bv þC ðt k Þ C v þC ðt k Þ MaþC ðt k Þ gM 2 1 F j A Bv C ðt k Þ C v C ðt k Þ MaC ðt k Þ sin yþa ðt k Þ ¼ gM
sin ya ðt k Þ ¼
ð31Þ ð32Þ
Finally, the interval of the slope angle is obtained by combining the þ predicted interval a computed from P ðt k Þ; aP ðt k Þ and the interval þ the acceleration measurements y a ðt k Þ; ya ðt k Þ .
ya ðt k Þ; yþa ðtk Þ and aP ðtk Þ; aþP ðt k Þ ! aC ðt k Þ; aþC ðtk Þ :
Prediction and correction steps are repeated as new measurements become available.
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
4.2.3. Numerical results The algorithm is illustrated with the simulated trajectory (Fig. 2) with the following characteristics: The initial intervals have the following values: [p(t0), p+(t0)] = [10 m, 10 m], [v(t0), v+(t0)] = [0 m/s, 0.1 m/s] and [a(t0), a+(t0)] = [0.05 rad, 0.05 rad]. Moreover, their confidence level is chosen as 1–104. The vehicle position, velocity and acceleration are measured with a uniform distribution ðUða; bÞÞ of errors and 1 s as sampling period. The size of the uniform window (b–a) is 5 m for the vehicle position, 1 m/s for the vehicle velocity and 0.1 m/ s2 for the acceleration. One sensor is available for position and velocity measurements whereas 3 accelerometers are considered. It is a reasonable configuration as the price of accelerometers (such as MEMS-based accelerometers) is competitive. Each measurement interval has 1–104 as confidence level. In the correction step, the integrity objective is 1–104 for the position and velocity, 1–106 for the acceleration and 1–104 for the slope angle. Fig. 3 shows the results of the estimation errors with the confidence interval observer. In these graphs, the error is represented with respect to the real trajectory. The computed intervals have the following properties: At each measurement instant, the chosen intervals for the vehicle position and velocity come from the sensors and satisfy the
547
integrity objective as the predicted interval has a lower confidence level than the measurement one. Consequently, the corrected position and velocity intervals have 1–104 as confidence level. The interval based on the combination of 3 accelerometers with 1–106 as integrity objective results from the consecutive application of an intersection and a union (c(3, 1, 3, 1), c(3, 1, 3, 2) and c(3, 1, 3, 3) according to the notation in Section 3.2) and a lower bound of the resulting confidence level is 1–2 108. The integrity of the slope angle interval at the correction step is obtained by the product of the corrected velocity interval (1– 104) and the corrected acceleration 1–2 108: (1–104)(1– 2 108) P 1–104–2 108 1–104. In the prediction step, the integrity of the position, the velocity and the slope angle are obtained from the integrity of the state vector at the initial instant, i.e., the integrity of the corrected position (1–104), the corrected velocity (1–104) and the corrected slope angle (1–104). As the intervals are not independent (the corrected slope angle is computed using the corrected velocity), formula (11) is used and one has the following lower bound: 1–3 104. The computational cost of the confidence interval observer is not prohibitive in the example as it uses a limited amount of intervals (4 intervals at most). Indeed, 2 intervals (from the prediction and from the sensor) are compared for the position and the velocity whereas 3 intervals are combined for the acceleration. Then, the intervals bounds are propagated by 2 sets of differential equations.
Fig. 3. Estimation with a confidence interval observer: error intervals. Red: lower and upper bounds. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
548
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
A large number of intervals would be unrealistic as it would imply a large number of sensors and an increase of the investment and operational costs. The following section compares the proposed confidence interval observer with two well-established techniques, Kalman filtering and guaranteed interval estimation. 5. Comparison with a guaranteed interval estimator and an extended Kalman filter In this section, the use of a guaranteed interval estimator and an extended Kalman filter is discussed as compared with the probabilistic observer. The guaranteed and probabilistic interval estimators both propagate intervals in the prediction step, assuming that the true values of state lie in the initial intervals. However, they differ in the correction step as the guaranteed interval observer only performs intersections. This procedure can cause failure of the algorithm in case of ‘‘unlucky’’ outliers implying an empty intersection. On the other hand, Kalman filtering [31] is a popular state estimation method for linear stochastic systems with a Gaussian distribution of the state and measurement noises. EKF, its extension to nonlinear systems, is widely used and generally gives satisfactory results. Based on a linearization along the estimated trajectory and the local application of Kalman filter equations, EKF can be considered as an exponential observer under some conditions [32] but it can fail because of large deviations in the initial conditions, model uncertainties or non-Gaussian noise [33] (uniform, multimodal, unsymmetrical noise distributions for instance).
In the sequel, two situations are considered: The measurement noise is described by a uniform distribution. An outlier is simulated at a given time instant. 5.1. Test 1: uniform measurement noise The first example considers a situation where the measurements are corrupted by noises described by uniform distributions (on given intervals). A typical example is the propagation errors in the ionosphere and troposphere modifying the measurement of a GNSS receiver. An extended Kalman filter is tested using the same set of measurements (given by intervals) as in Section 4.2.3 (Fig. 3). To describe the covariance matrix of the measurement noise (used as a tuning parameter in the EKF), a standard deviation is chosen for each measured variable in such a way that the 99.99%-confidence intervals computed with a Gaussian distribution are in the range of the measurement intervals given by the sensors according to a uniform distribution: 0.65 m for the position, 0.13 m/s for the velocity and 0.013 m/s2 for the acceleration. As in the previous example, three accelerometers are considered. Fig. 4 shows the estimation results with 99.99%-confidence intervals for the estimation errors. As compared to the confidence interval observer (Fig. 3), Kalman filtering generally provides narrower intervals, which is interesting at first sight. However these narrower intervals do not always contain the true values, which is detrimental in security applications. For instance, 8 position
Fig. 4. Estimation with an extended Kalman filter: error intervals. Red: lower and upper bounds. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
549
Fig. 5. Effect of an outlier on a position measurement. Left: EKF, right: confidence interval observer.
Fig. 6. Effect of an outlier on a velocity measurement. Left: EKF, right: confidence interval observer.
Fig. 7. Effect of an outlier on an acceleration measurement. Left: EKF, right: confidence interval observer.
intervals do not enclose the true position at a measurement instant (601 measurements in the simulation). Moreover, 10 velocity intervals and 32 acceleration intervals do not enclose the true values. To alleviate this problem the measurement noise covariance could be increased resulting in larger confidence intervals, still with no guarantee that they contain the true values.
5.2. Test 2: outlier The second example deals with the presence of an outlier. An aberrant measurement can happen according to environmental conditions or due to a sensor fault. For instance, in a GNSS receiver, the multipath effect can provide abnormal measurements because of the local geometrical configuration. Three cases are shown:
550
G. Goffaux et al. / Information Fusion 14 (2013) 541–550
An outlier is simulated for a position measurement at t = 25 s (Fig. 5). An outlier is simulated for a velocity measurement at t = 25 s (Fig. 6). An outlier is simulated for an acceleration measurement at t = 25 s (Fig. 7). The EKF is tested with measurements following a uniform distribution as in the previous test. One can notice that with the considered configuration of the confidence interval observer, recovery after the outlier is faster than with an EKF. It is explained by the fact that the measurement is selected at each measurement instant. Moreover, in the third case, the confidence interval observer is not sensitive to the outlier provided by an accelerometer thanks to the combination of 3 accelerometers. Furthermore, the presence of outliers can cause the failure of a guaranteed interval observer in the correction step. Indeed, the correction step proceeds with intersections and an outlier can lead to a void intersection. 6. Conclusion This paper deals with the design of confidence interval observers applied to vehicle positioning. The main objective is to satisfy security criteria which are critical in the railway or aeronautics area. The method uses a model-based interval predictor and a combination of intersection and union operations at each measurement instant between probabilistic intervals from hardware sensors and the interval predictor. The algorithm is compared to the extended Kalman filter. Even though Kalman filtering is a reference method in estimation, it is not directly applicable for mainly two reasons: (a) Kalman filtering assumes zero-mean white noise on the measurements and the state and does not allow uncertainties on the model parameters and system inputs to be accounted for, (b) Kalman filtering provides confidence intervals on the state estimate trajectories based on the assumption of a Gaussian distribution of the noise, which is not always verified in practice (particularly if GPS signals are used). Moreover, outliers can occur and a prior test should be performed before applying KF equations. Furthermore, the method is also compared with guaranteed interval observers, which are less robust to outliers and do not allow the consideration of security objectives specified in terms of probability. Acknowledgements This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its authors. References [1] M. Grewal, L. Weill, A. Andrews, Global Positioning Systems, Inertial Navigation and Integration, Second ed., John Wiley & Sons, Inc., New Jersey, 2007. [2] B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter – Particle Filters for Tracking Applications, Artech House, Boston, 2004. [3] T.I. Fossen, Marine Control Systems: Guidance, Navigation and Control of Ships, Rigs and Underwater Vehicles, Marine Cybernetics AS, Trondheim, 2002.
[4] A. Simsky, F. Wilms, J.-P. Franckart, GNSS-based fail-safe train positioning system for low-density traffic lines based on one-dimensional positioning algorithm, in: Proceedings of NAVITEC 2004 Workshop, Noordwijk, The Netherlands, 2004. [5] L. Jaulin, M. Kieffer, O. Didrit, E. Walter, Applied Interval Analysis, Springer, London, England, 2001. [6] M. Kieffer, E. Walter, Guaranteed nonlinear state estimator for cooperative systems, Numerical Algorithms 37 (1–4) (2004) 187–198. [7] M. Kieffer, E. Walter, Guaranteed nonlinear state estimation for continuoustime dynamical models from discrete-time measurements, in: Proceedings of the 5th IFAC Symposium on Robust Control Design, 2006. [8] L. Jaulin, Robust set-membership state estimation; application to underwater robotics, Automatica 45 (1) (2009) 202–206. [9] N. Ramdani, N. Meslem, Y. Candau, A hybrid bounding method for computing an over-approximation for the reachable set of uncertain nonlinear systems, IEEE Transactions on Automatic Control 54 (10) (2009) 2352–2364. [10] T. Raı¨ssi, G. Videau, A. Zolghadri, Interval observer design for consistency checks of nonlinear continuous-time systems, Automatica 46 (3) (2010) 518– 527. [11] O. Bernard, J.-L. Gouzé, Closed loop observers bundle for uncertain biotechnological models, Journal of Process Control 14 (7) (2004) 765–774. [12] M. Moisan, O. Bernard, J.-L. Gouzé, Near optimal interval observers bundle for uncertain bioreactors, Automatica 45 (1) (2009) 291–295. [13] G. Goffaux, A. Vande Wouwer, O. Bernard, Improving continuous–discrete interval observers with application to microalgae-based bioprocesses, Journal of Process Control 19 (7) (2009) 1182–1190. [14] N. Meslem, N. Ramdani, Y. Candau, Using hybrid automata for set-membership state estimation with uncertain nonlinear continuous-time systems, Journal of Process Control 20 (4) (2010) 481–489. [15] O. Bilenne, Integrity-directed sequential state estimation: assessing high reliability requirements via safe confidence intervals, Information Fusion 8 (1) (2007) 40–55. [16] E. Druet, O. Bilenne, M. Massar, F. Meers, F. Lefèvre, Interval-based state estimation for safe train positioning, in: 7th World Congress on Railway Research, Montreal, Canada, 2006. [17] Y. Zhu, B. Li, Optimal interval estimation fusion based on sensor interval estimates with confidence degrees, Automatica 42 (1) (2006) 101–108. [18] D.J. Berleant, M. Ceberio, G. Xiang, V. Kreinovich, Towards adding probabilities and correlations to interval computations, International Journal of Approximate Reasoning 46 (3) (2007) 499–510. [19] D.J. Berleant, O. Kosheleva, V. Kreinovich, H.T. Nguyen, Unimodality, independence lead to NP-hardness of interval probability problems, Reliable Computing 13 (3) (2007) 261–282. [20] G. Goffaux, A. Vande Wouwer, M. Remy, Vehicle positioning by a confidence interval observer – application to an autonomous underwater vehicle, in: Proceedings of the IEEE Intelligent Vehicles Symposium, Istanbul, Turkey, 2007. [21] O. Bilenne, Integrity-directed sensor fusion by robust risk assessment of faulttolerant interval functions, in: Proceedings of the 9th International Conference on Information Fusion, Florence, Italy, 2006. [22] M. Fréchet, Généralisation du théorème des probabilités totales, Fundamenta Mathematicae 25 (1935) 379–387. [23] M. Müller, Über das fundamentaltheorem in der theorie der gewöhnlichen differentialgleichungen, Mathematische Zeitschrift 26 (1926) 619–645. [24] H. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, American Mathematical Society, 1995. [25] M. Fruchard, O. Bernard, J.-L. Gouzé, Interval observers with guaranteed confidence levels. Application to the activated sludge process, in: Proceedings of the IFAC World Congress, Barcelona, Spain, 2002. [26] R. Kaller, J.-M. Allenbach, Traction Electrique, vol. 1, Presses Polytechniques et Universitaires Romandes, 1995. [27] C. Courtois, J. Coumel, Traction électrique ferroviaire: dynamique ferroviaire et sous-stations, Techniques de l’ingénieur (2009) 1–16. [28] M. Moisan, O. Bernard, Interval observers for non monotone systems. application to bioprocess models, in: 16th IFAC World Congress, Prague, Czech Republic, 2005. [29] F. Mazenc, O. Bernard, Interval observers for linear time-invariant systems with disturbances, Automatica 47 (1) (2011) 140–147. [30] N. Meslem, N. Ramdani, Interval observer design based on nonlinear hybridization and practical stability analysis, International Journal of Adaptive Control and Signal Processing 25 (3) (2011) 228–248. [31] R.E. Kalman, R.S. Bucy, New Results in Linear Filtering and Prediction Theory, Transactions of the ASME, Journal of Basic Engineering 83 (1) (1961) 95–108. [32] K. Reif, R. Unbehauen, The extended Kalman filter as an exponential observer for nonlinear systems, IEEE Transactions in Signal Processing 47 (8) (1999) 2324–2328. [33] E.L. Haseltine, J.B. Rawlings, Critical evaluation of extended Kalman filtering and moving-horizon estimation, Industrial & Engineering Chemistry Research 44 (8) (2005) 2451–2460.