Alexandru Cornel

Tadeo Brezoescu

Pedro Espinoza

Rogelio Castillo

Lozano

A Brezoescu

T Espinoza

P Castillo

R Lozano

Adaptive trajectory following control of a fixed-wing UAV in presence of crosswind

An adaptive backstepping approach to obtain directional control of a fixed-wing UAV in presence of unknown crosswind is developed in this paper. The dynamics of the cross track error with respect to a desired trajectory is derived from the lateral airplane equations of motion. Adaptation laws are proposed to estimate the parameters of the unknown disturbances and are employed in closed-loop system. The stability analysis is proved using Lyapunov theory. In addition, several simulations taking into account unknown wind gusts are performed to analyze the behavior and the robustness of the control scheme. A test platform has been developed in order to validate the proposed control law.

I. INTRODUCTION

Unmanned Aerial Vehicles (UAVs) represent an area of great interest in the automatic control community. The absence of the pilot renders them best suited to solve dangerous situations. However, it requires significant attention in the flight control design since the vehicle may experience large parameter variations and external disturbances. The largest use of the UAVs is within military applications but they are also used in a growing number of civil applications such as firefighting, digital mapping or monitoring. To increase the usefulness of UAVs, the capability of the autonomous controller to track a reference path is essential. Moreover, the robustness with respect to environmental disturbances must be considered. For example, small UAVs are significantly sensitive to wind since its magnitude may be comparable to the UAVs speed.

A wide range of control approaches for trajectory tracking purpose could be found in literature for both underwater and aerial vehicles. The problem of trajectory tracking for an underwater vehicle was formulated as a gain scheduling control problem in [START_REF] Silvestre | On the design of gainscheduled trajectory tracking controllers[END_REF] while Repoulias and Papadopoulos [START_REF] Repoulias | Trajectory planning and tracking control design of underactuated AUVs[END_REF] employed a method based on feedback linearization, backstepping and nonlinear damping design tools. In gain scheduling control the system dynamics are considered slowly varying [START_REF] Rugh | [END_REF] [START_REF] Shamma | Guaranteed Properties of Gain Scheduled Control of Linear Parameter-Varying Plants[END_REF] and that reduces the flight capabilities of an airplane. Trajectory linearization control (TLC) was used in [START_REF] Adami | 6DOF flight control of fixed-wing aircraft by Trajectory Linearization[END_REF] to avoid the use of gain scheduling and to enable operation across the full flight-envelope for a 6DoF fixedwing aircraft model. Nelson et al. proposed in [START_REF] Nelson | Vector field path following for small unmanned air vehicles[END_REF] a method based on the vector field approach for the case where the time dimension of the reference trajectory is removed. The algorithm was validated through simulations and real flight tests of a fixed-wing miniature air vehicle.

When accurate knowledge of the vehicle dynamics is not available, adaptive control design can be employed in order to estimate the unknown parameters. Many of the results in adaptive control are derived from Lyapunov stability theory [START_REF] Slotine | Applied Nonlinear Control[END_REF]- [START_REF] Krstic | Nonlinear and Adaptive Control Design[END_REF]. Several flight control algorithms which combines adaptation with other control tools, such as backstepping, neural networks or sliding mode control, can be found in the literature. For instance in [START_REF] Harkergard | Backstepping and Control Allocation with Applications to Flight Control[END_REF], flight control laws for two different control objectives were designed employing backstepping technique : maneuvering purpose and automatic control for the flight path angle. Also, two schemes based on adaptive backstepping and nonlinear observer design were proposed for estimating model errors. Likewise, a Lyapunovbased adaptive backstepping approach with online estimation of the uncertain aerodynamic forces and moments was used in [START_REF] Sonneveldt | Adaptive Backstepping Flight Control for Modern Fighter Aircraft[END_REF] to design a flight-path controller for a nonlinear highfidelity F-16 model. It was shown that trajectory control can still be accomplished with these uncertainties while good tracking performance is maintained. On the other hand, in [START_REF] Gavilan | Control of the longitudinal flight dynamics of an UAV using adaptive backstepping[END_REF] the authors introduced the design of an adaptive backstepping controller for longitudinal flight-path control when the aerodynamic coefficients are not known exactly. The system followed references in velocity and flight path angle and showed good performance in simulations.

Even if there are many adaptive approaches to flight control design, only few have been developed to realize airplane directional control in presence of unknown wind gusts. The goal of this work is to stabilize an airplane under crosswind and to realize the convergence to zero of the cross track error with respect to a desired trajectory. Moreover, the adaptive controller must be robust, by construction, with respect to external and unknown disturbances. We focus mainly in the lateral dynamic of the plane, for this, an analysis of this dynamic is presented in section II. Likewise, in this section we introduce the cross track error and the dynamic velocity of the plane with respect to the desired path. An adaptive control strategy is developed and presented in section III in order to follow the trajectory in presence of wind. Besides, the stability properties of the controller are discussed at the end of this section. The validation of the proposed control scheme is done in simulations and the main results are depicted in graphs in section IV. Additionally, a prototype of the airplane was developed to validate the control algorithm, the main characteristics of this prototype are described in section V. And finally in section VI, the conclusion and future work are discussed.

II. AIRCRAFT MODEL DESCRIPTION

The dynamic characteristics of an airplane strongly depend on many parameters such as altitude, speed, configuration or environmental disturbances. As a result, its complete dynamic is nonlinear, coupled and complex for control purposes. In this study we focus mainly in the trajectory following problem of an airplane flying in level flight in presence of crosswind. This kind of flight occurs when lift equals weight, thrust equals drag and the airplane flies at constant velocity, and this implies that the roll and pitch angle are small. In addition, level flight implies zero flight path angle, γ, and positive angle of attack, α, see [START_REF] Etkin | Dynamics of Atmospheric Flight[END_REF]. Therefore, using the relationship of the latter to pitch angle θ α = θγ it can be concluded that θ = α. Consequently, the airplane velocity, the angle of attack and, implicitly, the pitch angle vary slowly compared to the other parameters, and their time derivatives can be neglected in the flight dynamic.

In real conditions an airplane is generally exposed to crosswind. We define a crosswind like a wind which occurs perpendicular to the vehicle but parallel to the ground. If a plane is experiencing a crosswind, it will be pushed over or yawed away from the wind.

In order to derive the equations of motion two reference frames are used : the Earth-fixed frame, denoted by F E , and the body-fixed frame represented by F B . F E and F B have two dimensions since only the lateral dynamics are considered in this study. The origin of F B coincides with the vehicle's center of mass and the direction of its axes is according to Fig. 1. F E is employed as an inertial frame since the rotation of the Earth is neglected. Its origin is denoted by O E while O E x E is chosen northwards and O E y E points east. The classical relation of the aircraft velocity relative to the Earth is

V = V + W (1)
where

V = [u v]
T denotes the velocity of the aircraft relative to the local atmosphere and W represents the wind relative to F E . Besides, we only consider in this study a lateral wind having North, W N , and East velocity components, W E . Denote VB = u E v E T as the velocity of the aircraft relative to the Earth in the directions of the body frame axes. Thus, from (1) it follows

u E v E = u v + B B W N W E
where B B defines the complete transformation from F E to F B assuming constant pitch angle and it is given by

B B = c θ c ψ c θ s ψ s φ s θ c ψ -c φ s ψ s φ s θ s ψ + c φ c ψ
where s θ and c θ denote sin(θ) and cos(θ), respectively.

Then, the differential equations for the coordinates of the flight path in

F E are ẋ ẏ = B T B VB or ẋ = u E c θ c ψ + v E s φ s θ c ψ -v E c φ s ψ ẏ = u E c θ s ψ + v E s φ s θ s ψ + v E c φ c ψ with u E = u + W N c θ c ψ + W E c θ s ψ v E = v + W N s φ s θ c ψ -W N c φ s ψ + W E s φ s θ s ψ + W E c φ c ψ
where x and y represent the inertial position in the x-axis (North) and in the y-axis (East). Remember that the pitch and roll angles are small so that sin(θ, φ) ≈ 0 and cos(θ, φ) ≈ 1. Moreover, considering a symmetrical airplane with a rigid spinning rotor placed in the front of its body, it can then be considered, without loss of generality, V acting only in the x-axis, see Fig. 1. Hence, the following expression can be stated

v << 1 u ≈ V and consequently ẋ = V cos ψ + ω cos ψ ω ẏ = V sin ψ + ω sin ψ ω
where ω cos ψ ω = W N , ω sin ψ ω = W E , ω is the wind velocity and ψ ω describes the wind direction.

Notice that the above equations are relatively proportional to the variation of the yaw angle, and it can be controlled using the rudder deflection of the airplane. The differential equations describing this dynamics are

ψ ≈ r ṙ ≈ cτ ψ
where r stands for yaw rate, τ ψ represents the yawing moment and c is a constant related to the aircraft moment of inertia.

An analysis of the nonlinear model is presented for two different flight conditions (with and without wind). Notice from Fig. 2 the behavior of the plane when it flies in stable or moving atmosphere. The desired trajectory is plotted in thick dashed line while the solid path describes the real airplane trajectory. The crosswind has North and East velocity components of W N = -3 m/s and W E = 5 m/s, respectively. The airplane velocity relative to the surrounding air mass is 20 m/s.

III. CONTROL DESIGN

The main control objective is to obtain directional control in order to follow a desired trajectory even in presence of unknown crosswind. To simplify the analysis, let assume that the desired trajectory is aligned with the North axis of the reference frame, then, the desired path angle, ψ d , is equal to zero. Therefore, the amount of the trajectory deviation will depend on the velocity of the airplane and wind and also on the angle of the wind in relation to the airplane. In addition we consider, for control design, that the wind velocity changes slowly such that it can be considered quasiconstant. However, it will be proved in simulations that the closed-loop system remains stable even with no constant wind.

Thus, without loss of generality, the airplane dynamics for trajectory following purpose can be defined as

ḋ ≡ ẏ = V sin ψ + k ω (2) ψ = r (3) ṙ = cτ ψ (4)
where k ω = ω sin(ψ ω ) is considered quasi-constant and due to the wind perturbation, and d is the cross track error from the desired trajectory.

To stabilize the system resulted in (2), ( 3) and ( 4), the control law will be constructed using the Adaptive Backstepping approach. Then, we define the following error variable

e 1 = d -d min (5)
where d min is the minimum constant distance from the desired trajectory. Thus,

ė1 = V sin ψ + k ω (6)
A. Convergence of e 1 to zero Propose the following positive function

V L1 = 1 2 e 2 1 thus VL1 = e 1 (V sin ψ + k ω )
To stabilize e 1 we introduce ψ v as a virtual control in the following form

V sin ψ v = -c 1 e 1 -kω1
where kω1 is the estimate of k ω and c 1 > 0 is a constant. Evaluating VL1 when ψ → ψ v it follows that

VL1 | ψ=ψ v = -c 1 e 2
1 + e 1 kω1 where kω1 = k ω -kω1 . Notice from the above equation that if kω1 → k ω then VL1 ≤ 0. Thus, rewriting V L1 , it yields

V L1 = 1 2 e 2 1 + 1 γ 1 k2 ω1
where γ 1 > 0 is a constant adaptation gain. Then

VL1 | ψ=ψ v = -c 1 e 2 1 + e 1 - kω1
γ 1 kω1
Choosing the update law as kω1 = γ 1 e 1

It follows that

VL1 | ψ=ψ v = -c 1 e 2 1 B. Convergence of ψ to ψ v Define the error e 2 = V sin ψ -V sin ψ v = V sin ψ + c 1 e 1 + kω1 (8)
and rewrite [START_REF] Nelson | Vector field path following for small unmanned air vehicles[END_REF] in terms of e 1 and e 2 ė1 = e 2c 1 e 1 + kω1

This implies that

ė2 = V r cos ψ + γ 1 -c 2 1 e 1 + c 1 e 2 + c 1 kω1 (10)
Notice that cos ψ = 1 -(sin ψ) 2 . From ( 8)

sin ψ = e 2 -c 1 e 1 -kω1 V and assuming that -π 2 < ψ < π 2 it follows that (10) becomes ė2 = rR + γ 1 -c 2 1 e 1 + c 1 e 2 + c 1 kω1 (11)
with

R = V 2 -e 2 -c 1 e 1 -kω1 2 .
Introduce the following positive function

V L2 = V L1 + 1 2 e 2 2 = 1 2 e 2 1 + 1 γ 1 k2 ω1 + e 2 2
From ( 7), ( 9) and ( 11) the derivative reads VL2 = -c 1 e 2 1 + e 2 c 1 e 2 + e 1 (γ 1 + 1c 2 1 ) + c 1 kω1 + rR By selecting the virtual control as

r v R = -e 2 (c 1 + c 2 ) -e 1 (γ 1 + 1 -c 2 1 ) -c 1 ( kω2 + kω1 ) VL2 becomes when r → r v VL2 | r=r v = -c 1 e 2 1 -c 2 e 2 2 + c 1 e 2 kω2
where kω2 = k ω -kω2 , kω2 represents a new estimate for k ω and c 2 denotes a positive constant gain. Notice that if we had employed the existing estimate kω1 , we would have had no design freedom left to cancel the unknown parameter from VL2 . Additionally, kω2 could be seen as a factor correction for kω1 .

Notice from the above equation that if kω2 → k ω then VL2 ≤ 0. Thus, rewriting V L2 , it yields

V L2 = V L1 + 1 2 e 2 2 + 1 γ 2 k2
ω2 with γ 2 > 0 and constant. Hence VL2 becomes

VL2 | r=r v = c 1 e 2 1 -c 2 e 2 2 + kω2 c 1 e 2 - kω2 γ 2
Proposing the update law kω2 = γ 2 c 1 e 2 then, it follows

VL2 | r=r v = c 1 e 2 1 -c 2 e 2 2
C. Convergence of r to r v Let us define the third error variable

e 3 = rR -r v R = rR + L 2 e 2 + L 1 e 1 + c 1 ( kω2 -kω1 ) (12)
where

L 1 = 1 -c 2 1 + γ 1 , L 2 = c 1 + c 2 .
Rewriting the error system representation, we obtain

ė1 ė2 = -c 1 1 -1 -c 2 e 1 e 2
+ kω1

e 3 + c 1 kω2
thus, the derivative of e 3 yields

ė3 =cτ ψ R - r(e 2 -c 1 e 1 -kω1 )(e 3 -L 2 e 2 -L 1 e 1 ) R - r(e 2 -c 1 e 1 -kω1 )(c 1 kω1 -c 1 kω2 ) R + L 2 e 3 + L 3 e 2 + L 4 e 1 + c 1 L 2 kω2 + L 1 kω1 with L 3 = -c 1 c 2 -c 2 1 -c 2 2 + 1 + γ 1 + c 2 1 γ 2 and L 4 = -2c 1 -c 2 + c 3 1 -2c 1 γ 1 .
Finally, introduce the following Lyapunov function

V L = 1 2 e 2 1 + 1 γ 1 k2 ω1 + e 2 2 + 1 γ 2 k2 ω2 + e 2 3 then VL = -c 1 e 2 1 -c 2 e 2 2 + e 3 ( ė3 + e 2 ) (13)
Propose the control input as

cτ ψ = - e 3 (L 2 + c 3 ) + e 2 (L 3 + 1 -r 2 ) + e 1 (L 4 + c 1 r 2 ) R - kω3 (L 1 + c 1 L 2 ) -kω2 c 1 L 2 -kω1 (L 1 -r 2 )
R where kω3 = k ω -kω3 and c 3 is a positive constant gain. Notice that the unknown term k ω appears again in VL , thus we propose a correction factor in order to realize the convergence of the states.

Introducing the above into (13), we have

VL = -c 1 e 2 1 -c 2 e 2 2 -c 3 e 2 3 + e 3 (L 1 + c 1 L 2 ) kω3
Observe that VL ≤ 0 if kω3 → k ω . Therefore augmenting V L , it yields

V L = 1 2 e 2 1 + 1 γ 1 k2 ω1 + e 2 2 + 1 γ 2 k2 ω2 + e 2 3 + 1 γ 3 k2 ω3 and VL = -c 1 e 2 1 -c 2 e 2 2 -c 3 e 2 3 + kω3 e 3 (L 1 + c 1 L 2 ) - kω3 γ 3 Choosing kω3 = γ 3 (L 1 + c 1 L 2 )e 3 VL becomes VL = -c 1 e 2 1 -c 2 e 2 2 -c 3 e 2 3 ( 14
)
The error representation of the closed-loop adaptive system is summarized below

  ė1 ė2 ė3   =   -c 1 1 0 -1 -c 2 1 0 -1 -c 3     e 1 e 2 e 3   +   kω1 c 1 kω2 L 5 kω3       kω1 kω2 kω3     =   γ 1 0 0 0 c 1 γ 2 0 0 0 L 5 γ 3     e 1 e 2 e 3   ( 15
)
where

L 5 = c 1 c 2 + γ 1 + 1.
Rewriting the control input cτ ψ in terms of d, ψ, r we have

cτ ψ = tan ψ(r 2 -L 6 ) -L 7 r - L 8 d + L 9 kω1 + L 10 kω2 + L 11 kω3 V cos ψ (16)
with the updated parameters

kω1 = γ 1 d kω2 = γ 2 c 1 V sin ψ + c 1 d + kω1 kω3 = γ 3 L 11 V [r cos ψ + L 2 sin ψ] + +γ 3 L 11 dL 11 + c 1 kω2 + c 2 kω1
where

L 6 = 1 + L 2 c 3 + L 2 2 + L 3 L 7 = L 2 + c 3 L 8 = L 7 (L 1 + c 1 L 2 ) + c 1 (L 3 + 1) + L 4 L 9 = 1 -c 1 L 7 + L 3 -L 1 + L 2 L 7 L 10 = c 1 L 7 -c 1 L 2 L 11 = L 1 + c 1 L 2
Notice from ( 14) that VL ≤ 0 and it estates the global stability of the equilibrium (e i , kωi ) =(0, 0). From the LaSalle-Yoshizawa theorem, we have that e i and kωi ; i = 1, 2, 3; are bounded and go to zero as t → ∞. From (5) it follows that d → d min . ( 8) implies that kω1 is also bounded and

lim t→∞ ψ = arcsin - kω1 V (17)
Observe that from (12) r is bounded and r → 0. On the other hand, from (16) it follows that cτ ψ is bounded.

LaSalle's invariance principle assures that the state (e i , kωi ) converges to the largest invariant set M contained in {(e 1 , e 2 , e 3 , kω1 , kω2 , kω3 ) ∈ R 6 | VL = 0}. On this invariant set, we have e i ≡ 0 and ėi ≡ 0. From (15) it yields kωi = 0 and kωi = 0. Thus, the largest invariant set M is

M ={(e i , kωi ) ∈ R 6 |e i = 0, kωi = 0} ={(d, ψ, r, kω1 , kω2 , kω3 ) ∈ R 6 |(d, ψ, r, kω1 , kω2 , kω3 ) = (0, arcsin(- kω1 V ), 0, k ω , k ω , k ω )}
The manifold M is the single point d = 0, ψ = arcsin(-kω 1 V ), r = 0, kωi = k ω for i = 1, 2 and 3, which is globally asymptotically stable.

IV. SIMULATION RESULTS

The proposed control strategy was validated in closedloop system in simulations with various wind conditions. For simplicity we consider that the desired trajectory is aligned with the North axis of the inertial frame which makes the desired path angle ψ d = 0 • . In addition, the airplane is flying with a constant speed equal to 20 m/s and the crosswind has a direction West-East perpendicular to the desired path. For a smoother convergence we have used the following parameters in simulations :

c 1 = c 3 = 1.5 ; c 2 = 1.3 ; γ 1 = 1 ; γ 2 = 1.1 ; γ 3 = 1.4.

A. Case constant wind

Several simulations were performed to validate the controller and representative results are presented. The first simulations were carried out with a constant wind velocity of 7 m/s. The initial conditions are : d = 2 m/s, ψ = -10 • and r = 0 rad/s. For comparative control purpose, a standard nonlinear backstepping algorithm was developed to control the system (2)-( 4) and it is given by

cτ ψ b = -3r + tan ψ(r 2 -5) - 3d + 5k ω V cos ψ (18)
In Fig. 3 we show the time evolution of the aircraft deviation from the desired trajectory for constant wind when employing the controllers ( 16) and ( 18). The wind parameter, denoted by k ω , is not known and therefore considered zero in the simulation. Notice from this figure that the controller developed in ( 16) is able to provide cross track error regulation due to the adaptation laws presented in (15). For this case, the closed-loop adaptive system shows good response even in presence of unknown disturbance. Fig. 4 reveals the fact that to maintain alignment with the desired trajectory during a crosswind flight requires the controller to fly the airplane at a sideslip angle. Indeed, when the position error converges to zero, the yaw angle is stabilized around a constant value and the airplane keeps moving toward North. Notice that the yaw angle is nonzero unless the atmosphere is at rest. On the other hand, the proposed adaptation scheme guarantees the convergence of the unknown parameter estimates towards its true constant value, see Fig. 5. The Lyapunov function, plotted in Fig. 7, is semi-positive definite and continually decreasing which proves the stability properties of the system. Indeed, in Fig. 6 we illustrate the control input response.

B. Case variable wind gust

In order to demonstrate the robustness of the proposed control algorithm, some variations in the wind parameters are added. For this purpose, we assume that the wind varies in magnitude as shown in Fig. 8. Notice that, at time 20s, a sudden increase of 2 m/s is presented in speed of the wind. The initial conditions are the same as those for constant wind.

The main results are displayed in figures 9 -12. The wind deviates the airplane from the reference trajectory toward the wind direction but the controller ( 16) is able to recover the aircraft and to converge the position error to zero. The adaptation laws show relatively small convergence time and the estimated wind velocities are in agreement with the real values. When aligned with the reference trajectory, the airplane is flown at a sideslip angle to maintain directional control. The airplane used is the Multiplex TwinStar II whose technical characteristics are given in Table I. Its configuration is based on the classic aerodynamic layout and it is made of molded Elapor foam. Two brushless motors were mounted on the airfoil-shaped wings to power the airplane. A couple of ailerons, an elevator and a rudder are used as control surfaces and servo motors are attached to them as control surface actuators. The central processing unit, represented by the RabbitCore RCM4300 Microprocessor, collects the measurements of the IMU (Inertial Measurement Unit employed to estimate the airplane attitude and angular rates), of the airspeed sensor and of the GPS system, to compute the control law. The control responses are send to the servo signal generator/receiver unit and also to the two electric speed controllers to activate the brushless motors. Indeed, a modem is added to send and receive data from a base station. The developed prototype is presented in Fig. 13.

VI. CONCLUSIONS AND FUTURE WORK

An adaptive control algorithm based on the backstepping approach has been proposed in this paper. The control strategy was focused on reducing the position deviation of the airplane with respect to a desired path in the lateral dynamics in presence of unknown wind. The control scheme was derived considering adaptation laws to estimate the unknown wind parameters.

The closed-loop system was evaluated in several simulations and the main results, showing the good performance, were introduced by some graphs. An embedded control

Fig. 1 .

 1 Fig. 1. Tracking formulation problem

Fig. 2 .

 2 Fig.2. Earth-Relative Aircraft Location. First, the plane flies in stable atmosphere and it is capable to follow the desired path. When the atmosphere moves relative to the Earth, the airplane diverge from the path.

Fig. 3 .

 3 Fig. 3. Position Error for unknown wind. Solid line represents the proposed controller (16) whilst dashed line the standard backstepping control algorithm (18).

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 4567 Fig. 4. Yaw angle for kω = 7 m/s

  Fig. 8. Variable wind gust

Fig. 9 .

 9 Fig. 9. Position error for variable wind

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Yaw angle for variable wind

Fig. 12 .

 12 Fig. 12. Control input for variable wind

Fig. 13 .

 13 Fig. 13. The airplane prototype

TABLE I PARAMETERS

 I OF THE AIRPLANE.

	Parameter	Value
	Wingspan	1420 mm
	Fuselage length	1085 mm
	Wing area	43 dm 2
	Weight approx.	1500 g
	Wing loading	35 g / dm 2
	RC functions	Aileron, elevator, rudder, throttle

This work was supported in part by the DGA (French Arms Procurement Agency of the Ministry of Defence) and in part by the CNRS (National Center for Scientific Research)

A. Brezoescu,


A Hybrid Positioning Method Based on Hypothesis Testing Nicolas Amiot, Troels Pedersen, Mohamed Laaraiedh, and Bernard Uguen Abstract-We consider positioning in the scenario where only two reliable range estimates, and few less reliable power observations are available. Such situations are difficult to handle with numerical maximum likelihood methods which require a very accurate initialization to avoid being stuck into local maxima. We propose to first estimate the support region of the two peaks of the likelihood function using a set membership method, and then decide between the two regions using a rule based on the less reliable observations. Monte Carlo simulations show that the performance of the proposed method in terms of outlier rate and root mean squared error approaches that of maximum likelihood when only few additional power observations are available.

Index Terms-Location estimation, decision theory, estimation theory, time of arrival, received signal strength, set membership methods, interval analysis.

I. INTRODUCTION

W E consider the scenario where only two reliable range estimates and few less reliable observations are available. This situation occurs when only few links can provide time-of-arrival range estimates, but where a number of power observations can be achieved from elsewhere. In such scenarios, range-based positioning generally offers high accuracy, especially if considering ultra wideband observations. Power measurements in comparison, generally lead to inaccurate estimated ranges due to the log-normal relation between power and distance. The hybrid position estimator should thus fuse heterogeneous observations of very different accuracies. However, as reported in [START_REF] Laaraiedh | Hybrid data fusion techniques for localization in UWB networks[END_REF], this fusion is non-trivial for some hybrid positioning algorithm. Indeed, the introduction of additional, but less informative, power observations when a few accurate range observations are already available, may in fact lower the positioning accuracy. Thus, to take full advantage of the power information, suitable hybrid positioning algorithms are needed.

For the particular problem at hand, when considering only two reliable range estimates, the likelihood function can be dominated by two narrow peaks at the exact same height, leading to an ambiguous maximum likelihood (ML) estimator. This ambiguity is addressed in analysis of flip ambiguity problems, see e.g. in [START_REF] Kannan | Analysis of flip ambiguities for robust sensor network localization[END_REF]. Introducing extra power observations largely changes the heights, but not the widths, of these peaks. Albeit the ML estimation is in this case unambiguous, the required global numerical optimization of a likelihood Manuscript received February 21, 2012. The associate editor coordinating the review of this letter and approving it for publication was I. Guvenc.

N. Amiot, M. Laaraiedh, and B. Uguen are with IETR, University of Rennes 1, Rennes, France (e-mail: {nicolas.amiot, mohamed.laaraiedh, bernard.uguen}@univ-rennes1.fr).

T. Pedersen is with the Section Navigation and Communications, Department of Electronic Systems Aalborg University, Aalborg, Denmark (e-mail: troels@es.aau.dk).

Digital Object Identifier 10.1109/WCL.2012.051712.120131 function with multiple narrow peaks renders the ML too computationally demanding for many practical applications. The use of a local optimizer can be considered, but this optimization requires a good initialization to avoid to be trapped in a local maximum. An alternative approach is to solve the problem in two steps: First estimate the support regions of the peaks of the likelihood function, then calculate the final position estimate. The problem of estimating the regions can be solved among others by set-membership methods [START_REF] Jaulin | Applied Interval Analysis: With Examples in Parameter and State Estimation[END_REF], [START_REF] Meizel | Initial localization by set inversion[END_REF]. In a set-membership method, each observation defines a subset in space of possible position estimates. The support region is then be computed by intersecting these subsets by using, e.g. RSIVIA algorithm [START_REF] Drevelle | A set-membership approach for high integrity height-aided satellite positioning[END_REF]. In a heterogeneous positioning problem, as considered here, the high accuracy observations lead to small subsets, whereas the less informative observations lead to larger subsets, which include the smaller subsets. Thus, the support regions estimated by intersection of the small subsets alone do not shrink further by intersecting with the larger subsets. Consequently, in the problem at hand, the introduction of low accuracy information, such as the power observations, neither improves nor degrades the performance of the algorithm. In the case where only two range estimates are available, the set membership methods return one or two disjoint subsets corresponding to the supports of the peaks of the likelihood function. If the algorithm returns two subsets, a positioning ambiguity arises.

In the present contribution, we propose a method to obtain the final position estimate with the use of the two support regions returned by a set membership approach. We formulate the selection of the two disjoint subsets as a standard hypothesis test based on the less informative power observations. The final position estimate is afterward obtained as the centroid of the chosen subset. Simulations show that the performance of the proposed method are close to that of the ML estimator.

II. POSITION ESTIMATION BASED ON A DECISION CRITERION

A. Description of Scenario

The considered scenario is illustrated in Fig. 1. The position B of a blind node is estimated from the two noisy range observations r 1 , r 2 provided by the range nodes at known positions R 1 and R 2 ,

r i = R i -B + δ i ,i =1, 2, (1)
where δ i is the error in the range estimate. Given a probability model for δ i , it is possible to determine a confidence interval for the range estimate r i ,w h i c ha ss h o w ni nF i g .1 ,y i e l d s a confidence region shaped as an annulus with center R i .

The ambiguity problem occurs when the intersection of two 2162-2337/12$31.00 © 2012 IEEE Fig. 1: The blind node at position B receives two range estimates r1 and r2 from ranging nodes at position R1 and R2 respectively. Errors in range estimates are modeled by two confidence annuli centered at R1 and R2, respectively. The intersection of annuli are the two shaded regions C1 and C2, with their centroids C1 and C2 respectively. In addition, the blind node achieves a log power observation from each helping node at position

H k . The distance from H k to Ci is denoted by d k,i .
confidence annuli splits into two disconnected subsets C 1 and C 2 with centroids C 1 and C 2 , respectively. The two subsets and their centroids can be obtained using, e.g the algorithm proposed in [START_REF] Drevelle | A set-membership approach for high integrity height-aided satellite positioning[END_REF].

To provide an unambiguous position estimate, additional information is required. We assume that additional observations {P k } of the log power are available from the helping nodes at positions {H k }. We model the log power observation P k as a distance dependent mean µ(d k ) distorted by an additive error term X k , i.e.,

P k = µ(d k )+X k , (2)
with d k = B -H k . We assume the joint probability density function (pdf) of X 1 ,...,X K to be known. We gather the log power observations into a vector P =[P 1 ,...,P K ].

Practically, the log power information can be obtained from the received signal strength indicators.

B. Proposed Decision Rule

We approximate the conditional pdf for the power observation conditioned by the position of the blind node as:

f P|B,r1,r2 (p) ≈ ⎧ ⎪ ⎨ ⎪ ⎩ f P|C1,r1,r2 (p),B ∈ C 1 f P|C2,r1,r2 (p),B ∈ C 2 0, otherwise, (3)
with p ∈ R K . In (3) we neglect the probability of the event

B/ ∈ C 1 ∪ C 2 .
d k ≈ d k,1 = C 1 -H k ,B ∈ C 1 d k,2 = C 2 -H k ,B ∈ C 2 . ( 4
)
With the above approximations, the solution of the ambiguity problem can be phrased as a classical decision problem where λ is the likelihood ratio. The decision threshold γ can be defined to account with a priori information or costs [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection theory[END_REF]:

λ = f P|C1,r1,r2 (p) f P|C2,r1,r2 (p) C2 ≶ C1 γ. (5)
The ML decision rule is obtained for γ =1.

C. Special Case: Uncorrelated Gaussian Log Power Errors

In the special case where X 1 ,...,X K are independent Gaussian random variables with zero mean and variances σ 2 1 ,...,σ 2 K , yields for B in C i :

f P|C i ,r 1 ,r 2 (p)= K k=1 1 √ 2πσ k exp - (p k -µ k,i ) 2 2σ 2 k ,i=1, 2, (6)
with µ k,i = µ(d k,i ). Hence, the log likelihood ratio Λ=lnλ reads:

Λ= K k=1 (p k -µ k,2 ) 2 2σ 2 k - (p k -µ k,1 ) 2 2σ 2 k (7)
We obtain the ML decision rule upon insertion of ( 7) into ( 5) with γ =1:

K k=1 1 2σ 2 k µ 2 k,2 -µ 2 k,1 C1 ≶ C2 K k=1 1 σ k P k (µ k,1 -µ k,2 ). ( 8
)
It can be observed that for fixed centroids C 1 and C 2 ,t h e left hand terms are constants, while the right hand terms are a Gaussian random variable. Thus, the computation of error probability is well-known [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection theory[END_REF].

III. NUMERICAL ANAL Y S I SO FT H EPROPOSED METHOD

In this section the performance of the proposed method is compared to ML approaches via Monte Carlo simulations of the scenario described in Subsection III-A. We consider a true ML estimator relying on global optimization and an ML approximation (ML-WLS) in which a local optimizer is initialized with a weighted least squares solution [START_REF] Laaraiedh | Enhancing positioning accuracy through direct position estimators based on hybrid RSS data fusion[END_REF], both introduced in Subsection III-B.

A. Simulations Scenario

The performance of the proposed solution is assessed via Monte Carlo simulations. We consider the special case described in Subsection II-C for the setup given in Fig. 1 with the parameters settings as in Table I. We draw the positions B, R 1 ,R 2 independently and uniformly on the area S.T h e positions {H k } of the helping nodes are independently drawn according to an uniform distribution on the larger area L.T he range errors δ 1 and δ 2 assumed to be independent zero mean Gaussian random variables with variance σ 2 δi = σ 2 δ . The mean of the received log power is modeled according to the standard path loss model µ(d k ):

µ(d k )=P 0 -10n p log 10 (d k ), (9)
where P 0 is the power received at 1 meter and n p is the path loss exponent. The variance σ 2 k is chosen equal to σ 2 X for all k. Values for P 0 , n p , σ 2 δ and σ 2 X are chosen according to the measurements reported in [START_REF]Deliverable 4.1: Measurements of localtiondependent channel features[END_REF].

B. Maximum Likelihood Estimation

The ML estimator for the hybrid positioning problem reads B ∈ arg max z Λ Hybrid (z),w h e r eΛ Hybrid (z) denotes the log likelihood function for B based on {r i } and P. One approach is to find the maxima of the local extrema of the log likelihood function, obtained by equating the gradient to zero. For independent range estimates and power measurements, the gradient of the log likelihood function reads:

∇Λ Hybrid (z)=∇Λ Power (z)+∇Λ Range (z), (10)
where Λ Power and Λ Range are the log likelihood functions of the power measurements and of the range estimates, respectively. The two gradients read [START_REF] Laaraiedh | Enhancing positioning accuracy through direct position estimators based on hybrid RSS data fusion[END_REF], [START_REF] Laaraiedh | A maximum likelihood TOA based estimator for localization in heterogeneous networks[END_REF],

∇Λ Power (z)= K k=1 1 s 2 M k -s 2 -ln z -H k d k 2 (z -H k ), ∇Λ Range (z)= 2 i=1 1 σ 2 τ r i -z -R i r i (z -R i ), (11)
with the definitions

s = - σ X ln 10 10n p ,M k = (P 0 -P k )ln10 10n p +lnd 0 . (12)
Due to the non-linear relation (11), finding the roots of ( 10) requires global numerical optimization, which is not feasible for most applications. However, an approximate solution can be obtained by initializing a numerical local optimizer with an initial guess, e.g. a weighted least squares (ML-WLS) approach [START_REF] Laaraiedh | Enhancing positioning accuracy through direct position estimators based on hybrid RSS data fusion[END_REF].

C. Comparison of Performance

We compare the performances of the three algorithms in term of cumulative density functions (CDF), outlier rates and root mean square errors (RMSE). From the empirical CDFs shown in Fig. 2 it appears that for a low number of helping nodes, the performance of the proposed method outperforms ML-WLS and is close to that of ML. For high number of helping nodes, the proposed method and ML-WLS has similar performances, except in a large errors regime, where the proposed method prevails. These large errors are observed to be less frequent when the number of helping nodes is high. To inspect this difference, we consider the occurrence of outliers. We define an outlier as follows: if B ∈ C i , the estimate of B is called an outlier if it lies in the complement of C i .N o t e that for the proposed method, an outlier is equivalent to a decision error in [START_REF]Deliverable 4.1: Measurements of localtiondependent channel features[END_REF]. On Fig. 3 we observe that the outlier rate decreases with the number of helping nodes increases. Not surprisingly, the ML estimator yields the lowest outlier rate of the three methods. It also appears that the proposed method consistently outperforms the ML-WLS in terms of outlier rate. This is most significant when the number of helping nodes is less than four. These differences of performance are also reflected in the RMSEs reported in Fig. 4. In particular, for four or less helping nodes we observe that proposed method is close to the ML curve, compared to ML-WLS curve.

The above observations suggest that the RMSE for the proposed method approach can be attributed to two types of errors: large errors outliers due to decision error in [START_REF]Deliverable 4.1: Measurements of localtiondependent channel features[END_REF], and small errors resulting from the approximation in (4). The small errors occur since the centroids C 1 or C 2 are used to estimate the position of the blind node. Thus, we conjecture that the effect of these small errors can be reduced by improving this approximation, considering additional knowledge of the probability model for the range error δ i . This information could be included directly as a weighting function in the computation of centroids. Alternatively, the proposed method could be used to provide an initial guess for a numerical optimization of the likelihood function. We further conjecture that the outlier rate, which is equivalent to the rate of false decision in [START_REF]Deliverable 4.1: Measurements of localtiondependent channel features[END_REF], could be also reduced by improving the approximation (4).

IV. CONCLUSIONS

The proposed method yields a position estimate for situations where only two reliable range estimates are available along with a number of less informative observables, e.g. information on the received log power from other nodes.

In such a situation, the ML estimator necessitates numerical global optimization of an objective function with local maxima located at narrow peaks. The proposed method relies on an approximate ML decision rule with the hypothesis corresponding to the blind node residing in the support regions of each of the peaks of the likelihood function. The decision rule is formed using the less informative power observations. Finally, the position estimate is computed as the centroid of the selected peak's support regions. Monte Carlo simulations show that the performance of the proposed method in terms of outlier rate and root mean squared error, in a realistic scenario, approaches that of ML. This is in particular the case when only few additional power observations are available, i.e. when the errors due to outliers dominate. Furthermore, the proposed method outperforms an alternative procedure where a least squares approach provides initialization for numerical optimization of the likelihood function. Further improvement of the accuracy of the proposed method could be achieved by refining the estimate of the local maxima of the likelihood function.

Fig. 2 :

 2 Fig.2: CDF comparison of the absolute positioning error between the Cramer Rao Lower Bound (CRLB)[START_REF] Papakonstantinou | Cramer-Rao bounds for hybrid localization methods in LoS and NLoS environments[END_REF], the ML estimator (ML), the iterative ML estimator initialized with a weighted least squares solution (ML-WLS) and the proposed method with: 2(a) three helping nodes nodes, and 2(b) ten helping nodes.

Fig. 3 :

 3 Fig. 3: Outlier rate vs. number of helping nodes.

Fig. 4 :

 4 Fig. 4: RMSE vs. number of helping nodes.

TABLE I

 I B can be approximated byC 1 if B ∈ C 1 or by C 2 if B ∈ C 2 ,and thus:

		: Parameter Settings
	Parameter	Value
	S L σ δ	[-20, 20] × [-20, 20] m 2 [-80, 80] × [-80, 80] m 2 0.9m
	σX	4d B
	np	3
	P0	-40 dB

then

ACKNOWLEDGMENT

The work presented in this paper has been performed in the framework of the FP7 project ICT-248894 WHERE2 (Wireless Hybrid Enhanced Mobile Radio Estimators -Phase 2) which is funded by the European Union.


Nicolas Amiot
email: nicolas.amiot@univ-rennes1.fr

Mohamed Laaraiedh
email: mohamed.laaraiedh@univ-rennes1.fr

IETR Bernard Uguen
email: bernard.uguen@univ-rennes1.fr

Evaluation of a Geometric Positioning Algorithm for Hybrid Wireless Networks

Keywords: Localization, set membership method, interval analysis, hybrid wireless networks, range, difference of ranges, received power, maximum Likelihood

In this paper, we propose a geometric positioning method for hybrid wireless networks, based on a set membership method. Three common types of radio observables are considered for the position estimation: range, difference of ranges and received power. This paper details how to build geometric constraints from observables, and how to merge them to estimate the position. Given a realistic scenario, Monte Carlo simulation shows that the performance of the proposed method in terms of root mean squared error and cumulative density functions outperforms that of a numerically optimized maximum likelihood.

I. INTRODUCTION

Recently, the number of available radio access techniques (RATs) on a single terminal has drastically increased and thus simplified hybrid positioning. In the same time, the emergence of mechanic sensors embedded in mobile devices has provided a new source of exploitable information. One challenge for positioning applications is to properly merge all those information. In this situation, the standard algorithms based on least square or maximization of a likelihood function usually used to perform the position estimation have two major drawbacks: They are not convenient for using both radio and non-radio observables, and their linearization process makes difficult to approach non convex regions. This last limitation is especially an issue in positioning problems where non convex regions are often encountered. To cope with this problem, algorithms based on a geometric approach as set-membership and interval analysis have recently brought a solution [START_REF] Jaulin | A nonlinear set membership approach for the localization and map building of underwater robots[END_REF].

Contrary to classical algebraic methods, position estimation based on geometric method don't return a single position estimate but a set of intervals which contain the sought position. Recently used for addressing the problem of outdoor positioning, the geometrical approach has allowed to merge both GNSS observables and inertial sensors [START_REF] Drevelle | High integrity GNSS location zone characterization using interval analysis[END_REF]. The achieved positioning accuracy and the limited computation complexity have demonstrated the great interest of the method. As well, those geometric methods have advantageously show their performance for positioning in wireless sensor networks using range observables [START_REF] Mourad | Robust intervalbased localization algorithms for mobile sensor networks[END_REF] or received power observables [START_REF] Leger | Guaranteed robust distributed estimation in a network of sensors[END_REF].

Previous examples show that geometric positioning algorithms have always been envisaged with a single type of radio observable where a non-radio information can be added. To the best of our knowledge, this is the first attempt to apply a geometric method for hybrid positioning.

In this paper, the proposed geometric positioning algorithm for hybrid wireless networks is presented. The proposed algorithm allows to include the three most common radio observables: range, difference of ranges and observed power. The positioning accuracy is evaluated using Monte Carlo simulation and shows that the proposed method outperforms numerically optimized ML functions, for a given realistic scenario.

II. ASSUMED SCENARIO

The hybrid scenario of interest is described in Fig. 1. The blind node estimates its position β with the help of anchors providing three types of radio observables. The anchors at positions {A P } provide received power observations {P }, the anchors at positions {A D } provide difference of ranges observations {∆} and the anchors at positions {A R } provide range observations {r}. These three types of anchors are drawn on the edges of three different squares, thus :

{A P } ∈ H P , {A D } ∈ H D , {A R } ∈ H R with H R ⊂ H D ⊂ H P .
The blind node is assumed to search its position in H R .

III. CONSTRAINTS DESCRIPTION

The proposed geometric algorithm resolves the positioning problem by finding the region where all constraints are satisfied. In the following, a constraint designates both a simple mathematical expression which bounds a finite or infinite region of space, and its geometric representation. Note also that a constraint can be obtained from observables from very different nature such as a layout of a building [START_REF] Raspopoulos | Location-dependent information extraction forpositioning[END_REF], inertial motion data [START_REF] Drevelle | High integrity GNSS location zone characterization using interval analysis[END_REF], or radio observables. In the following only radio observables are considered. Three types of constraint are built: the range constraint, the difference of ranges constraint and the power constraint.

A. Range Constraint

A range constraint is defined from a range observation r evaluated between an anchor at position A R and the blind node at position β as such:

r = A R -β + δ r , (1)
where δ r is the error in the range estimate. Given a probability model for δ r , it is possible to determine a confidence interval for the range estimate r which yields a confidence region shaped as an annulus in two dimensions (2D) or a shell in three dimensions (3D), with center A R .

B. Difference of Ranges Constraint

Given a difference of ranges ∆ from anchors at position A D and A D ′ to the blind node at position β, a difference of range constraint can be expressed as:

∆ = A D -β -A D ′ -β + δ ∆ , (2)
where δ ∆ is the error in the difference of ranges estimate. Given a probability model for δ ∆ , it is possible to determine a confidence interval for ∆ which yields a confidence region contained between two hyperbolas in 2D or two hyperboloids in 3D with their common focal point D:

D = D if ∆ > 0 D ′ if ∆ < 0 (3)

C. Power Constraint

Building spatial constraints requires to formalize mathematically a distance information. Hence, to build the power constraint we propose to model the log received power observation according to the standard path loss model :

P = P 0 -10n p log 10 (d), (4)
where P 0 is the power received at 1 meter and n p is the path loss exponent. Thus, according to [START_REF] Laaraiedh | Hybrid data fusion techniques for localization in UWB networks[END_REF], the distance can be estimated as:

d = exp(M -S 2 ) + δ P (5)
with δ P the error in distance estimate, and M = log(10)(P0-P ) 10np and S = -log(10)σ X 10np

, where σ 2 X is the variance of the received power observation perturbation. Once this distance is obtained and given a probability model for δ P it is possible to determine a confidence interval for the distance d which yields a confidence region shaped as an annulus in 2D or a shell in 3D with center A P . Practically, the log received power information can be obtained from the received signal strength indicators.

D. Confidence Interval Determination

Due to the error in the observations, the constraint is associated to a limited region of the space. The extension of this region is proportional to the confidence interval chosen for the probability models of the considered observation. Practically, the probability model chosen for δ r , δ ∆ and δ P assumed to be zero mean Gaussian with σ 2 r , σ 2 ∆ and σ 2 P their variances respectively. Thus, we can build a constraint interval [I], and in particular, [I r ], [I ∆ ] [I P ], the constraint interval of the range, of the difference of ranges constraint and of the power constraint respectively:

[I r ] = [r -γσ r , r + γσ r ] (6) [I ∆ ] = [∆ -γσ ∆ , ∆ + γσ ∆ ] (7) [I P ] = [d -γσ P , d + γσ P ] (8) (9)
with γ, an adjustment factor. Without prior information we first set γ = 3 to ensure a 99% confidence interval for all the constraints.

IV. GEOMETRIC ALGORITHM DESCRIPTION Table I summarizes the 5 steps of the geometric algorithm. The detailed operation of this algorithm is presented step by step in this Section.

Table I: The Proposed Geometric Method: Algorithm Description 1) Build the constraints (Fig. 2), 2) Box the constraints (Fig. 3), 3) Merge the constraints to obtain a merged box (Fig. 4), 4) Approximate the merged box with a Kd-Tree algorithm to obtain an approximated region (Fig. 5), 5) Estimate the position from the approximated region (Fig. 6).

A. Build the Constraints

2:

The three constraints are built thanks to three different types of radio observables. Fig. 2 shows the three types of radio constraints built from the radio observables as described in Section III. Then, for the ease of computation, the constraints need to be boxed as shown in Fig. 3. Boxing a constraint consists in projecting the constraint interval on each axis in order to obtain a box B defined as [START_REF] Jaulin | Applied Interval Analysis[END_REF]:

B. Box the Constraints

B = [I] = [I 1 ] × [I 2 ] × . . . × [I k ], (10)
with

[I k ] = proj k ([I]),
where k is the axis dimension, set at k = 2 for a 2D problem, and proj k ([I]) returns the projection of the interval [I] on axis k. In order to find the smallest constraints intersection B m , the adjustment factor value γ of each constraint is reduced until min γ (B m ). Algorithm 1 describes the procedure.

C. Merge the Constraints

Once all constraint boxes have been obtained, they are merged as illustrated in Fig. 4. This merging step allows to obtain a box which, if one assumes the absence of biases, necessarily encloses the blind node position. Practically, this merged box B m is obtained by finding the intersection of all the boxed constraints B n of constraint n:

B m = N n B n ( 11
)
where N is the total number of constraints. Once the smallest merged box has been obtained, the blind node position needs to be approximated as shown in Fig. 5. For that purpose, this merged box is approximated by a recursive Kd-Tree algorithm (a.k.a. quadtree algorithm in 2D, or octree algorithm in 3D). The Kd-Tree algorithm puts a single box at the input and returns 2 n sets of boxes {B}, where n is the space dimension. Practically, this partitioning method consists in splitting all intervals of a box into two complementary intervals for each dimension of the box. Each box returned by the Kd-Tree algorithm is intersected with each boxed constraint. If the intersection is not void, the box is candidate for a new Kd-Tree iteration, otherwise the box is removed. This whole process is repeated until at least µ sets of enclosed boxes {B e } are obtained, whereupon the process is stopped. Algorithm 2 describes the complete procedure of the implemented Kd-Tree algorithm. for q in {Q} do if q ∩ ∀ B n then {B} = {B} + q : the box q obtained from the Kd-Tree is remained if q intersect all constraints boxes end if end for end for end while

D. Approximate the Merged Box

E. Estimate the Position

Fig. 6 illustrates the position estimation step. It consists in estimating the true position from the center of mass of the set of enclosing boxes {B e }. Fig. 6b is a screenshot of an estimated blind node position using our demonstrator. Note that in case where the set of enclosing boxes {B e } are disjoints, the position estimate would take advantage of an advanced estimation procedure based on hypothesis testing decision as described in [START_REF] Amiot | A hybrid positioning method based on hypothesis testing[END_REF].

V. RESULTS AND DISCUSSIONS

A. Simulation setup

The performance of the proposed geometric method is compared to a ML approximation and to the Cramer-Rao lower bound (CRLB) via Monte Carlo simulation based on the realistic scenario described in Section II. The ML approximation uses a Nelder-Mead simplex optimizer initialized with a weighted least square solution (ML-WLS) [START_REF] Laaraiedh | Enhancing positioning accuracy through direct position estimators based on hybrid RSS data fusion[END_REF]. Multidimensional likelihood functions corresponding to the given scenarios are described in [START_REF] Laaraiedh | Contribution on Hybrid Localization Techniques For Heterogeneous Wireless Networks[END_REF]. As mentioned in Section III-D, the perturbation of the observation of the range constraint, the difference of range constraint and the power constraint are supposed zero mean Gaussian, with their respective variances σ 2 r , σ 2 ∆ and σ 2 P . The parameter settings in Table II have been chosen compliant with the WHERE2 measurement campaign [START_REF]Deliverable 4.1: Measurements of locationdependent channel features[END_REF]. Finally, the last version of the entire framework used to perform those simulations can be obtained on the github website https://github.com/niamiot/RGPA. We compare the performances of the three algorithms in terms of cumulative density function (CDF), root mean square error (RMSE) and computation speed for the hybrid cases for four hybrid configurations:

B. Comparison of Performances

• Powers + difference of ranges, using 4 received power and 3 difference of ranges observables (Fig. 7a), • Powers + ranges, using 4 received power and 4 ranges observables (Fig. 7b), • Ranges + difference of ranges, using 4 ranges and 3 difference of ranges observables (Fig. 7c), • Full hybrid, using 4 received power, 4 ranges and 3 difference of ranges observables (Fig. 7d). The non-hybrid cases using a unique type of observation (only range, difference of ranges, or received power) are not considered here. From the empirical CDF shown in Fig. 7 (ad) it appears that the proposed geometric method prevails on ML-WLS. This increased accuracy of positioning is especially significant on Fig. 7a and Fig. 7b. Those two cases using received power observables allow a 1 m gain for all blind nodes. Other cases based only on time based observables as shown in Fig. 7c, or using all type of observables as shown in Fig. 7d, also show a better accuracy in terms of position estimation. Those results are confirmed by the RMSE values shown in Table III. The most significant improvement is observed for the hybrid scheme mixing powers and ranges.

In average, the proposed geometric method ensures a 30 cm increase of positioning accuracy for blind nodes drawn in a 20 × 20 m 2 room. Obviously these improvements come out at the cost of extra computation complexity. In spite of providing a complete complexity study, Fig. 8 shows some preliminary results based on an average of computation speed for each method. On those histograms, it can be observed that the proposed method is generally slower than the ML-WLS excepted when the received power observables and range observables are used. It also shows that the difference of ranges constraint is the worst in term of speed. A further investigation would be to improve the speed of the difference of ranges constraint. Moreover, between both methods is unfair, because the ML-WLS numerical optimization is based on an optimized compiled Fortran code, whereas the proposed geometric method is based on an interpreted code in Python. Considering that difference of implementation, a geometrical method as fast as the ML-WLS could be feasible. Moreover, the geometrical method is highly parallel and involves only elementary operations and could probably be very efficiently implemented in dedicated hardware.

VI. CONCLUSION

This paper has presented and evaluated a geometrical method for the positioning problem in case of hybrid observables. We have considered three cases of radio observables: range, difference of ranges and received power. Those radio observables are used to build constraints, which are merged to obtain a position estimate. Monte Carlo simulation have been computed in a realistic hybrid scenario. This simulation shows that the performance of the proposed geometrical method in terms of RMSE and CDF globally outperforms numerically optimized ML functions. In average, a 30 cm improvement of position accuracy has been observed for a blind node drawn randomly in a 20 × 20 m 2 room. As well, compared to the Nelder-Mead simplex optimized ML, the method shows promising results in term of computation speed, considering the current stage of its development. Our current work consists in the development of a dynamic and cooperative version of the algorithm. We also investigate the impact of additional non-radio constraints on the position estimation.

Fig. 1 :

 1 Fig. 1: A blind node at position β receives range observations {r} from anchors at positions {AR}, difference of ranges observations {∆} from anchors at positions {AD} and received power observations {P } from anchors at positions {AP }.

Fig. 3 :

 3 Fig. 3: The three constraints are boxed.

Algorithm 1

 1 iterative interval reduction γ=3: to ensure a 99% confidence interval B m = N n B n : compute the box intersection of all constraints while B m = ∅ do γ = γα : reducing γ and the confidence interval with α ∈ R + B ms = B m : save the previous value of B m B m = N n B n : compute the box intersection of all constraints end while B m = B ms : keep the smallest non void intersection of constraints.

Fig. 4 :

 4 Fig. 4: The three boxed constraints are merged. It results in a merged box which contains the true blind node position.

  approximated region is composed of a cluster of boxes.

Fig. 5 :

 5 Fig.5: The quadtree approximation allows to enclose the blind node position.

Algorithm 2

 2 interval approximation by boxes {B e } = {B m } : initialization with the merged box B m while card({B}) <= µ do for b in {B} do {Q} = KdT ree(b) : for each box b from the merged box B m , apply a Kd-Tree

Fig. 6 :

 6 Fig. 6: The position is estimated by computing the center of mass of all boxes.

Fig. 7 :

 7 Fig.7: CDFs of positioning error using the proposed geometric method, ML-WLS, and CRLB applied on hybrid positioning technique.

Fig. 8 :

 8 Fig.8: Speed computation comparison between the proposed geometric method and a ML-WLS using a numerical optimizer.

Table II :

 II Parameters settings

	Parameter	Value
	HP	[-1, 1] × [-1, 1] km 2
	HD	[-100, 100] × [-100, 100] m 2
	HR	[-10, 10] × [-10, 10] m 2
	σr	2.97 m
	σ∆	3.55 m
	σX	4.34 dB
	np	2.64
	P0	-40 dB

Table III :

 III RMSE vs Method

	Hybrid mode	Geometric ML-WLS CRLB
		(m)	(m)	(m)
	Power + diff. of ranges	2.46	2.91	2.23
	Power + ranges	2.51	2.81	1.78
	Ranges + diff. of ranges	1.68	1.93	0.96
	Full Hybrid	1.65	1.92	0.95

ACKNOWLEDGMENT

The work presented in this paper has been performed in the framework of the FP7 project ICT-248894 WHERE2 (Wireless Hybrid Enhanced Mobile Radio Estimators -Phase 2) which is funded by the European Union.


P-O Vandanjon

A Coiret

T Lorino

Viability Theory and Road safety

Keywords: viability theory, dynamic model, road safety, experimental tests, calm driving

Viability theory proposes geometric metaphors in addition to classical ordinary dierential equation analysis. This theory is applied on a vehicle cornering. The issue is to determine if a state of the system is compatible with a soft controls strategy which leads the car to the end of the corner. The viability theory oers tools which enable to answer to this question without computing in real time these controls. The issue resumes to verify that this state belongs to a set called the viability kernel. The construction and the use of the viability kernel for a vehicle system dynamic is proposed by using Support Vector Machines algorithm. The applicability of this theory is demonstrated through experimental tests. It is the rst time that the viability theory is applied to vehicle dynamics with road safety concerns.

1 Introduction ESP and ABS system try to nd optimal controls in order to track the desired trajectory. These systems are triggered in dangerous situations. Before to reach these dangerous situations, it would be interesting to warn the driver that it will no longer be possible to follow the desired trajectory with soft controls compatible with a calm driving.

It is dicult to compute in real time a control only to check if a soft control can keep the car on the road. The challenge is : how to warn the driver that there is no more soft control to maintain the car on the road without computing in real time such a control ?

The viability theory oers a theoretical frame suitable to this challenge. The viability kernel is precisely the set of vehicle states from which there exist a soft trajectory. The viability kernel is determined a priori, i.e. before driving. To warn the driver is equivalent to verify that the vehicle no longer belongs to the viability kernel.

The rst two sections set the scene by presenting the viability theory and the vehicle dynamic modeling. The next section presents the link between these two elds by describing the viability constrains for the vehicle system dynamics. From these constrains, the computation and the use of the viability kernel for the system are described in the section 5. In the last section of the paper, the applicability of this theory is demonstrated through experimental tests.

Viability theory

In this section, some aspects of the viability theory are presented. The reader can refer to the rst book on viability [START_REF] Aubin | Viability Theory[END_REF] for a more detailed presentation and to the last publication [START_REF] Aubin | Saint-Pierre Viability Theory: New Directions[END_REF] for the last developments. This paper focuses on concepts that will be used later.

Viability theory denes mathematical metaphors illustrating the evolution of a system:

• Whose dynamics is not deterministic,

• Which obeys to the viability constrains,

• Which is controlled by feedback loops.

One of the goals of the theory is to nd mechanisms for selecting regulons (or controls) to regulate the system. One of the selection mechanisms is based on the principle of inertia: the controls remain constant as the viability of system is not involved.

This theory has reached a level of abstraction that allows it to oer metaphors in very dierent elds, from economics to biology through Image Processing and cognitive sciences. Rather than presenting the theory from a general but complicated point of view, a didactic example is presented. This example will illustrate the main concepts that will be used afterwards.

A didactic example

In the following,a classical example of a controlled system with constrains on states and controls is presented.

ẋ = f (x(t), u(t)) (1)
u(t) ∈ U (x(t)), ∀t ≥ 0 (2) x(t) ∈ K, ∀t > 0 (3)
with u(t) is the control of the system.

submitted to Vehicle System Dynamics 2.1

A didactic example 2 VIABILITY THEORY U (x(t))
is the set of admissible controls. This set can depend on the state x(t). K is the set on states which verify constrains (they are called viability constrains). f (x, u) is the equation of the dynamic system.

Viability theory deals with set-valued map which means that the image of a point is a set. Starting from this point on view, the following set-valued map is dened :

F (x) = {f (x, u) | u ∈ U (x(t))} (4)
Then, the initial ordinary dierential equation ( 1) becomes a dierential inclusion :

ẋ ∈ F (x) (5) Denition 2.1 A state x 0 is viable for the system F if a control u(t) ∈ U (x(t)
exists such that the trajectory starting from x 0 controlled by u(t) through the

equation 3 remains in K. Denition 2.2 K is viable for the system F if ∀x ∈ K, x is viable. K T K (x 1 ) ∩ F (x 1 ) = ∅ f (x 1 , u 2 ) f (x 1 , u 1 ) T K (x 2 ) ∩ F (x 2 ) = ∅ f (x 2 , u 2 ) f (x 2 , u 1 ) T K (x 3 )
Figure 1: viability constrains, contingent cones, viability domain One of the issues of viability theory is to verify that K is viable without the need to solve for each point x an ordinary dierential equation. A geometric intuition is that to check that K is viable, it seems sucient to focus on the submitted to Vehicle System Dynamics In fact, while the trajectory is in the interior of K, there is no problem of viability. A problem emerges when the trajectory reaches the boundary of K. At this point, the viability is ensured if a control enabling the system to return in K exists. This intuition is illustrated in the gure 1 where viability constrains are represented by an ellipse. In this gure, the green and red arrows indicate the directions that can take the system starting from the origins of the arrows. Green directions point from the border to the interior of K, so the starting point is viable. In contrast, red directions do not return to K. If these directions are all possible commands, point to the origin of these arrows is not viable.

This geometric intuition is formalized mathematically by dening the tangent of a set. In viability, it is adopted the denition of Bouligand which introduced the concept of contingent cone. Denition 2.3 The contingent cone T K (x) from K at x is dened by :

T K (x) = {v ∈ X | lim inf h→0 + d K (x + hv) h = 0} (6)
• d K (x) is the distance between the set K and the point x. A deeper mathematical analysis of this distance is out of the scope of this paper.

The brown cone on the gure 1 is a graphical example of the contingent cone. The regulatory function R K (x) is dened as the set of viable controls, i.e. leading to a direction belonging to the contingent cone. This function is a set-valued map (the image of a point is a set). Denition 2.4 The Regulation map R K (x) from K in x is dened by :

R K (x) = {u ∈ U (x) | f (x, u) ∈ T K (x)} (7)
K is a viability domain if for all x in K, the regulation map R K (x) is not empty, i.e. a control u exists such that f (x, u) belongs to

T K (x) . Denition 2.5 K is a viability domain if ∀x ∈ K, ∃u ∈ U (x) | f (x, u) ∈ T K (x)
or shortly : K is a viability domain if

T V iab(K) (x) ∩ F (x) = ∅
Main theorems of viability are about the minimal assumptions for an equivalence between :

• K is a viability domain,

• and K is viable.

The interest of these theorems is that they replace a property on the dynamics, the viability of a system, by a geometrical property, the viability domain.

In many applications, the set K is not viable. Then, the viability kernel V iab(K) is dened as the set of viable states. submitted to Vehicle System Dynamics Denition 2.6 The viability kernel V iab(K) is dened by :

V iab(K) = {x 0 ∈ K | ∃u(t), x(t) ∈ K, ∀t} (8)
The anatomy of K is studied by breaking it down between the viability kernel and the complement of its kernel.

In the section 5, an algorithm for building a viability kernel is proposed. This algorithm is based on the property that a viability kernel is necessary a viability domain.

3 Vehicle dynamic modeling In this part, the model of vehicle used for the application of the viability theory to road safety is presented. The application of this theory is not specic to this model and can be adapted to others models as a more realistic 4 wheels model.

3.1 Modeling a 1 a 2 -→ X v -→ Y v -→ Z v -→ X r -→ Y r -→ Z r θ β -→ X r δ 1 δ 2 - → V α
A common model : the two-wheels vehicle model in a cornering maneuver was chosen. A not symmetric braking system is added to study the instability. submitted to Vehicle System Dynamics The gure 2 presents a top view of the model's geometry. This model comes from [3, p. 27]) . The geometric variables are the following :

• θ is the yaw angle. The angle between the x abscissa of the road frame -→ X r and the longitudinal axis of the vehicle -→ X v .

• α is the vehicle sideslip angle. The angle between the longitudinal axis of the vehicle and the velocity vector of the centre of gravity -→ V . In the following, this angle is supposed to be small (for clarity's sake, small angles are exaggerated on the gure)

• δ 1 (resp. δ 2 ) is the front (resp. rear) wheel sideslip angle. tan(δ 1 ) is the ratio of the lateral and the forward velocity of the front wheel. This angle depends on the transversal elasticity of the tyre.

• β is the front steer angle. This angle is supposed to be the same for the two front wheels. This angle is controlled by the driver.

The geometric parameters are :

• a 1 (resp. a 2 ) is the distance from the front wheels (rear wheels) to the centre of gravity G. L (= a 1 + a 2 ) is the wheelbase of the vehicle.

Model's equations are given below. The complete explanation is given in [START_REF] Ellis | Vehicle Handling Dynamics[END_REF].

The so-called small angle approximation is applied, i.e, sinus of small angle are replaced by the angle itself, cosinus of small angle are approximated by 1.

M V ( θ + α) = -D 1 δ 1 -D 2 δ 2 (9)
I z θ = -D 1 a 1 δ 1 + D 2 a 2 δ 2 - b 2 iM γ (10) (11)
• M is the mass of the vehicle ;

• D 1 and D 2 are the cornering stiness of tyres.

• V is the velocity, it is controlled by the driver according to the following equation :

V = γ ( 12
)
• γ is controlled by the driver by the braking system.

• The expression b 2 iM γ is the irregularity momentum due to longitudinal actions.

• b is the front track width.

submitted to Vehicle System Dynamics • i r is the dissymmetrical coecient. It represents the unbalance of the baking system between the right wheels and the left wheels. In the case of acceleration, this coecient is set to zero.

For clarity's sake, a simple model of tyre is presented. In fact, the cornering stinesses D 1 et D 2 vary during braking, they cannot be supposed to be constant. In the algorithm implementation (section 6), the transversal force delivered by the Pacejka model ( [START_REF] Pacejka | Magic formula tyre model[END_REF]) is used.

Cinematic considerations enable to compute δ 1 and δ 2 from the geometric variables α, β, θ and V :

δ 1 = α + a 1 θ V -β (13)
δ 2 = α - a 2 θ V ( 14
)

State space representation

In control engineering, dynamic systems are commonly described according to a state space representation. This representation is the following.

α θ = -D1+D2 M V -(1 + D1a1-D2a2 M V 2 ) -D1a1-D2a2 Izz - D1a 2 1 +D2a 2 2 V Izz α θ + -D1β M V D1βa1 Izz + biM γ Izz V = -γ (15)
The controls of the driver are β, the steering angle, γ, the acceleration (or deceleration).

The previous representation is rewritten in order to emphasize the controls u :

ẋ = A(x) • x + B • u (16)
The state x is :

x =   V α θ   (17)
The matrix A is :

A =    0 0 0 0 -D1+D2 M V -(1 + D1a1-D2a2 M V 2 ) 0 -D1a1-D2a2 Izz - D1a 2 1 +D2a 2 2 V Izz    (18)
The matrix B is :

submitted to Vehicle System Dynamics B =   0 1 -D1 M V 0 D1a1 Izz -biM Izz   (19)
It is not a linear model. The matrix A vary according to the state x : some coecients of the matrix vary with V . Moreover, the cornering stinesses D 1 and D 2 vary with x.

The controls of the driver are :

u = β γ (20) 4 Constrains
The vehicle dynamic equation 17 can be interpreted as the ordinary dierential equation 1 of the viability theory. Now, the dierent constrains on the vehicle dynamic which denes a calm driving are explored. These constrains are directly interpreted in the viability theory as :

• constrains on the states : the set K in the equation 3

• constrains on the controls : the set U (x) in the equation 2. The speed is, obviously, limited :

V min < V < V max (24)
submitted to Vehicle System Dynamics One obvious constrain is to avoid lane departure. This constrain is tackled by computing the Cartesian coordinates, x 1 , y 1 , and orientation (the yaw) θ from the Eulerian variables : α, the vehicle sideslip angle, and V , the velocity.

Then, the belonging of these Cartesian coordinates to the road is veried. For the algorithm implementation, these Cartesian and orientation coordinates are added to the state x.

Stability

Another constrain is that the grip needed by the car have to be compatible with the grip oered by the road. These constrains are directly taken into account in the Pacejka model in the sense that the car slips and does not stay on road.

A more subtle constrain is the stability of the system. According to [START_REF] Brossard | PPUR presses polytechniques[END_REF], the stability constrain when driver brakes cornering is :

L V 2 + K > 0 (25)
• L is the wheelbase ;

• K = M 1 2D1 -M 2
2D2 ,the understeer gradient ;

• M 1 (resp. M 2) is the static mass on the front (resp. rear) axle ;

M 1 = M a2 a1 + a2 ; M 2 = M a1 a1 + a2 (26)

Remarks

• The viability set K includes the following constrains : the car does not slip, and the stability is ensured. These two conditions are also required for the modeling. The proposed model is no more valid when these conditions are violated. For example, the small angle approximation is no more valid.

The viability constrains ensure that the model is still valid.

• Some constrains can be added, others can be modied or cancelled according to the objectives of the application. The viability theory takes naturally into account constrains on the state of the system.

Viability kernel

Starting from a state of the vehicle, i.e. the vector x of the equation 17, the issue is : is-it possible to drive along the corner in a calm way or, in a more mathematical sense, does one smooth trajectory exists without computing explicitly this trajectory ? The viability theory assets that the existence of this trajectory is submitted to Vehicle System Dynamics This section presents algorithms in order to set up and to use a viability kernel. Before this presentation, the system has to be discretized.

Discretization

In time

The discretization of the equation 16 gives :

x n+1 = x n + (A(x n ) • x n + B • u n ) ∆T (27)
• ∆T is the sampling time.

In space

Every dimension of the state space are discretized. The state space is 3 (see the length of x in equation 17), but 3 dimensions for the spatial localization are added, and the two dimension of the control space are discretized too. Finally, the dimension of the proper space is 8.

The implementation of the viability theory faces the dimensionality curse :

each additional variable, discretized in, for example, ten values, multiply by ten the number of states to be evaluated. This point is one of the current limit of this theory.

Algorithm Principle

The algorithm and the proof of its convergence were initially presented in [START_REF] Saint-Pierre | Approximation of the viability kernel[END_REF].

It is a backward recurrence algorithm. It starts from the end of the trajectory (called the target ). The target of a viable trajectory corresponds to the car at the end of the corner with an adapted speed, and a proper direction. The states of this target belong to the viability kernel. Then, states belonging to K from which, it is possible to reach the viability kernel with a constant control during ∆T , are added to the viability kernel. This step is applied recursively.

This algorithm is based on the link between viability kernel and viability domain which is described in the section 2.

The gure 3 illustrates this procedure in two dimension (the number of dimension is 8 for the given application). Road limits are plotted by black arc circles. Blue points belong to the viability kernel, red points do not belong to it.

The red car does not belong to the viability kernel. Indeed, driver control which leads the vehicle into the viability kernel does not exist :

• either the control brings the car outside the road, so the constrains on the state are violated, submitted to Vehicle System Dynamics For this car, softcontrols, which drive the car to the end of the corner, do not exist On the contrary, a control drives the blue car to a blue point which belongs to the viability kernel. Therefore, it exist controls which drive the car to the target in a calm way. The blue car belongs to the viability kernel.

This paper focuses on the existence of at least one control to maintain the car on the road in a way compatible with calm driving. It could be of interest to count the number of dierent strategies of soft control. Higher is this number, more robust will be the given state.

Practical use : classication of a state

Knowing a state x 0 belonging to K, the issue of this subsection is to determine if x 0 belongs to the viability kernel of K : V iab(K). V iab(K) had been discretized

with the algorithm described in the previous subsection. The previous algorithm gives two discretized sets. A set of states which belongs to the viability kernel (blue points) and a set of states which do not belong to the viability kernel (red points). To know if x 0 belongs to a set or to its complement is typically the phase of generalization for the machine learning methods.

One of the most fruitful and, relatively new, technique is this eld is the Support Vector Machines (SVM). Support Vector Machines (SVM) were issued from the Vapnik Chervonenkis theory (the initial article is [START_REF] Cortes | Support-vector networks[END_REF]). Due to their generalization performance on classication problems, SVM have obtained much success in a variety of domain. Unlike other machine learning methods, generalization errors of SVMs are not related to the input dimensionality of the problem, but to the margin with which a SVM separates data. The application of SVM to build viability kernel had been initially carried out in the ecology eld ( [START_REF] Chapel | Dening yield policies in a viability approach[END_REF] and [10]). In the following, the principles of the method are recalled.

The viability kernel algorithm, described in the previous section, gives a set of example (x i , y i ) ; where x i is a state of K and y i = 1 if x i ∈ V iab(K) and

y i = -1 if x i / ∈ V iab(K).
d(x) = n i=1 α i y i k(x i , x) + b (28)
• k(y, z) is a kernel function which denes a scalar product between the projection of y and z in an other space than the initial space to which y and z belongs ;

• b is a constant.

• α i is a scalar which is non null only for the vectors on the margins of the set. These vectors are called the Support Vector.

submitted to Vehicle System Dynamics This section aims at proving the applicability of the viability theory to a real example in road safety. The scenario is the following : a car is in cornering maneuver at two dierent speeds. At the middle of the clothoid, the driver turns suddenly the wheel in order to simulate an inappropriate reaction. At low speed, this reaction still compatible with a calm driving. At high speed, this reaction is no more compatible and the driver will not be able to maintain the car on the lane with soft controls.

In the following, the geometry of the infrastructure, the test vehicle,the maneuver and the algorithm implementation are described. The last subsection is devoted to the analysis of this maneuver with viability tool.

Infrastructure

The application is done on a real corner which is illustrated by the gure 4. This corner, called the controllability curb, is a part of the Ifsttar experimental track (located in Nantes, France). It is consisted of :

• a straight alignment (its length is 50m)

• a rst clothoid (87m)

• an arc circle (Radius : 110m, its length is 134m) • a second clothoid which is the symmetric of the rst one,

submitted to Vehicle System Dynamics 6.2 Vehicle 6 APPLICATION M I z i r b a 1 a 2 D 10 D 20 (kg) (kg • m 2 ) (m) (m) (m) (N ) (
• a second straight alignment (50m)

Vehicle

The vehicle used during our trials is a sedan. It is a front wheel drive car equipped with a V12 engine coupled with a 5-speed automatic gear box. Ground link is composed by a pseudo Mac Pherson at the front and a exible H-form torsion beam with trailing arms at the rear. The tires are classical summer tires.

Parameters of the vehicle needed by the algorithm are given in the table 1.

D 10 and D 20 are the nominal tyre lateral stiness. These gures are given in order to have a physical understanding of the system.

Maneuvers

1. the rst maneuver consist in driving along the controllability curb at a constant speed of 60 km/h (37 mph) ;

2. the second maneuver is based on the rst one, except that the driver turns the wheel suddenly along the rst clothoid ; after this action, he has to maintain his vehicle in the lane ;

3. the third maneuver consist in driving along the controllability curb at a constant speed of 80 km/h (50 mph); 4. the fourth maneuver is a repetition of the second one at a speed of 80 km/h. In particular, the driver has to turn the wheel with the same angle.

The gure 5 illustrates the state of the test vehicle during these maneuvers.

Algorithm implementation

The code has been written with the Free and Open Source software scilab ( [START_REF]Enterprises Scilab: Free and Open Source software for numerical computation[END_REF]).

Viability kernel construction

The algorithm described in the subsection 5.2 has been applied to the controllability curb. The direct constrains are the followings :

• the minimal velocity is 60 km/h, the maximal velocity is 90 km/h ; submitted to Vehicle System Dynamics maneuver 1 : speed 60 km/h maneuver 2 : speed 60 km/h and sudden turn maneuver 2 : speed 80 km/h maneuver 4 : speed 80 km/h and sudden turn • the acceleration has to be included between -0.1 g and 0.1 g (it is the frame of a smooth driving) ;

• the maximal steering velocity is +/-3 deg/s ;

• the maximal jerk +/-0.2 g/s .

The discretization is

• 4 m along the curvilinear abscissa ;

• 0.5 m for the transversal axis ;

• 0.2 s for the time It takes 51 mn to build the viability kernel (the computer is a classical old PC scoring a performance of 10.3 gigaops according [START_REF] Dongarra | Performance of Various Computers Using Standard Linear Equations Software[END_REF]) . 759 states belong to the viability kernel and 291471 do not belong to the viability kernel. Thus, about 3 % of the states belongs to the viability kernel. It includes small variation around nominal trajectories, i.e. trajectories at constant speed in the middle of the road. Every states with large deviation with these nominal trajectories belong to K but not to V iab(K).

This application shows that applying Viability theory on real case needs computing power. It has been one of the barriers for its application on a real road until now.

submitted to Vehicle System Dynamics 6.5

Results

APPLICATION

The algorithm can be parallelized. Therefore, by using Graphics Processing Unit (GPU), this computing time could be dramatically reduced or the addressed model could be more complex with the same computing time.

This computation has to be done a priori, i.e. before the practical use (i.e. on board) of the viability kernel.

Viability kernel classication

Among the dierent available algorithms for SVM, the algorithm described in [START_REF] Tohme | A simple and fast training algorithm soft margin Support Vector Classication[END_REF] is chosen because it is suitable to large data set problem. The kernel function is the Gaussian function which is a classical choice when the problem is nonlinear.

It takes 27 mn to get the decision function (the computer performance is 10.3 gigaops). The regularization parameters C is set up in order to get 19 examples which are misclassied (19 states of the viability kernel are considered as not belonging to the viability kernel and all the states of the complementary of the viability kernel are well classied), and to get 1482 support vectors.

It means that to know if a state belongs to the viability kernel, it is enough to compute the equation 28 with n = 1482. The results is immediate with a modern computer.

Results

According to the feeling of the driver, the maneuvers 1, 2 and 3 are compatibles with calm driving but not the maneuver 4.

During the maneuver, the instant when the driver is turning the wheel is considered.

These four situations are modeled by the states given in the table 2 for the maneuvers at 60 km/h and table 3 for those corresponding to 80 km/h. The interpretation of this table is the following.

• The rst three columns : x, y, θ are the Cartesian and orientation coordinates.

• The dierence between the two states of the same speed is, mainly, the steering angle β.

• The last column is the value of the decision function which classes a state according to the equation 28. The sign of the value indicates if the state belongs or does not belong to the viability kernel.

According to this last information, the maneuvers 1,2,3 are detected as viable by the algorithm but not the maneuver 4. These results are coherent with the driver's feeling.

As expected, it means that the tolerance to an inappropriate control in the steering angle is more important at a low speed than at a higher speed. It implies that the frontiers of the viability kernel is nearer to the nominal trajectory at high speed than at lower speed. In this example, the value of the decision function d(.) ranks the situation in a way coherent with our intuition. This suggests that this function can be analyzed as a distance to the frontiers of the viability kernel so as a ranking function. Thus, this function could be used to warn gradually the driver.

This algorithm can be implemented in roadside system camera. In this case, the dynamic model has to be simplied because the system can not access to the vehicle parameters, dierent viability kernels have to be computed according to weather conditions and viability kernels have to evolve with the road in order to take into account the road polishing.

It can also be implemented directly on the vehicle. In this case, the dynamic model can be richer by using all available information on CAN-bus. However, as the available grip between the tyre and the road cannot be known accurately by classical embedded measurement ( [START_REF] Andrieux | New results on the relation between tyreroad longitudinal stiness and maximum available grip for motor car[END_REF]), this model cannot be very accurate.

The most ecient solution is to mix the two systems. The car updates its viability kernel thanks to information delivered by road systems.

Conclusion -Perspective

A viability kernel has been built for a road-car-driver system. This set includes trajectories compatible with a calm driving style. The system is modelled as a classical two-wheels model in a cornering maneuver including braking. A decision function is build by using the Support Vector Machine (SVM) technique.

This decision function indicates if a state of the vehicle belongs or not to the viability kernel. The calculation of the viability kernel and the construction submitted to Vehicle System Dynamics of the decision function are done before to drive. During driving, the decision function verify if the state of the vehicle still belongs to the viability kernel. This operation consumes little computing time and can be computed on board.

Experimental tests were carried out to prove the applicability of the proposed method.

It is the rst time that the viability theory is applied to vehicle dynamics with road safety concern. This was possible through the evolution of computation power and recent progresses in classication techniques. Development of computation with Graphics Processing Units will make applications of this theory easier.

Beyond this immediate application, viability theory provides a new point of view on road safety: the issue is not to calculate optimal or robust controls but it is to know if theses controls exist, without calculating them in real time.

Perspectives are multiple. The decision function is used in a binary mode but it could also be analyzed as a distance to the frontier of the viability kernel and, thus, helps to deliver gradually warnings to the driver. The number of driving strategy to maintain the car on the road can be computed, this denes a system robustness. The most promising issue is the anatomy of the viability kernel. The shape of the viability kernel would dene the critical states.

Figure 2 :

 2 Figure 2: Two-wheels vehicle model in a cornering maneuver

  that the state belongs to the viability kernel. This property replaces the resolution of a dierential equation with the belonging of a state to a set.

Figure 3 :

 3 Figure 3: Algorithm principle

  From this list, SVM dene the decision function d(.) which separates the given examples :

Figure 4 :

 4 Figure 4: Controllability corner

Figure 5 :

 5 Figure 5: Photos of the dierent maneuvers

Table 1 :

 1 Vehicle parameters

	N )

Table 2 :

 2 Comparison of two states at 60 km/h

	submitted to Vehicle System Dynamics

Table 3 :

 3 Comparison of two states at 80 km/hThis rst investigation of the viability kernel shows that the study of the shape of the viability kernel is fruitful for analysing the road-driver-vehicle sys-

tem. Critical states belongs to parts of the viability kernel which are narrow in a geometrical sense to the boundary. Viability theory replaces analysis of differential equation by analysis of a set. It replaces dierential equations analysis by geometry analysis.

[10] L. Chapel, Maintenir la viabilité ou la résilience d'un système : les machines à vecteurs de support pour rompre la malédiction de la dimensionnalité ?, Université Blaise Pascal -Clermont-Ferrand II, 2007.submitted to Vehicle System Dynamics