André Benine-Neto
email: andre.benine-neto@u-bordeaux.fr

Stefano Scalzi
email: stefano.scalzi@uniroma2.it

Saïd Mammar
email: said.mammar@ibisc.univ-evry.fr

Mariana Netto
email: mariana.netto@ifsttar.fr

Benoit Lusetti
email: benoit.lusetti@ifsttar.fr

Model Reference Based Vehicle Lateral Control for Lane Departure Avoidance

Keywords: Invariant Sets, BMI: Bilinear Matrix Inequality, Driving Assistance Systems, Vehicle Lateral Control, Lane Departure Avoidance

This paper presents the design and practical implementation of a lane departure avoidance assistance for passenger vehicles based on a state feedback dynamic controller. The road curvature is taken into account as an internal model to ensure convergence of the lateral offset to zero at steady state, even in curvy roads. Lyapunov theory and Bilinear Matrix Inequalities, including bounds in the control input and constraints for poles clustering are used to minimise the reachable set of the vehicle after activation of the assistance. The proposed control strategy is simulated in CarSim environment and successfully tested on a prototype vehicle.

Introduction

Many highway accidents occur due to unintended lane departure, therefore intelligent vehicles and automated highway systems have attracted a growing attention in the last years. An overview of the research activity presented by [START_REF] Shladover | Path at 20 -history and major milestones[END_REF] and the progress done in the last decade shown by [START_REF] Urmson | From automated highways to urban challenges[END_REF]. The entirely or partially automation of driving task aims to increase safety by reducing accidents and the driver's workload.In order to prevent vehicle crash due to inattention, illness or sleepiness, vehicles have increasingly been equipped with electronic control systems that provide active safety. [START_REF] Kullack | Reflektas: lane departure prevention system based on behavioural control[END_REF] have proposed a lane departure warning system. Such systems are also available in some commercialised passenger vehicles. Other works address the problem of combining the driver and the automatic lane keeping control action, either by sharing or overriding the driving task with the driver. It can be noticed from the literature concerning the lane keeping problem that several control algorithms are employed. [START_REF] Mammar | Lane keeping automation at tire saturation[END_REF], [START_REF] Hessburg | A fuzzy rule-based controller for automotive vehicle guidance[END_REF] and [START_REF] Wu | The heterogeneous systems integration design and implementation for lane keeping on a vehicle[END_REF] make use of fuzzy rule-based controllers. Model predictive control is applied by [START_REF] Falcone | Predictive active steering control for autonomous vehicle systems[END_REF] using a nonlinear vehicle model to predict the system evolution for the development of an active steering system for autonomous vehicles. H ∞ techniques are explored in the design of an assistance system [START_REF] Raharijaona | h ∞ controller synthesis and analysis with application to lateral driving assistance[END_REF] and an automatic lane keeping system for highway scenario [START_REF] Cerone | Experimental results in vision-based lane keeping for highway vehicles[END_REF]. The use of PID is employed in many papers. In the framework of autonomous vehicles some techniques involve: a gain scheduling based proportional feedback from the lateral offset [START_REF] Broggi | The argo autonomous vehicle's vision and control systems[END_REF]; a PID based on a weighted sum of the heading error and the lateral offset with feedforward term from curvature [START_REF] Cremean | Alice: An information-rich autonomous vehicle for high-speed desert navigation[END_REF]; and the yaw angle and a nonlinear term proportional to the lateral offset are used as measurements to design the steering controller [START_REF] Thrun | Stanley: The robot that won the darpa grand challenge[END_REF]. As for driving assistance systems, [START_REF] Takahashi | Introduction of honda asv-2 (advanced safety vehicle-phase 2)[END_REF] applies a feedforward on the road curvature and feedback PD on vehicle lateral position and heading angle; and a comparison of PID, fuzzy, H ∞ , and adaptive controllers is presented by [START_REF] Chaib | H-infinity, adaptive, pid and fuzzy control: a comparison of controllers for vehicle lane keeping[END_REF]. Optimal preview control theory is explored by [START_REF] Macadam | Application of an optimal preview control for simulation of closedloop automobile driving[END_REF] and [START_REF] Kang | Design and testing of a controller for autonomous vehicle path tracking using gps/ins sensors[END_REF], who also makes use of a feedforward term. Linear and nonlinear sliding mode controllers with feedback of the lateral displacement and yaw rate are designed by [START_REF] Ackermann | Linear and nonlinear controller design for robust automatic steering[END_REF]. Control synthesis based on Lyapunov techniques are presented by [START_REF] Rosseter | Lyapunov based performance guarantees for potencial field lane keeping assistance system[END_REF] in which a potential field control structure for a vehicle with steer-by-wire is introduced . Kritayakirana and Gerdes (2012) carried out tests on an autonomous vehicle equipped with this controller. Minoiu-Enache et al. (2009a,b) proposed to compute closed-loop invariant sets for the lane keeping system by means of LMI (and BMI) optimisation based on Lyapunov stability conditions. [START_REF] Benine-Neto | Piecewise affine control for lane departure avoidance[END_REF] applied similar techniques to design lane keeping system based on piecewise affine controllers in order to take into account the nonlinear behaviour of the lateral tyre forces.

In the framework of advanced driving assistance systems (ADAS), many possibilities are shown in the literature concerning the activation strategies of lane keeping systems. [START_REF] Eidehall | A new approach to lane guidance systems[END_REF] proposes a system that perceives the intention of the driver based on the vehicle dynamics and tracks the environment (road geometry, obstacles/vehicles on current and adjacent lanes) to compute the risk of the manoeuvre using the time to lane crossing to predict a collision. [START_REF] Fujiwara | Steering assistance system for driver characteristics using gain scheduling control[END_REF] defines a time-varying parameter, in order to reduce the conflicts between the driver and the system assistance, so that the control action increases when the driver steering wheel torque decreases; [START_REF] Leelavansuk | Study on cooperative control of driver and lane-keeping assistance system[END_REF] make use of a weight (between zero and one) which multiplies the proposed LQ steering control in order to reduce the interference with the driver. In the ADAS proposed by [START_REF] Cerone | Combined automatic lane-keeping and driver's steering through a 2 dof control strategy[END_REF] and [START_REF] Rosseter | Lyapunov based performance guarantees for potencial field lane keeping assistance system[END_REF] an automatic lane keeping is combined with the driver's steering with no need of switching strategies between the driver and the lane keeping control. A feedforward and a feedback controllers are designed by [START_REF] Ishida | Development evaluation and introduction of a lane keeping assistance system[END_REF] to support the driver in order to reduce his/her workload without reducing the situation awareness and motivation of the driver. The lane keeping system proposed by [START_REF] Minoiu-Enache | Driver steering assistance for lane departure avoidance[END_REF] uses a switching strategy during diminished driving capability or driver inattentiveness to help the controller to fully replace the driver.

Such activation strategy seems to be the most adequate as according to ergonomic and psychology studies [START_REF] Hoc | Cooperation between drivers and automation: implications for safety[END_REF] as the driver usually does not have time to think and plan in emergency situations. There is only time to react in order to re-establish a safe state, therefore sharing the control with the driver is not really an advantage and in these cases it is better that the controller takes completely the responsibility during the control duration. As a consequence, this paper addresses the problem of designing a steering controller that switches between the driver and the assistance system, based on the attentiveness of the driver and the risk of lane departure during lane keeping.

The controller synthesis proposed in Minoiu-Enache et al. (2009a) is able to compute by means of Bilinear Matrix Inequalities (BMI), a minimal ellipsoidal invariant set that ensures the stability of the closed loop system at its activation, even in the presence of road curvatures (considered as disturbance) and taking into account bounds in the input control. This paper enhances these results by using a dynamic state feedback controller, which ensures zero lateral offset in steady state, even in curved road section characterised by nonzero curvature. Simulations results of the dynamic state feedback have been presented by the authors in [START_REF] Benine-Neto | Dynamic controller for lane keeping and obstacle avoidance assistance system[END_REF], whereas in this paper, the controller have been implemented in a prototype vehicle and results include practical experimentation on a test track. A second extension of the previous work consists in the improvement of the damping characteristics of the closed-loop response by additional constraints for pole clustering.

Section 2 presents the vehicle model used for the control synthesis, as well as a model for the disturbance, which is taken into account in the internal model design. In section 3 the control objectives and synthesis are presented. The activation strategy of the steering assistance system is also presented in this section. Experiments of the system on Carsim environment and on a test track are shown in section 4 confirming the performance of the controller and activation strategy as a lane keeping assistance. The practical experiments involve sharper curves than those specified in the design phase of the controller while even more severe conditions such as the test standard ISO/DIS 3888 are evaluated in Carsim environment. These tests also allow the evaluation the controller performance in situations not taken into account in the design phase, such as lateral tyre force saturation. Section 5 concludes the work providing some perspectives for future work.

Vehicle model

According to [START_REF] Ackermann | Linear and nonlinear controller design for robust automatic steering[END_REF], a single track vehicle model is considered to capture the essential vehicle lateral dynamics. It is developed under the assumption that wheel sideslip angles remain small, therefore considering a linearised model of the lateral tyre forces. The proposed controller is intended to be used in situations of unintended lane departure, in which the lateral tyre forces are not saturated. Hence, the single track vehicle model is adequate for the control design phase.

The front and rear wheels are represented each by a single wheel located at the centre of the axle and only the planar motion of the vehicle is considered, as depicted in Figure 1. The involved model variables are the vehicle sideslip angle β and the vehicle yaw rate r. The control input is the steering angle, denoted by δ f .

β ṙ = a 11 a 12 a 21 a 22 β r + b 11 b 21 δ f , (1)
with:

a 11 = - 2(c f +cr) mv , a 12 = -1 - 2(c f l f -crlr) mv 2 , b 11 = 2c f mv , a 21 = - 2(c f l f -crlr ) J , a 22 = - 2(c f l 2 f +crl 2 r ) Jv , b 21 = 2c f l f J .
(2)

The definition and the numerical values of the above parameters are given in Table 1. Two more equations have to be included in the model for lane keeping objectives. The corresponding variables are: the relative yaw angle, denoted by ψ L which is the difference between vehicle heading angle ψ and desired heading angle ψ d , as shown in Figure 1; and the lateral offset y L at a look-ahead distance l s . The equations describing the evolution of these variables are:

ψL = r -vρ, ẏL = vβ + l s r + vψ L , ( 3
)
where ρ is the road curvature, defined as ρ = 1/R, with R the curvature radius. This exogenous input is considered to be bounded by a maximum value of road curvature, |ρ| ≤ ρ max . A normalized variable w that permits a parameterization of the road curvature ρ can be defined as:

w = ρ ρ max . (4)
Considering equations ( 1) and (3) a linear model of the vehicle lateral motion can be written as:

ẋ = Āx + Bu u + Bw w ȳ = x, x ∈ R 4 , u ∈ R 1 , w ∈ R 1 , (5)
where

x ∆ = [β, r, ψ L , y L ]
T is the state vector, u is the control input corresponding to the steering angle δ f and the matrices describing the linear system are:

Ā = ⎡ ⎢ ⎢ ⎣ a 11 a 12 0 0 a 21 a 22 0 0 0 1 0 0 v l s v 0 ⎤ ⎥ ⎥ ⎦ , Bu = ⎡ ⎢ ⎢ ⎣ b 11 b 21 0 0 ⎤ ⎥ ⎥ ⎦ , Bw = ⎡ ⎢ ⎢ ⎣ 0 0 -vρ max 0 ⎤ ⎥ ⎥ ⎦ .
One of the main limitations of the linear model ( 5) is that the road adherence is considered constant (implicitly in the cornering stiffness parameters c f and c r ). Despite of this shortcoming, this simplified model is adequate for the control synthesis, whereas the validation of the control strategy is carried out on a experimental vehicle CarSim model which includes non-linearities, as shown in section-4.

In order to design a controller that is able to avoid lane departure despite road curvature, it is required to consider the exogenous input within the vehicle model, building an internal model. Such control technique among with the control objective and synthesis are detailed in the following section.

Control Strategy

In this section, the design requirements of the controllers for the lane departure avoidance system are presented. The development of an internal model in order to incorporate the perturbation (road curvature, ρ) into the vehicle model is also discussed. This contribution enables the convergence of the lateral offset to zero in steady state. Minoiu-Enache et al. (2009a) proposed a synthesis of state feedback controller by solving an optimisation problem subjected to Lyapunov-based BMI constraints, which are used to guarantee invariant set properties. As an extension to this work, part of this section is dedicated to show enhancement to this controller synthesis by considering pole clustering constraints in order improve damping characteristic of the closed loop system.

Control Objectives for Lane Departure Avoidance

In order to design a lane departure avoidance system that performs satisfactorily even in curves, it is essential to take the road curvature into account in the control synthesis. The controller is intended to operate within the linear behaviour of the lateral tyre forces. The following requirements are imposed on the control synthesis:

• The controller must ensure asymptotic stability of the closed loop system.

• The vehicle state trajectories should be bounded in order to avoid vehicle instability and excursions out of the lane.

• In order to ensure satisfactory performance, the vehicle must be driven to the centre of the lane as fast as possible without detriment of the vehicle stability.

• The vehicle must present zero lateral offset in steady state conditions, even while negotiating curves.

• The control input must be bounded in order to avoid damage of the actuators.

Based on these requirements, the design of a state feedback controller for lane departure avoidance is shown in the subsequent part.

Internal model

According to [START_REF] Francis | The internal model principle of control theory[END_REF], an adequate strategy for the design of control laws, considering parameter uncertainties or perturbations, consists in integrating the perturbations in the controller by means of an internal model, in which the disturbances are represented as a family of known functions.

Road building standards define that clothoids are to be used to smoothly interconnect straight roads to constant curvature roads [START_REF] Rajamani | Vehicle Dynamics and Control, chapter 2 -Lateral Vehicle Dynamics[END_REF]. Clothoids are based on Fresnel Integral with shape defined by a single parameter c. Considering l the distance measured along the clothoid, the road curvature is defined by :

ρ(l) = 1 c 2 l (6)
At constant speed, l(t) = vt, which implies that the road curvature increases linearly with respect to curviligne abscissa and thus to time:

ρ(t) = v c 2 t (7)
Thus, the disturbances on the model used in the design of the lane departure avoidance system consist of constant values and ramps. Noting where L is the Laplace transform operator, the transfer function L{y L }/L{ρ} of closed loop system formed by (5) and state feedback controller of the form u = K x is given by :

L{y L } L{ρ} = (s -z 1 )(s -z 2 ) (s -p 1 )(s -p 2 )(s -p 3 )(s -p 4 ) (8)
where the zeros z i and poles p i for {i = 1, 2, 3, 4} depend on the vehicle parameters. Such control action is not able to reject the disturbances due to the road curvature (such as constant values and ramps) and ensure a zero lateral offset in steady state. For that, the transfer function L{y L }/L{ρ} must present a double zero at the origin of the complex plane, which is obtained by a double integration of the lateral offset to the control input, i.e.

u = K x + k I1 t 0 y L dτ + k I2 t 0 τ 0 y L dτ dν.
As the control synthesis applied in this work (shown in details in section 3.3) is based on a state feedback controller, a state space representation of the internal model (dynamics of integrators) is more adequate. Thus the internal model consists of the states α = [α 0 , α 1 ] T , which uses the lateral offset y L as input, as follows:

α = 0 1 0 0 α + 0 y L (9)
These additional variables introduced by the internal model are included in the state vector of (5), which yields x ∆ = [β, r, ψ L , y L , α 0 , α 1 ] T , and the augmented system becomes,

ẋ = Ax + B u u + B w w y = x, x ∈ R 6 , u ∈ R 1 , w ∈ R 1 . ( 10
)
Accordingly the system matrices become:

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
a 11 a 12 0 0 0 0 a 21 a 22 0 0 0 0 0 1 0 0 00 v l s v 0 0 0 0 0 0 0 01 0 0 0 1 00

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , B u = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ b 11 b 21 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , B w = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 -vρ max 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
.

Applying in system (10) a state feedback of the form u = Kx , with K = [k 1 , k 2 , k 3 , k 4 , k 5 , k 6 ], the closed loop system becomes:

ẋ = (A + B u K)x + B w w. (11)
and its equilibrium point is:

x eq = -wρ max ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (b1a22-a12b2) (-b2a11+b1a21) v -v (b1lsa21-b1a22v-b2a11ls+b2a12v) -b2a11+b1a21 v 0 x eq5 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (12)
in which

x eq5 = -(vb 21 k 3 a 12 -b 21 k 3 l s a 11 + va 11 b 21 k 2 -vb 21 k 1 a 12 -vb 11 k 3 a 22 + b 11 k 3 l s a 21 +vb 11 k 1 a 22 -va 21 a 12 + va 11 a 22 -va 21 b 11 k 2 )/(k 5 v(-b 21 a 11 + b 11 a 21 ))v
From (12) it can be observed that in steady-state the yaw rate is equal to vρ and the lateral offset is equal to zero. Therefore, by choosing adequately the state feedback gains, it is possible to stabilise (11) and steer the vehicle to a zero lateral offset condition even in the presence of disturbances. Thus, the control strategy using the internal model fulfils the control objective of converging to the centre of the lane while negotiating curves.

The subsequent part presents the synthesis of the dynamical state feedback gain in order to ensure that the other control objectives are achieved.

Controller synthesis

Since it is intended to apply a control law that enables the vehicle to perform curves, ρ ̸ = 0, the equilibrium point ( 12) is no longer fixed at the origin (as in the case of straight lines). Therefore, to make use of the Lyapunov theory of stability, it is necessary to ensure that the state trajectories converge to an invariant set containing all the equilibrium points.

As presented by Minoiu-Enache et al. (2009a), it is possible, by means of BMI optimisation procedure, to find simultaneously a state feedback gain and a quadratic Lyapunov function, which guarantee the invariant set properties even in the presence of disturbances.

BMI contraints for Lyapunov invariant set

An ellipsoidal invariant set of the form:

ε(P ) ∆ = {x ∈ R n : x T P x ≤ 1}, with P = P T and P ≻ 0, (13)
with n = 6 for model (10), can be obtained from a quadratic Lyapunov function:

V (x) = x T P x, with P ≻ 0 and V (x) < 0. ( 14
)
The invariant set must contain all trajectories starting from a subset which represents the states at the activation of the assistance system. This subset of initial conditions is a polytope L(T N ) defined as:

x ∈ L(T N ) ∆ = {x ∈ R 6 : |T N x| ≤ 1}, (15)
where T N ∈ R 6x6 is a diagonal matrix with the following elements:

T N ∆ = diag[(β N ) -1 , (r N ) -1 , (ψ N L ) -1 , (y N L ) -1 , (α N 0 ) -1 , (α N 1 ) -1 ]. ( 16
)
where the values β N and r N can be set as the equilibrium point of (11) for ρ = ρ max and the remaining parameters are chosen based on safety conditions such as maximum vehicle positioning on the road. Minimising the ellipsoid that contains the polytope (15) at the same time as satisfying the conditions V < 0, V (x) ≥ 1, w T w ≤ 1 can be expressed in terms of a bilinear matrix inequality optimisation problem, by using the S-procedure [START_REF] Boyd | Linear Matrix Inequalties in System and Control Theory[END_REF] as follows:

Minimise:

tr(Q), det(Q) or λ max (Q) Subjected to: η ≻ 0, Q ≻ 0, QA T + Y T B T u + AQ + B u Y + ηQ B w B T w -η ≼ 0, 1 q T i q i Q ≽ 0, i = 1, ..., g/2, ( 17
)
where η is a real constant, Q = P -1 , Y = KP -1 , and q i represents each of the g = 64 vertices of the polytope (15). The ellipsoidal invariant set obtained form (17), delimits the bounds of the state trajectories despite the disturbance satisfying the control strategy.

Another control objective requires bounds on the control input, which can be achieved by adding a constraint to the optimisation problem. This constraint is important not only to limit the control input with respect to the actuator physical limitations but also to keep the vehicle dynamics within the linear behaviour, avoiding excessive lateral accelerations and saturation of the lateral tire forces. The constraint can be expressed as:

⎛ ⎝ 1 1 δ max f Y 1 δ max f Y T Q ⎞ ⎠ ≽ 0, (18)
where δ max f is the maximum steering input used in the control design phase. Whereas the contribution of the additional dynamics from the internal model improve the steady state performance of the lane keeping system, as shown in section 3.2, the transient performance can also be enhanced by including a constraint for pole clustering, as an extension of the control synthesis proposed by Minoiu-Enache et al. (2009a).

Linear matrix inequality constraint for pole clustering

In order to enhance the damping characteristic of the closed loop response, constraints can be included to specify regions where the poles can be located. A constraint for pole clustering defining a conic sector of angle θ centred at the origin [START_REF] Chiali | Hinf design with pole placement constraints: An lmi approach[END_REF], as shown in Figure 2, can be described as:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ sin θ AQ + QA T +B w Y + Y T B T w cos θ AQ -QA T +B w Y -Y T B T w cos θ -AQ + QA T -B w Y + Y T B T w sin θ AQ + QA T +B w Y + Y T B T w ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≺ 0 (19)
If the optimisation problem composed of ( 17), ( 18) and ( 19) is feasible, the resulting controller is able to reject the disturbances of road curvature and it ensures that all the state trajectories are contained in the invariant set ε(P ). Moreover, the pole clustering condition (19) ensures that all eigenvalues have negative real parts, therefore the linear system (11) is asymptotically stable even in the presence of disturbances. Thus, by taking this constraint into account, all requirements of the the controller, as listed in 3.1, are met.

Steering Assistance Activation Law

The human interactions play an important role in the driving assistance systems, therefore an adequate activation strategy is required. In order to avoid an intrusive assistance that may not be well accepted by drivers, it is important to consider the driver awareness, as well as the risk of lane departure or collision accidents.

The assessments of driver's awareness can be found for instance in the works of [START_REF] Tran | Towards a vision-based system exploring 3d driver posture dynamics for driver assistance: Issues and possibilities[END_REF] and [START_REF] Trivedi | Holistic sensing and active displays for intelligent driver support systems[END_REF], while [START_REF] Mccall | Video-based lane estimation and tracking for driver assistance: Survey, system, and evaluation[END_REF] present a survey of methods to estimate the risk of lane departure. Other sophisticated strategies can be used in the design of the driver assistance, but they are beyond the scope of this work, since it focuses the control synthesis and its performance.

Due to the simplicity of implementation, the strategy chosen for this work is similar to the one used in [START_REF] Minoiu-Enache | Driver steering assistance for lane departure avoidance[END_REF], which takes into account the driver's attentiveness as well as the risk of either lane departure during lane keeping manoeuvre.

The driver's attentiveness is estimated by measuring the driver's input torque on the steering wheel, denominated τ . If τ < 5N m the driver is considered inattentive. The driver has priority and control of the vehicle whenever his awareness is recovered, which is corresponds to τ > 2N m.

The risk of lane departure is estimated by the position of the front wheels with respect to a fixed strip of width 2d in the centre of the lane.

The positions of the left (right) front wheels, y l,(r) , can be described by the following:

y l = y L + (l f -l s )ψ L + a 2 , y r = y L + (l f -l s )ψ L -a 2 , ( 20
)
where a represents the vehicle width. Enforcing that the front wheels must remain inside the centred lane strip ±d, yields:

- 2d -a 2 ≤ y L + (l f -l s )ψ L ≤ 2d -a 2 (21)
This set, corresponding to the region between two parallel hyperplanes, can be written with respect to the state vector as:

T ∆ = {x ∈ R 6 : |F x| ≤ 1} ( 22
)
where F = (0, 0,

2(l f -ls)
2d-a , 2 2d-a , 0, 0). In order to ensure the invariant set property, the activation must take place inside the ellipsoid (13), therefore the risk of lane departure is considered when the vehicle leaves the set:

G ∆ = {ϵ(P ) ∩ T } (23)
The activation and deactivation of the assistance system follows:

• Activate if (τ < 5N m) & (|F x| ≥ 1) & (x ∈ ϵ(P )) • Deactivate if (τ ≥ 2N m),
where & is the Boolean operation AND.

Control Law Computation

The software PENBMI was used to solve the BMI problem ( 17), ( 18) and ( 19). The maximum road curvature was considered ρ max = 0.005m -1 . The polytopic region T N was defined with the following values:

T N = [0.013, 0.174, 0.017, 0.2, 0.005, 0.005] .

The bound for control input was defined as |δ max f | = 5°and the angle for pole clustering was chosen θ = 30°. The computed gain from the optimisation problem ( 17) is ensured to stabilise asymptotically the system for bounded disturbances |ρ| ≤ ρ max , since constraint ( 19) is satisfied and consequently all real parts of the eigenvalues of (11) are negative. For v = 15m/s the computed eigenvalues are:

λ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -6.7218 + 1.3347i -6.7218 -1.3347i -2.1680 -1.5181 -0.4520 -0.2470 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (25)
and the corresponding gain is given by: K = [-0.1813, -0.0955, -0.9418, -0.0781, -0.0045, -0.0341]

The implementation and performance analyses of the obtained control law in Carsim environment and in a prototype vehicle are discussed in the sequel.

Results

In order to verify the performance of the obtained state feedback gain (26), simulations have been carried out in CarSim environment taking into account a disturbance rejection situation. The implementation on a prototype vehicle and results obtained on the test track are discussed in section 4.2. The controller is also tested in an obstacle avoidance situation as shown in section 4.3.

Simulation Results on a CarSim Vehicle

CarSim vehicle model MacAdam (1981) uses detailed nonlinear tire models according to combined slip theory and takes into account the major kinematics and compliance effects of the suspensions (nonlinear spring models) and steering systems which were neglected at the control design stage. The standard CarSim F-Class vehicle has a nonlinear second order speed depending rack and pinion ratio steering system; a realistic actuator with a bandwidth of 10Hz is considered for the active steering.

In order to verify the control objective of zero steady state error while negotiating curves, the first simulation consists of a vehicle driving at 22m/s on a road with the geometrical trajectory as show in the first plot of Figure 3. The driver is considered to be inattentive throughout the simulation.

The other two plots of Figure 3 depict respectively the vehicle sideslip angle and yaw rate during the lane keeping manoeuvre.

The top plot of Figure 4 shows the lateral offset which clearly converges to zero despite of the increasingly road curvature, whereas the limited control input can be seen in the bottom plot.

Since the simulations results satisfactorily confirm the achievement of the control objectives, the same controller has been implemented in a prototype vehicle and the performance of the lane departure avoidance system has been evaluated in a test track as shown in the sequel.

Practical Implementation

The experimental vehicle used in the practical implementation (Figure 5) is equipped with a CORREVIT sensor to measure the vehicle sideslip angle β; an Inertial Navigation System to measure the vehicle yaw rate r and an odometer used to compute the vehicle longitudinal speed. A DC-motor mounted on the steering column provides the assistance torque through a PID controller which is used to track the control input in terms of steering angle δ f . An optical encoder is used to measure the steering angle, and load cell sensors integrated on the steering wheel are used to measure the driver torque.

A video camera is used to detect the lane markers by the vision algorithms of [START_REF] Labayrade | A reliable and robust lane detection system based on the parallel use of three algorithms for driving safety assistance[END_REF] permitting the measurements of lateral offset at look-ahead distance y L and relative yaw angle ψ L . The whole control algorithm, including the image capturing and processing, runs with a period of 0.04s.

The tests took place on a test track located in Satory-Versailles, 20km west of Paris, France. The track is 3.5km long consisting of typical road profiles as straight lines and bends of various radii. Figure 5 depicts an aerial view of the test track.

Two manoeuvres are presented, in both of them the vehicle longitudinal speed is 15 ± 2m/s and the driver is considered inattentive in the beginning of the manoeuvre, but retakes control in the end.

The first manoeuvre consists of on an S-shaped turn, with first turn to the left, as it can be noticed with the bottom subplot of Figure 8. The positioning of the vehicle on the lane is shown on Figure 6. The first subplot shows the relative yaw angle while in the second it is shown the position of the front wheels with respect to the lane, which is shown by the dashed-dot line. The activation occurs at 51s approximately when the right wheel crosses the fixed centred strip of ±d, depicted by the dashed line. The front wheels remain in the lane during the whole control action. The driver recovers attention at 63s, and the controller is deactivated, since he/she apply a steering torque greater then 2N m.

The dynamics of vehicle during the intervention of the lane keeping system can be seen in Figure 7. The vehicle sideslip angle, yaw rate and longitudinal speed are depicted. The required control input is shown in the top subplot of Figure 8. In the bottom subplot it can be seen the road curvature for the corresponding part of the track. Even though the values of the road curvature are larger than the bound chosen in the design phase (ρ max = 0.005 m -1 ), the controller performs satisfactory keeping the vehicle within the lane limits.

Another manoeuvre is carried out in order to show the lane keeping controller performance. As it can be seen by the road curvature of the test track depicted in Figure 11, this scenario corresponds to an inattentive driver on a curve to the right side, which is also sharper than the considered bound ρ max in the design stage.

The activation, shown in the bottom subplot of Figure 9, takes place due to the positioning of the right front wheel, which reaches the centred strip on the interior side of the curve. Even though the correction induces a longer excursion of the vehicle near the exterior side of the lane, the performance is satisfactory, as the front wheels do not exit the lane. The relative yaw angle is plotted on the top subplot of Figure 9.

The dynamics of vehicle during the control action are illustrated on Figure 10. The subplot shows the vehicle sideslip angle, yaw rate and longitudinal speed respectively. Finally the required control input is shown in the top subplot of Figure 11.

Due to the satisfactory performance of the lane keeping system, obstacle avoidance manoeuvres with the same controller are performed in the CarSim simulation environment on a vehicle model similar to the experimental vehicle. The results are shown in the next section.

Obstacle avoidance manoeuvre

As an attempt to test the controller in different emergency situations, simulations on Carsim environment are carried out considering a double lane change manoeuvre from ISO/DIS3888 standard. This test is well known for being severe, as the vehicle is submitted to lateral tyre forces which enter the non-linear behaviour (which has been neglected in the design phase). It is usually applied for the validation of advanced driver assistance systems, commercially known as ESC (electronic stability control), which may include active steering, independently wheel braking and active suspension.

The path considered in the test is shown in Figure 12. The test is defined by a predetermined placement of cones on the road and the manoeuvre is carried out on a dry surface. The vehicle longitudinal speed is 25m/s.

The resulting vehicle sideslip angle and the yaw rate are depicted in Figure 12. The lateral offset y l and the corresponding control action are shown in Figure 13. The controller achieves the path following with a reduced path following error in the lateral direction (0.4m as can be observed in the first subplot of Figure 13.

Conclusions

This work has proposed a lane keeping assistance system which takes into account simultaneously the attentiveness of the driver and the risk either of unintended lane departure. Moreover, it was shown that it is possible to reject the disturbances, steering the vehicle to zero lateral offset, by modelling the exogenous input (road curvature) into an internal model with two integration of the lateral offset signal. The compensation of the disturbances is guaranteed by the dynamic state feedback controller. The design of the controller also takes into account limitation in the control input.

The problem of finding a feedback gain that stabilises the system with bounded control input, can be casted as an BMI optimisation problem. It has been shown that pole-clustering constraint can been added to the BMI optimisation problem in order to improve the damping characteristic of the closed-loop response.

Simulations on CarSim environment showed that the proposed controller is able to ensure zero lateral offset in steady-state and satisfactorily transient response while rejecting a ramp disturbance corresponding to an increasingly road curvature.

The proposed control law has been implemented in a Peugeot 307 prototype vehicle and tested on the 3.5Km long test track located in Satory, 20Km west of Paris. The satisfactory performance of the practical tests confirm the use of the assistance as a lane departure avoidance system.

Moreover, simulations of a double lane-change manoeuvre have been performed with the CARSIM environment, showing satisfactorily performance of the controller when lateral tyre forces are beyond the linear behaviour in addition to its potential as a controller for obstacle avoidance. Such use is still to be evaluated on the test track, as well as testing the lane departure avoidance system on degraded road adherence.

Also as a perspective from this work, efforts should be concentrated in designing an output feedback controller in order to use only on-board vehicle sensors rendering the controller financially accessible for car manufactures.

List of Tables and Figures

Figure 1 Figure 2 Figure 3 Figure 4

 1234 Figure 1 Dynamics and positioning variables of a single track vehicle model

Figure 5

 5 Figure 5 Aerial View from Satory-Versailles test track. Photo from: TGLR, Atelier Gérard Leroux, 13, rue Pierre Clavillier, 78800 Houilles, France and experimental vehicle.

Figure 12 Figure 13

 1213 Figure 12 (a) Standard CarSim Moose test path, (b) vehicle sideslip angle β and (c) yaw rate r for manoeuvre performed at a longitudinal speed of v = 25m/s.

Table 1

 1 Values of the vehicle parameters.

Research and Safety, Versailles, France, where he was involved in research on traffic network control. From 1994 to 2002, he was an Assistant Professor with the Evry val d'Essonne University, French Embassy Hague, The Netherlands. has been supervisor of 30 PhD students and several of Master students. He has been leader and involved in many industry, government and EU-funded research projects. as a researcher, in


Trung-Dung Vu

Julien Burlet

Olivier Aycard

Grid-based Localization and Online Mapping with Moving Objects Detection and Tracking: new results

In this paper, we present a real-time algorithm for local simultaneous localization and mapping (SLAM) with detection and tracking of moving objects (DATMO) in dynamic outdoor environments from a moving vehicle equipped with a laser scanner. To correct vehicle location from odometry we introduce a new fast implementation of incremental scan matching method that can work reliably in dynamic outdoor environments. After a good vehicle location is estimated, the surrounding map is updated incrementally and moving objects are detected without a priori knowledge of the targets. Detected moving objects are finally tracked by a Multiple Hypothesis Tracker (MHT) coupled with an adaptive IMM (Interacting Multiple Models) Filter. The experimental results on datasets collected from different scenarios such as: urban streets, country roads and highways demonstrate the efficiency of the proposed algorithm on a Daimler Mercedes demonstrator in the framework of the European Project PReVENT-ProFusion2.

I. INTRODUCTION

Perceiving or understanding the environment surrounding of a vehicle is a very important step in driving assistant systems or autonomous vehicles. The task involves both simultaneous localization and mapping (SLAM) and detection and tracking of moving objects (DATMO). While SLAM provides the vehicle with a map of static parts of the environment as well as its location in the map, DATMO allows the vehicle being aware of dynamic entities around, tracking them and predicting their future behaviors. It is believed that if we are able to accomplish both SLAM and DATMO in real time, we can detect every critical situations to warn the driver in advance and this will certainly improve driving safety and can prevent traffic accidents.

Recently, there have been considerable research efforts focusing on these problems [START_REF] Hähnel | Mobile robot mapping in populated environments[END_REF] [START_REF] Wang | Simultaneous Localization, Mapping and Moving Object Tracking[END_REF]. However, for highly dynamic outdoor environments like crowded urban streets, there still remains many open questions. These include, how to represent the vehicle environment, how to obtain a precise location of the vehicle in presence of dynamic entities, and how to differentiate moving objects and stationary objects as well as how to track moving objects over time.

In this context, we design and develop a generic architecture to solve SLAM and DATMO in dynamic outdoor environments. This architecture (Fig. 2) is divided into two main parts: the first part where the vehicle environment is mapped, fusion between different sensors is performed and moving objects are detected; and the second part where previously detected moving objects are verified and tracked. This architecture is currently used in the framework of the European project PReVENT-ProFusion 1 . The goal of this project is to design and develop generic architectures to perform perception tasks (ie, mapping of the environment, localization of the vehicle in the map, and detection and tracking of moving objects). In this context, our architecture has been integrated and tested on two demonstrators: a Daimler-Mercedes demonstrator and a Volvo truck demonstrator. In previous paper [START_REF] Td | Online localization and mapping with moving objects tracking in dynamic outdoor environments[END_REF], a description of the first level is reported. In this paper, we detail the description of the sensor data fusion part of the first level and the second level and show some results on the Daimler-Mercedes demonstrator moving at high speed.

The rest of the paper is organized as follows. In the next section, we present the Daimler Mercedes demonstrator. A brief overview of our architecture is given in section III. Description of first level of architecture is summarized and Sensor Data Fusion is described in Section IV. Second level is detailed in section V. Experimental results are given in Section VI and finally in Section VII conclusions and future works are discussed.

II. THE DAIMLER MERCEDES DEMONSTRATOR

The DaimlerChrysler demonstrator car is equipped with a camera, two short range radar sensors and a laser scanner (Fig. 1). The radar sensor is with a maximum range of 30m and a field of view of 80 • . The maximum range of laser sensor is 80m with a field of view of 160 • and a horizontal Fig. 2. Architecture of the perception system resolution of 1 • . In addition, vehicle odometry information such as velocity and yaw rate are provided by the vehicle sensors. The measurement cycle of the sensor system is 40ms. Images from camera are for visualization purpose.

III. GENERAL ARCHITECTURE

We design and develop a generic architecture (Fig. 2) to solve SLAM and DATMO in dynamic outdoor environments.

In the first part of the architecture, to model the environment surrounding the vehicle, we use the Occupancy Grid framework developed by Elfes [START_REF] Elfes | Occupancy grids: a probabilistic framework for robot percpetion and navigation[END_REF]. Compared with featurebased approaches, grid maps can represent any environment and are specially suitable for noisy sensors in outdoor environments where features are hard to define and extract.

In general, in order to perform mapping or modelling the environment from a moving vehicle, a precise vehicle localization is essential. To correct vehicle locations from odometry, we introduce a new fast laser-based incremental localization method that can work reliably in dynamic environments. When good vehicle locations are estimated, by integrating laser measurements we are able to build a consistent grid map surrounding of the vehicle. Finally by comparing new measurements with the previously constructed local vehicle map, dynamic objects then can be detected.

Finally, sensor data coming from different sensors are fused.

In the second part, detected moving objects in the vehicle environment are tracked. Since some objects may be occluded or some are false alarms, multi objects tracking helps to identify occluded objects, recognize false alarms and reduce mis-detections.

IV. FIRST LEVEL

In this section, we first summarized the description of the first level of our architecture: Environment Mapping & Localization, Moving Objects Detection. More details on the two first parts could be found in [START_REF] Td | Online localization and mapping with moving objects tracking in dynamic outdoor environments[END_REF]. In the last subsection, we describe the fusion between objects detected by laser and radar data.

A. Environment Mapping & Localization

To map the environment and localize in the environment, we propose an incremental mapping approach based on a fast laser scan matching algorithm in order to build a consistent local vehicle map. The map is updated incrementally when new data measurements arrive along with good estimates of vehicle locations obtained from the scan matching algorithm. The advantages of our incremental approach are that the computation can be carried out very quickly and the whole process is able to run online.

1) Environment mapping using Occupancy Grid Map: Using occupancy grid representation, the vehicle environment is divided into a two-dimensional lattice M of rectangular cells and each cell is associated with a measure taking a real value in [0, 1] indicating the probability that the cell is occupied by an obstacle. A high value of occupancy grid indicates the cell is occupied and a low value means the cell is free. Suppose that occupancy states of individual grid cells are independent, the objective of a mapping algorithm is to estimate the posterior probability of occupancy P(m | x 1:t , z 1:t ) for each cell of grid m, given observations z 1:t = {z 1 , ..., z t } from time 1 to time t at corresponding known poses x 1:t = {x 1 , ..., x t }.from time 1 to time t.

2) Localization of the vehicle in the Occupancy Grid Map: In order to build a consistent map of the environment, a good vehicle localization is required. Because of the inherent error, using only odometry often results in an unsatisfying map.

To solve this problem, we used a particle filter. We predict different possible positions of the vehicle (ie, one position of the vehicle corresponds to one particle) using the motion model and compute the probability of each position (ie, the probability of each particle) using the laser data and a sensor model.

B. Moving Objects Detection

After a consistent local grid map of the vehicle is constructed, moving objects can be detected when new laser measurements arrive by comparing with the previously constructed grid map. The principal idea is based on the inconsistencies between observed free space and occupied space in the local map. If an object is detected on a location previously seen as free space, then it is a moving object. If an object is observed on a location previously occupied then it probably is static. If an object appears in a previously not observed location, then it can be static or dynamic and we set the unknown status for the object in this case.

C. Fusion with radars

After moving objects are identified from laser data, we confirm the object detection results by fusing with radar data and provide the detected objects with their velocities. For each moving object detected from laser data as described in the previous section, a rectangular bounding box is calculated and the radar measurements which lie within the box region are then assigned to corresponding object. The velocity of the detected moving object is estimated as the average of these corresponding radar measurements. Figure 4 shows an example of how the fusion process takes place. Moving objects detected by the Laserscanner are displayed in red with green bounding boxes. The targets detected by two radar sensors are represented as small circles in different colors along with corresponding velocities. We can see in the radar field of view, that two objects detected by the Laserscanner are also seen by two radars so that they are confirmed and their velocities are estimated. Radar measurements that do not correspond to any dynamic object or fall into another region of the grid are not considered.

V. SECOND LEVEL

In general, the multi objects tracking problem is complex: it includes the definition of tracking methods, but also association methods and maintenance of the list of objects currently present in the environment [START_REF] Bar-Shalom | Tracking and Data Association[END_REF] [START_REF] Thrun | Probabilistic Robotics (Intelligent Robotics and Autonomous Agents)[END_REF]. Regarding tracking techniques, Bayesian filters [START_REF] Arulampalam | A tutorial on particle filter for online nonlinear/non-gaussian bayesian tracking[END_REF] are generally used. These filters require the definition of a specific motion model of tracked objects to predict their positions in the environment. Using this prediction and some observations, in a second stage, an estimation of the position of each object present in the environment is computed.

In this section, we describe the four different parts of our architecture (figure 2) to solve the different parts of multiobjects tracking:

• The first one is the gating. In this part, taking as input predictions from previous computed tracks, we compute the set of new detected objects which can be associated to each track. • In a second part, using the result of the gating, we perform objects to tracks association and generate association hypothesis, each track corresponding to a previously known moving object. Output is composed of the computed set of association hypothesis. • In the third part called tracks management, tracks are confirmed, deleted or created according to the association results and a pruned set of association hypothesis is output. • In the last part corresponding to the filtering step, estimates are computed for 'surviving' tracks and predictions are performed to be used the next step of the algorithm. In this part, we use an adaptive method based on Interacting Multiple Models (IMM). More details about these different parts are outlined next.

A. Gating

In this part, taking as input predictions from previous computed tracks and newly detected objects, a gating is performed. It consists in, according to an arbitrary distance function, determine the detected objects which can be associated with tracks. Also during this stage, clustering is performed in order to reduce the number of association hypothesis. It consists in making clusters of tracks which share at least one detected object. In the next stage, association can be performed independently for each cluster decomposing a large problem in smaller problems which induce generation of less hypothesis.

B. Association

In this part, taking as input clusters of tracks and detected objects validated by the gating stage, association hypothesis are evaluated. By considering likelihood of objects with tracks, new track apparition probability and non-detection probability, an association matrix is formed. Let be L(o i ,t j ) the function giving the likelihood of object i with track j, P NT the new track apparition probability and P ND the non detection probability. Taking as an example the situation in the Figure 5, the association matrix is written:

    L(o 1 ,t 1 ) -∞ P NT L(o 2 ,t 1 ) L(o 2 ,t 2 ) P NT -∞ L(o 3 ,t 2 ) P NT P ND P ND -∞    
Thus a possible association hypothesis corresponds to a valid assignation in the matrix of detected objects with tracks i.e one unique element in each row and each column is chosen to compose the assignation. In order to reduce the number of hypothesis, only the m-best association hypothesis are considered. The m-best assignment in the association matrix are computed using the Murty method [START_REF] Murty | An algorithm for ranking all the assignments in order of increasing costs[END_REF] which computes the m-Best assignations in the matrix and by this way we obtain the m-Best Hypothesis.

C. Track management

In this third stage, using the m-Best Hypothesis resulting of the association stage, the set of tracks, is maintained i.e tracks are confirmed, deleted or created.

New tracks are created if a new track creation hypothesis appears in the m-best hypothesis. A new created track is confirmed if it is updated by detected objects after a fixed number of algorithm steps (three in our implementation). Thus spurious measurements which can be detected as objects in the first step of our method are never confirmed.

If a non-detection hypothesis appears and so to deal with non-detection cases (which can appear for instance when an object is occulted by an other one, tracks without associated detected objects are updated according to their last associated objects and next filtering stage becomes a simple prediction. But if a track is not updated by a detected object for a given number of steps, it is deleted.

D. Adaptive Filtering using Interacting Multiple Models

In this filtering stage, according to previously computed predictions, estimations are performed for each association of all hypothesis and new predictions are computed for the gating stage. Regarding filtering techniques, there exists several kinds of filters, the most classical is the well known Kalman filter. But in all kinds of filters, the motion model is the main part of the prediction step.

To deal with these motion uncertainties, Interacting Multiple Models (IMM) [START_REF] Li | A survey of maneuvering target tracking-part v: Multiple-model methods[END_REF] have been successfully applied in several applications [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF]. The IMM approach overcomes the difficulty due to motion uncertainty by using more than one motion model. The principle is to assume a set of motion models as possible candidates of the true displacement model of the object at one time. To do so, a bank of elemental filters is ran at each time, each corresponding to a specific motion model, and the final state estimation is obtained by merging the results of all elemental filters according to the distribution probability over the set of motion models (in the next part we note µ this probability). By this way different motion models are taken into account during filtering process. As the quality of gating relies directly on the quality of filtering and especially the prediction step, we have chosen Interacting Multiple Models (IMM) [START_REF] Li | A survey of maneuvering target tracking-part v: Multiple-model methods[END_REF] to deal with motion uncertainties in this filtering part.

Besides, we developed an efficient method in which critical parameter of the IMM is on-line adapted according to the most probable trajectories formed by tracks. Thus as Fig. 6 shows our filtering stage is composed of three parts : an IMM filtering part, a part in which most probable trajectories are computed and a last part in which we adapt the IMM filter. These three parts are described in the next paragraphs. Fig. 7. The sixteen chosen motion models in the vehicle's frame 1) Definition of our IMM: Nevertheless, to apply IMM on real applications a number of critical parameters have to be defined, for instance the set of motion models and the transition probability matrix(TPM). To cope with this design step which can no match the reality, we propose an efficient method in which the TPM is on-line adapted.

In our specific application, different objects such as cars or motorcycles can move in any directions and can often change their motions. Thus in our aim we choose various IMM's motion models to cover the set of possible directions and velocities. As each filter corresponds to a specific motion model, we have to define each motion model. So, assuming we have different possible velocities defined according to the vehicle velocity and eight directions in the set of possible directions an object can follow, we obtain sixteen motion models (Fig. 7).

Hence, according to the definition of these sixteen motion models, our IMM is composed of sixteen kalman filters.

The TPM is initially chosen to be uniform.

2) Computation of the most probable trajectories: Once estimates are performed for all hypothesis, the most probable trajectory is computed for each track. This step permits to give users more readability on what is happening during tracking process and also permits us to adapt on-line the IMM parameter according to these trajectories.

3) Adaptation of the IMM: To adapt the TPM in our specific situation i.e tracking detected objects, most probable trajectories are considered. Taking as input the set of trajectories computed during filtering process, we will adapt online the TPM of the IMM filter in order to obtain a better transition between motion models close to the real behavior of tracked objects.

The principle is the following. For a given number N of trajectories we build sequences of associated motion models probabilities. And then, using these motion models probabilities, the TPM is adapted and reused in the IMM filters for the next estimations. In more details, algorithm 1, given in pseudo-code, is the algorithm defined to compute one adaptation of the TPM. Return TPM in IMM An adaptation of the TPM is done after a given number N of trajectories obtained from tracks, to update TPM using a window on trajectories (cf. loop line 3-19 of algorithm 1). Moreover trajectories are processed one by one in three steps:

1-Models' probabilities are collected by travel through the computed most probable sequence 2-Most probable models' sequence is computed 3-Most probable models' transitions are quantified Collection of models' probabilities : : For each part of a given most probable trajectory computed in last stages of the filtering process, we collect the distribution over models(lines 7). Thus a model probabilities' sequence S n obtained in such a way and is stored to be processed (line 8).

Computation of the most probable model sequence : : In a next step, the most probable models' sequence of S n is computed (line 11). More precisely, considering the actual TPM and a set S n = µ 0 ...µ K of model probabilities through time 0 to K, we aim to obtain the most probable models' sequence knowing the estimates computed by the IMM:

Max P(µ 0 µ 1 ...µ k | x 0 x 1 ... x K ) (1)
We just need to obtain the maximum of the distribution P(µ 1 µ 2 ...µ K | x 0 x 1 ... x K ), thus the inference is made using the Viterbi Data Algorithm [START_REF] Forney | The viterbi algorithm[END_REF]. As complexity of this algorithm is in O(KM 2 ), we efficiently obtain the most probable models' sequence.

Quantification of most probable model transitions : : Using this most probable models' sequence, the number of transitions from one model to an other is quantified (lines 13 to 17). To do so a frequencies matrix is considered. This matrix models the number of transitions which have occurred from one model to an other. We note F this matrix and so F i j gives the number of transitions which has occurred from model i to j. Using the most probable models' sequence corresponding to a specific trajectory and computed by the Viterbi algorithm, the update of F is directly obtained by counting transitions in this sequence. Furthermore, F is kept in memory to be used in next adaptation and before the first update all its elements are set to 1.

Finally, when N trajectories have been treated, the new TPM is obtained by normalization of the frequencies matrix F. Thus the TPM is re-estimated using all model sequences S 1 ...S N and is reused in the IMM for next executions (lines 21 and 22). In practice, before the first run, the TPM is chosen uniform (according to F initialization) as we do not want to introduce a priori data.

VI. EXPERIMENTAL RESULTS

The detection and tracking results are shown in Fig. 8. The images in the first row represent online maps and objects moving in the vicinity of the vehicle are detected and tracked. The current vehicle location is represented by blue box along with its trajectories after correction from the odometry. The red points are current laser measurements that are identified as belonging to dynamic objects. Green boxes indicate detected and tracked moving objects with corresponding tracks displayed in different colors. Information on velocities is displayed next to detected objects if available. The second row are images for visual references to corresponding situations.

In Fig. 8, the leftmost column depicts a scenario where the demonstrator car is moving at a very high speed of about 100 kph while a car moving in the same direction in front of it is detected and tracked. On the rightmost is a situation where the demonstrator car is moving at 50 kph on a country road. A car moving ahead and two other cars in the opposite direction are all recognized. Note that the two cars on the Fig. 8. Experimental results show that our algorithm can successfully perform both SLAM and DATMO in real time for different environments left lane are only observed during a very short period of time but both are detected and tracked successfully. The third situation in the middle, the demonstrator is moving quite slowly at about 20 kph in a crowded city street. Our system is able to detect and track both the other vehicles and the motorbike surrounding. In all three cases, precise trajectories of the demonstrator are achieved and local maps around the vehicle are constructed consistently. In our implementation, the computational time required to perform both SLAM and DATMO for each scan is about 20 -30 ms on a 1.86GHz, 1Gb RAM laptop running Linux. This confirms that our algorithm is absolutely able to run synchronously data cycle in real time. More results and videos can be found at http: //emotion.inrialpes.fr/ ∼ tdvu/videos/.

VII. CONCLUSIONS AND FUTURE WORKS

We have presented an approach to accomplish online mapping and moving object tracking simultaneously. Experimental results have shown that our system can successfully perform a real time mapping and moving object tracking from a vehicle at high speeds in different dynamic outdoor scenarios. This is done based on a fast scan matching algorithm that allows estimating precise vehicle locations and building a consistent map surrounding of the vehicle. After a consistent local vehicle map is built, moving objects are detected and are tracked using an adaptive Interacting Multiple Models filter coupled with an Multiple Hypothesis tracker.

Fig. 1 .

 1 Fig. 1. The Daimler Mercedes demonstrator car.

Fig. 3 .

 3 Fig. 3. Moving object detection example. See text for more details.

Fig. 4 .

 4 Fig. 4. Moving object detected from laser data is confirmed by radar data.

Fig. 3

 3 Fig. 3 illustrates the described steps in detecting moving objects. The leftmost image depicts the situation where the vehicle is moving along a street seeing a car moving ahead and a motorbike moving in the opposite direction. The middle image shows the local static map and the vehicle location with the current laser scan drawn in red. Measurements which fall into free region in the static map are detected as dynamic and are displayed in the rightmost image. After the clustering step, two moving objects are identified (in green boxes) and correctly corresponds to the car and the motorbike.

Fig. 5 .

 5 Fig. 5. Example of association problem

Fig. 6 .

 6 Fig. 6. Principle of our adaptive filtering program

www.prevent-ip.org/profusion

VIII. ACKNOWLEDGMENTS

The work is supported by the European project PReVENT-ProFusion, and partially by the Délégation Générale pour L'Armement (DGA).


G Chmaycem
email: ghada.chmaycem@gmail.com

M Jazar
email: mjazar@laser-lb.org

R Monneau
email: monneau@cermics.enpc.fr

A new contraction family for porous medium and fast diffusion equations

Keywords: AMS Classification: 35K55, 35K65 contraction, porous medium equation, fast diffusion

In this paper, we present a surprising two-dimensional contraction family for porous medium and fast diffusion equations. This approach provides new a priori estimates on the solutions, even for the standard heat equation.

Introduction

In this paper, we answer a long standing open question about the existence of new contractions for porous medium type equations. For m > 0 and d ≥ 1, we consider nonnegative solutions U (t, x) of the following normalized equation (1.1) mU t = ∆U m on (0, +∞) × Ω = Q with Ω = R d .

For m > 1, this equation is called the porous medium equation (PME), while for m < 1, it is called the fast diffusion equation (FDE). We choose the following standard definition of solution:

Definition 1.1 (Notion of solution) Let m > 0. We say that U ≥ 0 is a solution of (1.1) with intial data 0 ≤ U 0 ∈ L 1 (Ω), if U, U m ∈ L 1 loc (Q) and for every function η ∈ C ∞ c (Q) with compact support, we have

Q U m m ∆η + U η t dxdt + Ω U 0 (x)η(0, x) dx = 0.
We recall that for m > 0 and any initial data satisfying

(1.2) U 0 ∈ L 1 (Ω) ∩ L ∞ (Ω) with U 0 ≥ 0,
there is exactly one solution U ∈ C([0, +∞); L 1 (Ω)) ∩ L ∞ (Q) in the sense of Definition 1.1 (see [START_REF] Vázquez | Smoothing and Decay Estimates for Nonlinear Diffusion Equations: Equations of Porous Medium Type[END_REF][START_REF] Vázquez | The Porous Medium Equation. Mathematical Theory, Oxford mathematical monographs[END_REF]). In the special case where m > m c (d) := max 0, d-2 d , it is well known that it is possible to remove the condition U 0 ∈ L ∞ (Ω) to get existence and uniqueness in some weaker space. To present our results in a concise and unified way for m > 0, we will stick under assumptions (1.2).

In any dimensions, it is known (see [START_REF] Vázquez | The Porous Medium Equation. New contractivity results[END_REF] and the references therein), that this equation has at least three types of contractions: in L 1 , in H -1 and finally for the 2-Wasserstein distance if m > m c (d). This last contraction property has been discovered by Otto in [START_REF] Otto | Double degenerate diffusion equations as steepest descent[END_REF] and published in [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF] (see also later [START_REF] Carrillo | Contractions in the 2-Wasserstein Length Space and Thermalization of Granular Media[END_REF]). It is also known for the PME, that for the p-Wasserstein distance, there is no contraction for d ≥ 2 and p > p 1 (m, d), while there is contraction for any p ∈ [1, +∞] in dimension d = 1 (see [START_REF] Vázquez | The Porous Medium Equation. New contractivity results[END_REF]). Note also that in [START_REF] Vázquez | The Porous Medium Equation. Mathematical Theory, Oxford mathematical monographs[END_REF] (Theorem A.5 p. 583), the author proves that the PME is not contractive with respect to L p norm for m ≥ 2 if p > p 2 (m, d).

In this paper, for 0 < m < 2, we present a new family of contractions for this equation in any dimensions, which extends the L 1 contraction properties. Our contraction can be seen as the fourth known contraction for this equation. Even for the case m = 1, our approach leads to new results for the standard heat equation. More precisely, for U and V two nonnegative solutions of (1.1), we show that the following quantity

(1.3) Ω |U α -V α | p
is a Lyapunov functional which is nonincreasing in time for all (α, p) in some admissible set. For convenience, we will work in the whole paper with

n = m -1.
For 0 < |n| < 1, we define the admissible (convex) set (which is skeched on Figure 1)

(1.4) K |n| = (α, p) ∈ [|n|, 1] × R; P -(α) ≤ p ≤ P + (α), with P ± (α) = 1 + 2 n 2 (1 -α)(α ± √ α 2 -n 2 ) .
For n = 0, we also set K 0 = {(α, p) ∈ (0, 1] × R; αp ≥ 1}. One may wonder if contraction (1.3) is related or not to some gradient flow structure of the equation. Indeed, for a given positive solution U of (1.1), we can set u = U α which solves (1.5)

u t = u γ-1 u∆u + γ|∇u| 2 with γ = n α
and γ = γ -1 + α -1 . In the special case where γ = γ/2, it can be seen that equation (1.5) is the negative L 2 -gradient flow of some energy, i.e. it solves

u t = -∇ L 2 E(u) with E(u) = Ω u γ |∇u| 2 2 .
Moreover, it is easy to check (computing the hessian of E) that E is convex if γ ∈ [-1, 0]. This corresponds exactly to the points (α, p) = (1 + n 2 , 2) ∈ K |n| that we capture for n ∈ [- 2 3 , 0]. This interpretation is similar to the derivation of Yamabe flow (see for instance [START_REF] Ye | Global existence and convergence of Yamabe flow[END_REF]) with the difference here that the exponent of the fast diffusion is not directly related to the space dimension. Nevertheless, except this very exceptional case, contraction (1.3) does not seem to be related to any gradient flow structure. The paper is organized as follows. In Section 2, we present our main results. In Section 3, we prove the new contraction property. Finally in Section 4, we prove several gradient decay estimates.

Main results

We state precisely our main results. With the previous notation, and for u, v > 0, we define the symmetric matrix for p > 1 and α ∈ (0, 1]: (2.6)

Q α,p (v, u) = (p -1)    v γ 1 + Γ u v -1 - 1 2 (v γ + u γ ) - 1 2 (v γ + u γ ) u γ 1 + Γ v u -1    with Γ = 1 -α α(p -1)
.

and set Q 1,1 (v, u) = 0.
We can check that this matrix is nonnegative for all u, v > 0 if and only if (α, p) ∈ K |n| (see Lemma 3.1). For y ∈ R, we set y + = max(0, y). Then our main result is the following

≤ U 0 ∈ L 1 (Ω) ∩ L ∞ (Ω) (resp. 0 ≤ V 0 ∈ L 1 (Ω) ∩ L ∞ (Ω)). Let (α, p) ∈ K |n| where K |n| is defined in (1.4). We set u = U α and v = V α .
Then for w = u and v, we have

w ∈ C(Q) ∩ L ∞ (Q) ∩ C([0, +∞); L p (Ω)) and w ∈ C ∞ (Q ∩ {w > 0}). Moreover, with the notation v(τ ) = v(τ, •), u(τ ) = u(τ, •), we have for all t ≥ 0 (2.7) Ω (v(t) -u(t)) p + dx + p t 0 dτ {v-u>0}∩{u>0} e[v, u] dx ≤ Ω (v(0) -u(0)) p + dx with (2.8) e[v, u] = |v -u| p-2 ∇v ∇u T Q α,p (v, u) ∇v ∇u ≥ 0, for v > u > 0.
Note that it is possible to see that there is equality in (2.7), for instance if u and v are smooth and positive on the torus T d instead of the whole space. We now define a second symmetric matrix M 1,1 = 0 and for (α, p)

∈ K |n| \ {(1, 1)} (2.9) M α,p = (p -1)     1 -Γ + γ 2 -Γ + γ 2 Γ(1 -γ)     ,
which turns out to be also nonnegative (see Lemma 4.1).

Then, we have also the following result Theorem 2.2 (Decay of the gradient) Let U be the solution of equation (1.1) with initial data U 0 , under the assumptions of Theorem 2.1. Assume that ∇U α 0 ∈ L p (Ω). Setting u = U α , we have ∇u ∈ L ∞ ([0, +∞); L p (Ω)). Moreover, with the notation u(τ ) = u(τ, •), we have for almost every t ≥ 0:

(2.10)

Ω |∇u(t)| p dx + p t 0 dτ {∇u =0}∩{u>0} ē[u, ∇u, D 2 u]dx ≤ Ω |∇u(0)| p dx. with (2.11) ē[u, w, A] = |w| p-2 u γ-2 uA 0 w 2 T M α,p uA 0 w 2 + |w| p-2 u γ |A 2 | 2 + p|A 1 | 2 ≥ 0,
where we write a general symmetric matrix A as

A = A 0 A 1 A T 1 A 2 by blocks in the vector space Rb ⊕ b ⊥ with b = w |w| .
Again, note that it is possible to see that there is equality in (2.10), if U 0 is smooth and positive with Ω = T d .

Let us mention that our results will be used in a future work to get information on the anomalous exponents associated to self-similar solutions like the one of Aronson-Graveleau (see [START_REF] Aronson | A Self-similar solution to the focusing problem for the porous medium equation[END_REF]), and also the one of King-Peletier-Zhang (see [START_REF] King | Self-similar behaviour for the equation of fast nonlinear diffusion[END_REF][START_REF] Peletier | Self-similar solutions of a fast diffusion equation that do not conserve mass[END_REF]).

The new contraction property

This section is devoted to the proof of Theorem 2.1. To this end, we first need the following result whose proof is elementary.

Lemma 3.1 (Properties of the matrix Q α,p ) For |n| < 1, let (α, p) ∈ Int K |n| where K |n| is defined in (1.4). i) Positive definite matrix Let η > 0. Then for v > u > 0 with 1 + η -1 > v/u ≥ 1 + η, we have (3.12) Q α,p (v, u) ≥ v γ νI 2 ,
where I 2 is the 2 × 2 identity matrix, Q α,p is defined in (2.6), and ν = ν(n, α, p, η) > 0.

ii) Bound from below for the top diagonal term Let v ≥ u > 0, then the top diagonal term of the matrix satisfies

(3.13) (Q α,p (v, u)) 11 ≥ v γ ν 0 ,
for some ν 0 = ν 0 (n, α, p) > 0.

Proof of Lemma 3.1

We set w = u/v ∈ (0, 1] and start with α ∈ (0, 1) and p > 1.

Step 1: on the determinant A direct computation shows that det Q α,p (v, u) > 0 if and only if

(3.14) G γ (w) := Γ (1 -Γ) f 1 (w) - 1 4 f γ (w) > 0 with f γ (w) := w γ -2 + 1 w γ .
Note that G γ (1) = G γ (1) = 0. We now show that G γ is strictly convex for w ≤ 1, which will imply (3.14) for w = 1. Indeed, we have

G γ (w) = -w -3 F γ (w) with F γ (w) := 1 4 γ (γ -1)w γ+1 + (γ + 1)w 1-γ -2Γ (1 -Γ) ≤ F γ (1),
where we have used |γ| < 1 and w ≤ 1 to get the inequality. We see that F γ (1) < 0 if and only if Γ is in between the two distinct roots Γ

± := 1 ± 1 -γ 2 2
, which is equivalent to (α, p) ∈ Int K |n| . This implies the strict convexity of G γ .

Step 2: conclusion

Recall that for (α, p) ∈ Int K |n| , we have p > P -(α) ≥ 1 α . We deduce that each diagonal term of Q α,p is positive (which shows ii)) and then tr Q α,p (v, u) > 0. With Step 1, this implies the result i).

Proof of Theorem 2.1

We only deal with the case (α, p) ∈ Int K |n| . Note that the border case can easily be recovered, by a passage to the limit. The proof is divided into several steps.

Step 1: Estimate for smooth positive solutions Given two nonnegative and smooth initial data U 0 , V 0 with compact support, we consider modified initial data for ε ∈ (0, 1)

U 0ε = U 0 + ε and V 0ε = V 0 + ε.
According to the standard parabolic theory (see [START_REF] Ladyženskya | Linear and Quasilinear Equations of Parabolic Type[END_REF]), we can consider the smooth positive functions U ε , V ε , which are the classical solutions of (1.1) respectively associated to initial data U 0ε and V 0ε . We have moreover the obvious bounds

(3.15) ε ≤ U ε , V ε ≤ |U 0 | L ∞ (Ω) + ε. For u ε = U α ε , v ε = V α ε , for a cut-off function ϕ ∈ C ∞ c
(Ω) satisfying ϕ ≥ 0 and for p > 1, we will compute d dt Ω ϕ 2 (v ε -u ε ) p + . We will use the positivity properties of the matrix Q α,p , to control the error term created by the cut-off. To this end, it is useful to introduce for δ ∈ [0, ε), the function Ψ δ,p (w) which is an approximation of Ψ 0,p (w) = w p + , given by

Ψ δ,p (w) = w p p -δ p-1 w -δ p 1 p -1 • 1 {w>δ} .
Then a direct computation gives (3.16)

d dt Ω ϕ 2 Ψ δ,p (v ε -u ε ) dx = Ω ϕ 2 Ψ δ,p (v ε -u ε ) {v γ ε ∆v ε -u γ ε ∆u ε } dx + Ω ϕ 2 Ψ δ,p (v ε -u ε ) γ v γ-1 ε |∇v ε | 2 -u γ-1 ε |∇u ε | 2 dx = {vε-uε>δ} f dx - {vε-uε>δ} ϕ 2 e[v ε , u ε ] dx + δ p-1 {vε-uε>δ} ϕ 2 g dx with the notation e[•, •] defined in (2.8) and f = -Ψ δ,p (v ε -u ε ) 2ϕ∇ϕ•{v γ ε ∇v ε -u γ ε ∇u ε } and g = (α -1) α |∇v ε | 2 v γ-1 ε -|∇u ε | 2 u γ-1 ε .
Here in (3.16) we have used equation (1.5) to get the second line, and have done an integration by parts on the first term of the second line, using in particular for δ > 0 the chain rule ∇Ψ δ,p (w) = Ψ δ,p (w) ∇w with equality almost everywhere with the convention that the right hand side is zero if ∇w = 0 irrespective of whether Ψ δ,p is defined. The case δ = 0 can be recovered, passing to the limit δ → 0. Using bounds (3.15), we see that 0 < v γ ε , u γ ε ≤ M for some constant M > 0. Moreover for v ε -u ε > δ > 0, we have by Young inequality

(3.17) |f | ≤ (v ε -u ε ) p-2 2M εν (v ε -u ε ) 2 |∇ϕ| 2 + ϕ 2 εν(v γ ε |∇v ε | 2 + u γ ε |∇u ε | 2 )
where ν = ν(δ) > 0 is given in (3.12) for δ > 0.

Step 2: First integral estimate on the gradient We now choose u ε ≡ ε and δ = 0. Then estimate (3.17) still holds true, but with ν replaced by ν 0 given in (3.13). This implies

d dt Ω ϕ 2 Ψ 0,p (v ε -ε) dx ≤ Ω C 0 |∇ϕ| 2 Ψ 0,p (v ε -ε) dx -(1 -ε) {vε-ε>0} ϕ 2 e[v ε , ε] dx with C 0 = 2pM εν 0 .
We then take a sequence of functions ϕ converging towards ϕ λ (x) = e -λ|x| for λ > 0. Integrating in time, and passing to the limit λ → 0, and using again (3.13), we get

(3.18) Ω Ψ 0,p (v ε (t) -ε) dx + (1 -ε)ν 0 t 0 dτ Ω v γ ε |∇v ε | 2 ≤ Ω Ψ 0,p (v ε (0) -ε) dx
Step 3: Refined estimate We come back to general v ε , u ε as in Step 1 and now consider δ > 0. Using (3.17), we get

d dt Ω ϕ 2 Ψ δ,p (v ε -u ε ) dx ≤ Ω C|∇ϕ| 2 Ψ 0,p (v ε -u ε ) dx -(1 -ε) {vε-uε>δ} ϕ 2 e[v ε , u ε ] dx + δ p-1 {vε-uε>δ} ϕ 2 g dx with C = 2pM εν
. Note that estimate (3.18) for both v ε and u ε , controls uniformly in time

Ω Ψ 0,p (v ε -u ε )
dx and the time integral of Ω δ p-1 |g| dx. Therefore, we can again apply the choice of ϕ = ϕ λ , integrate in time, and take the limit λ → +∞ as in Step 2, and then conclude in the limit δ → 0 (3.19)

Ω Ψ 0,p (v ε (t) -u ε (t)) dx+(1-ε) t 0 dτ {vε-uε>0}∩{uε>0} e[v ε , u ε ] dx ≤ Ω Ψ 0,p (v ε (0) -u ε (0)) dx.
Step 4: The limit ε → 0 For a point P 0 = (t 0 , x 0 ), let us denote the open parabolic cylinder Q r (P 0 ) = (t 0 -r 2 , t 0 ) × B r (x 0 ), where B r (x 0 ) is the open ball of center x 0 and radius r > 0. It is known by Theorem 1.1 in Sacks [START_REF] Sacks | Continuity of solutions of a singular parabolic equation[END_REF] that for any smooth solution U ε of (1.1) on Q 2r (P 0 ), there exists a modulus of continuity ω of U ε on Q r (P 0 ), depending only on r, d, n and |U ε | L ∞ (Q 2r (P 0 )) . This property is automatically transfered to u ε = U α ε with the modulus of continuity ω α . This implies, by Ascoli-Arzela theorem, that U ε → U , where U is still a solution of (1.1) with initial data U 0 in the sense of Definition 1.1. We also note that (with αp ≥ 1)

Ω |u(t, •) -u(s, •)| p ≤ Ω |U (t, •) -U (s, •)| αp ≤ (2|U | L ∞ (Q) ) αp-1 |U (t, •) -U (s, •)| L 1 (Ω)
which shows in particular that U ∈ C([0, +∞); L 1 (Ω)) implies u ∈ C([0, +∞); L p (Ω)).

From the standard parabolic theory [START_REF] Ladyženskya | Linear and Quasilinear Equations of Parabolic Type[END_REF], we know that U ∈ C ∞ (Q ∩ {U > 0}), with corresponding quantitative estimates on solutions locally bounded from above and below. Similarly, we have V ε → V , and we call u = U α , v = V α . Using the C 1 convergence of (v ε , u ε ) to (v, u) on compact sets inside {v -u > 0} ∩ {u > 0}, we can pass to the limit in (3.19), and get the same inequality for ε = 0 with (v ε , u ε ) replaced by (v, u). This means (2.7). We finally conclude to the result for general initial data U 0 , V 0 , by a standard approximation argument in L 1 (Ω).

Decay of the gradient

This section is divided into two subsections. In the first subsection, we give the proof of the gradient decay Theorem 2.2. In the second subsection, we give a directional derivative estimate (Theorem 4.2) as a corollary of our contraction estimate.

Decay of the gradient

We start with the following simpler analogue of Lemma 3.1, whose proof follows from an elementary computation. where I 2 is the 2 × 2 identity matrix, M α,p is defined in (2.9), and ν 1 = ν 1 (n, α, p) > 0.

Proof of Theorem 2.2

We follow the lines of proof of Theorem 2.1. For δ > 0, a direct computation gives with Here we have integrated by parts the term ∆u ε in the second line of (4.21). We then pass to the limit δ → 0. We replace the fine property (3.12) of the matrix Q α,p , by (4.20). We finally conclude as in the proof of Theorem 2.1.

b ε = ∇u ε |∇u ε | (4.21) d dt Ω ϕ 2 Ψ δ,p (|∇u ε |) dx = Ω ϕ 2 Ψ δ,p (|∇u ε |) b ε • ∇ u γ-1 ε u ε ∆u ε + γ|∇u ε | 2 dx = {|∇uε|>δ} f dx -

Decay of the directional derivative

For ξ = (ξ t , ξ x ) ∈ R × R d , we define the directional derivative (4.22) D ξ u = ξ t u t + ξ x • ∇u.

In this subsection, we prove the following result as a consequence of our contraction estimate.

Figure 1 :

 1 Figure 1: The set K |n| .

Theorem 2 . 1 (

 21 The new contraction family) We recall that Ω = R d with d ≥ 1. Let |n| < 1 and U (resp. V ) be the solution of (1.1) in the sense of Definition 1.1, with initial data 0

Lemma 4 . 1 (

 41 Properties of the matrix M α,p ) For |n| < 1, let (α, p) ∈ Int K |n| where K |n| is defined in (1.4). Then we have (4.20) M α,p ≥ ν 1 I 2 ,

{|∇uε|>δ} ϕ 2 1 ε

 21 ē[u ε , ∇u ε , D 2 u ε ] dx + δ p-1 {|∇uε|>δ} ϕ 2 ḡ with the notation ē[•, •, •] defined in (2.11) and f = d i=1 -2ϕ∇ i ϕΨ δ,p (|∇u ε |) b ε • ∇ {u γ ε ∇ i u ε } and ḡ = (α -1) α b ε • ∇ u γ-|∇u ε | 2 .

Theorem 4 . 2 (

 42 Decay of the directional derivative)Let U be the solution of equation (1.1) with initial data U 0 , under the assumptions of Theorem 2.1. We also set u = U α . Let ξ ∈ R × R d , and let us assume that

( 4 .

 4 23) C ξ := lim sup η→0 + Ω u((0, x) + ηξ) -u(0, x) η p dx < +∞ With the notation D ξ u defined in (4.22), we have D ξ u ∈ L ∞ ([0, +∞); L p (Ω)). Moreover, with the notation u(τ ) = u(τ, •), we have for almost every t ≥ 0: (4.24) Ω |D ξ u(t)| p dx + p t 0 {D ξ u =0}∩{u>0} ē[u, D ξ u, ∇u]dx ≤ C ξ .

In particular (4.23) is satisfied if ξ t = 0 and D ξ U α 0 ∈ L p (Ω). We will use the following result. Lemma 4.3 (Limit of the quadratic term) For |n| < 1, let (α, p) ∈ K |n| where K |n| is defined in (1.4).

Let Ω 1 be an open set such that u, w, z ∈ C 1 (Ω 1 ) with u > 0 on Ω 1 . For η ∈ (0, 1), let v η be a function satisfying

Then we have on Ω 1

where we recall that the matrices Q α,p and M α,p are respectively defined in (2.6) and (2.9).

Proof of Lemma 4.3

The proof is done using a simple Taylor expansion argument. For Q α,p (v, u), we set

which implies the result.

Proof of Theorem 4.2 Let v η (t, x) = u((t, x) + ηξ). We simply apply Theorem 2.1 to v η -u and u -v η and use the symmetry e[v, u] = e[u, v], in order to control

We pass to the limit as η → 0 using Lemma 4.3 which gives (4.23).

Aknowledgements

The authors thank B. Andreianov for helpful indications on the literature.