Gabriele Perozzi
email: gabriele.perozzi@onera.fr

Denis Efimov
email: denis.efimov@inria.fr

Jean-Marc Biannic
email: jean-marc.biannic@onera.fr

Laurent Planckaert
email: laurent.planckaert@onera.fr

Trajectory Tracking For A Quadrotor Under Wind Perturbations: Sliding Mode Control With State-Dependent Gains

Keywords: Quadrotor modeling, sliding mode control, wind rejection

The problem of position tracking of a mini drone subject to wind perturbations is investigated. The solution is based on a detailed unmanned aerial vehicle (UAV) model, with aerodynamic coefficients and external disturbance components, which is introduced in order to better represent the impact of the wind field. Then, upper bounds of wind-induced disturbances are characterized, which allow a sliding mode control (SMC) technique to be applied with guaranteed convergence properties. The peculiarity of the considered case is that the disturbance upper bounds depend on the control amplitude itself (i.e. the system is nonlinear in control), which leads to a new procedure for the control tuning presented in the paper. The last part of the paper is dedicated to the analysis and reduction of chattering effects, as well as investigation of rotor dynamics issues. Conventional SMC with constant gains, proposed first order SMC, and proposed quasi-continuous SMC are compared. Nonlinear UAV simulator, validated through in-door experiments, is used to demonstrate the effectiveness of the proposed controls.

aerodynamic coefficients.

Introduction

In the last decades there has been a growing interest for small unmanned aerial vehicles both for military and commercial applications. They are useful in different domains such as monitoring, inspection, and other actions, especially in urban areas or nearby buildings and dangerous interiors, for natural calamities such as earthquakes. These machines are often required to move in unfamiliar environments in terms of geography and in terms of aerological conditions. In addition, the low mass of such units (comparing to the forces generated by the air disturbances) reduces significantly the domain of stable flight, then additional constraints have to be considered in the control design. Thus, it appears inevitable that, if we want to let UAVs operate in urban environments, inside turbulent air flow patterns for which accurate prediction is not possible a priori with limited resources, we need to focus on detailed aerodynamic models and sophisticated control laws. This paper is part of the four years project "Small drones in the wind" at ONERA Lille, which aims to use the drone itself as a wind sensor. First year was dedicated to the identification of the UAV model affected by wind velocity, which is used in this article and explained in [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF]. This control design is carried out in parallel with development of a wind estimation tool in [? ].

The considered problem consists in design of a nonlinear robust control law, which ensures a stable and efficient navigation of a small UAV under unpredictable wind perturbations. The model described in [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF] can also be used to estimate wind velocity, in this way the estimated values can themselves be used as inputs to the control to properly adapt the gains on-line and use this estimation in path planning and trajectory control to smartly avoid collisions.

There exist many control design techniques to counteract the effects of wind perturbations on flight of small UAVs, among which SMC plays a keyrole. Many methods have been proposed in the literature, for instance, some principal SMCs with their relative sliding surfaces and Lyapunov functions are illustrated in [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF], [START_REF] Polyakov | Stability notions and Lyapunov functions for sliding mode control systems[END_REF]. Its insensitivity to the model errors, parametric uncertainties and other disturbances and its ability to globally stabilize the underactuated systems are two advantages of the sliding mode controller [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF]. Sliding mode algorithms are extensively applied to dynamic systems and optimal algorithms are discussed in [START_REF] Basin | Optimal sliding mode algorithms for dynamic systems[END_REF]. Dozens of articles have applied SMC methodology to UAVs in order to solve the position and the attitude tracking problems ensuring robustness against external disturbances. They have been compared extensively with respect to other controls in in-door experiments such as in [START_REF] Bouabdallah | Backstepping and Sliding-mode Techniques Applied to an Indoor Micro Quadrotor[END_REF], where SMC is compared with backstepping control for micro quadrotor. Chattering-free sliding mode controller is proposed in [START_REF] González | Chattering-Free Sliding Mode Altitude Control for a Quad-Rotor Aircraft: Real-Time Application[END_REF] by replacing a sign function with a high-slope saturation function, in [START_REF] Sumantri | Least square based sliding mode control for a quad-rotor helicopter and energy saving by chattering reduction[END_REF] the energy saving effect because of chattering reduction is also evaluated. Second order sliding mode is used in [START_REF] Zheng | Second order sliding mode control for a quadrotor UAV[END_REF] where two different sliding manifolds are defined for fully actuated and underactuated subsystems. In the paper [START_REF] Ramirez-Rodriguez | Robust Backstepping Control Based on Integral Sliding Modes for Tracking of Quadrotors[END_REF] a robust backsteppingbased controller is proposed that induces integral sliding modes for the Newton-Euler underactuated dynamic model of a quadrotor subject to smooth bounded disturbances. The trajectory tracking of uncertain underactuated nonlinear dynamic systems is tackled by an adaptive fuzzy hierarchical slidingmode control in [START_REF] Hwang | Adaptive Fuzzy Hierarchical Sliding-Mode Control for the Trajectory Tracking of Uncertain Underactuated Nonlinear Dynamic Systems[END_REF]. In the work [START_REF] Xiong | Position and attitude tracking control for a quadrotor UAV[END_REF] the controller of the fully actuated subsystem using a robust terminal sliding mode control algorithm is designed. A controller to provide robust position and attitude of the vehicle while relying only on knowledge of the limits of the disturbances is proposed in [START_REF] Besnard | Quadrotor vehicle control via sliding mode controller driven by sliding mode disturbance observer[END_REF]. Sliding mode techniques are used also as observers estimators of the effect of the external disturbances such as wind and noise, and the whole observer-estimator-control law approach is presented in [START_REF] Benallegue | High-order sliding-mode observer for a quadrotor UAV[END_REF]. A SMC is proposed to stabilize a class of cascaded under-actuated systems, in which the UAV system is divided, in [START_REF] Xu | Sliding Mode Control of a Quadrotor Helicopter[END_REF]. Famous super-twisting algorithm, which is able to ensure robustness with respect to bounded external disturbances, is discussed in [START_REF] Derafa | Super twisting control algorithm for the attitude tracking of a four rotors UAV[END_REF]. In [START_REF] Perozzi | On sliding mode control design for UAV using realistic aerodynamic coefficients[END_REF], the authors propose control algorithms similar to the ones described in this note, however they do not consider the rotors dynamics, the 3D trajectory, the first and second derivative of the UAV position and angles, making their approach interesting only in a mathematical point of view but impracticable in real 3D trajectory flight and far from the hover condition.

The quadrotor model itself is well established in the literature under various simplifications, in the present note we investigate the influence of the wind taking into account an identification work for the aerodynamic coefficients, which is useful for the control design. Two robust nonlinear SMC law design are described, which consider realistic assumptions on external disturbances of quadrotors. In the considered case the upper bound of matched disturbances depends nonlinearly on the control itself, the system state vector and wind disturbances. The closed-loop system stability is ensured for a selected maximum admissible value of the wind speed. The control strategy proposed in this article can be equipped with an additional wind estimator algorithm, as in [START_REF] Perozzi | Wind estimation algorithm for quadrotors using detailed aerodynamic coefficients[END_REF], reducing automatically the control effort on the rotors when it is not necessary. However the peculiarity of the proposed SMC approach design is that the control allows the UAV to remain stable even without any coupled external disturbance observer.

The paper outline is as follows. In Section 2 the considered UAV drone is described and the flight dynamics model is derived. The control and disturbance bounds are calculated in Section 3. The two control designs are presented in Section 4, where the rotors issue, the chattering, the physical constraints and tools used for the implementation are also illustrated. The results of numerical experiments are shown in Section 5. The remarks and discussion conclude the paper in Section 6.

Quadrotor model

The presented work is based on the commercial Parrot Ar. Drone 2.0 having configuration as it is shown in Fig. 1. This section presents a detailed model of the UAV dynamics used to develop the controller and to estimate the wind velocity, making the drone as a wind sensor.

Quadrotor flight dynamics

According to the identification work for the Parrot AR Drone 2.0 in [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF] performed at maximal speeds in forward, lateral, and vertical directions respectively ±5, ±5 ± 1 m/s (as in Table 3), rotor gyroscopic effects and inertial counter torques can be neglected since they are rather small. Then the translation dynamics of the drone in the body frame yield

F Xaero F Y aero F Zaero 3 5 + mR T 2 4 0 0 g 3 5 , (1)
which can be represented in the earth frame

m 2 4 ẍ ÿ z3 5 = R 2 4 F Xaero F Y aero F Zaero 3 5 + m 2 4 0 0 g 3 5 , ( 2
)
where m is the mass of the UAV, u, v, w are its linear velocities expressed in body frame, p, q, r are its angular velocities in body frame, F Xaero , F Y aero , F Zaero are the external aerodynamic forces in body frame, g is is the gravity acceleration in inertial frame, R is the rotational matrix defined as

R = 2 4 c c ✓ s c + c s ✓ s s s + c s ✓ c s c ✓ c c + s s ✓ s c s + s s ✓ c s ✓ c ✓ s c ✓ c 3 5 ,
where c = cos( ), s = sin( ) and similar for the other angles. Then the passage from the earth frame (R 0 ) to the body one (R) is

[X T ] R = [X T ] R 0 R. (3)
The rotational dynamics of the drone with respect to inertial earth frame are

I 2 4 ṗ q ṙ3 5 = 2 4 p q r 3 5 ⇥ I 2 4 p q r 3 5 + 2 4 L aero M aero N aero 3 5 , (4)
with the inertia matrix

I = 2 4 I xx 0 0 0 I yy 0 0 0 I zz 3 5 ,
where L aero , M aero , N aero are the external aerodynamic moments in the body frame.

The relation between angular velocities and Euler angles ˙ = p + tan ✓(q sin + r cos ), ✓ = q cos r sin , ˙ = q sin + r cos cos ✓ ,

are considered avoiding the singularity at ✓ = ⇡ 2 , which is a reasonable assumption in our case since the topic of this article is not to achieve extreme maneuvers.

Hence the full model of the system is presented by the equations ( 2), ( 4), [START_REF] Basin | Optimal sliding mode algorithms for dynamic systems[END_REF]. The aerodynamic forces F Xaero , F Y aero , F Zaero , moments L aero , M aero , N aero , and related coefficients are derived below using a combination of momentum and blade element theory in helicopters, well explained in [START_REF] Bramwell | Bramwell's Helicopter Dynamics[END_REF], [START_REF] Johnson | Helicopter Theory[END_REF], [START_REF] Leishman | Principles of Helicopter Aerodynamics with CD Extra[END_REF].

Quadrotor aerodynamics

Aerodynamic forces and moments for each rotor, where subscript j indicates the j th rotor, are derived as

F X j = ⇢AR 2 u j u w p (u j u w ) 2 + (v j v w ) 2 C H j ! 2 j , F Y j = ⇢AR 2 v j v w p (u j u w ) 2 + (v j v w ) 2 C H j ! 2 j , F Z j = ⇢AR 2 C T j ! 2 j , L j = sign ! j ⇢AR 3 u j u w p (u j u w ) 2 + (v j v w ) 2 C Rmj ! 2 j , M j = sign ! j ⇢AR 3 v j v w p (u j u w ) 2 + (v j v w ) 2 C Rmj ! 2 j , N j = sign ! j ⇢AR 3 C Q j ! 2 j
, where ⇢ is the air density, A is the rotor area, R is the rotor radius, (u w , v w , w w ) are the wind velocities with respect to the earth in body frame respectively in (x, y, z) directions, C H j is the hub force coefficient, C T j is the rotor thrust coefficient, ! j is the rotor angular speed, C Q j is the rotor drag moment coefficient, C Rmj is the rotor rolling moment coefficient.

Total aerodynamic forces are

F Xaero = 4 X j=1 F X j , F Y aero = 4 X j=1 F Y j , F Zaero = 4 X j=1 F Z j . (6)
Total aerodynamic moments are

L aero = 4 X j=1 (L j + F Z j ls j hF Y j ), M aero = 4 X j=1 (M j F Z j lc j + hF X j ), N aero = 4 X j=1 (N j + F Y j lc j F X j ls j ), (7)
where h is the distance between rotors plane and the center of gravity of the UAV, l is the arm length, and with

c j = cos ⇡ 2 (j 1) + ✏ , s j = sin ⇡ 2 (j 1) + ✏ ,
where in our UAV configuration we have ✏ = ⇡ 4 . Thus, for vectors c j and s j we have cosines and sinus of the angles

[ ⇡ 4 , 3 4 ⇡, 5 4 ⇡, 7 4 ⇡].
The blade flapping can be modeled as a fictitious hinge spring.

The aerodynamic coefficients (C T coefficient of the thrust, C Rm coefficient of the rolling moment, C H coefficient of the hub force, C Q coefficient of the drag moment, µ advance ratio, inflow ratio) are defined as follows

C Rmj = a µ j 8 ( j 4 3 ✓ 0 + ✓ tw ) + b 1 16 (1 µ 2 j 2 ) , C T j = a (1 + 3 2 µ 2 j ) ✓ 0 6 (1 + µ 2 j ) ✓ tw 8 j 4 , j = a (1 + 3 2 µ 2 j ) ✓ 0 6 (1 + µ 2 j ) ✓tw 8 j 4 2 q µ 2 j + 2 j + w w w j + qlc j pls j R|! j | , µ j = 1 R|! j | q (u j u w ) 2 + (v j v w ) 2 , C Q j = C Q P j + C Q ij , C Q P j = 1 8 (C D0 + C Di ✓ 2 0 )(1 + µ 2 j ) C Di ✓ 0 ✓ tw ( 1 5 + µ 2 j 6 ) + C Di ✓ 2 tw ( 1 12 + µ 2 j 16
)

C Di j ( ✓ 0 3
✓ tw 4 ) + C Di ⇣ µ 2 j 8 (a 2 0 + a 2 1 4 + 3b 2 1 4 ) + 1 16 (a 2 1 + b 2 1 ) + 2 j 4 + µ j 6 a 0 a 1 j µ j 4 b 1 ⌘ , C Q ij a = j ( ✓ 0 6 ✓ tw 8 j 4 ) µ 2 j 8 (a 2 0 + a 2 1 4 + 3b 2 1 4 ) 1 16 (a 2 1 + b 2 1 ) µ j 6 a 0 a 1 + j µ j 4 b 1 , C H j = C H P j + C H ij , C H P j = µ j 4 (C D0 + C Di ✓ 2 0 ) + C Di µ j 24 (3✓ 2 tw 8✓ 0 ✓ tw ) + ✓ 0 24 (3µ 2 j b 1 12 j µ j 4b 1 ) ✓ tw 16 (µ 2 j b 1 4 j µ j 2b 1 ) + µ 2 j 8 a 0 a 1 + µ j 16 (a 2 1 b 2 1 ) + 1 4 j b 1 , C H ij a = ✓ 0 4 ( j µ j + 2 3 b 1 ) ✓ tw 8 ( j µ j + b 1 ) + µ j 8 (a 2 0 + b 2 1 ) 3 8 j b 1 + 1 12 a 0 a 1 , (8)
where w w is the wind velocity with respect to the earth in body frame in z direction, is the rotor solidity ratio, a is the lift curve slope of the blade section, C D0 is the drag coefficient of the blade section, C Di is the induced drag coefficient of the blade section, ✓ 0 is the angle of attack of the root profile, a 0 , a 1 , b 1 are the coefficients of the blade flapping equation, is the inflow ratio, µ is the advance ratio, subscript stat indicates the value in stationary phase. The shape of above UAV coefficients can be explained recalling aerodynamic science: the thrust is the resultant of the vertical forces acting on all the blade elements. The inflow ratio is the ratio between the component of UAV velocity perpendicular to the rotor disk with respect to the blade tip speed. The advance ratio indicates the component of the UAV velocity parallel to the rotor disk with respect to the blade tip speed. The rolling moment of a propeller exists in forward flight when the advancing blade is producing more lift than the retreating one and it is the integration over the entire rotor of the lift of each section acting at a given radius. The hub force is the resultant of the horizontal forces acting on all the blade elements.

The drag moment about the rotor shaft is caused by the aerodynamic forces acting on the blade elements, the horizontal forces acting on the rotor are multiplied by the moment arm and integrated over the rotor. Identification results thanks to indoor experiments on Parrot drone at low/moderate velocity, validates the previous UAV model and allows some simplifications to be accepted:

✓ tw = 0. C Di = 0. j = stat 4 a K z w j ww R|! j | , K z , stat 0
, where K z (as in Table 1 and by [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF]) comes from the approximation of the equation in vertical flight.

C T j = C T stat + K z w j ww R|! j | , C T stat
0, such models of j and C T j are rather precise in the climbing phase, but less accurate in descent phase, since the model tends to slightly overestimate the propulsion in the descent phase. C H j = K D µ j , K D 0, the UAV drag is modeled as ⇢AR 2 P C H j ! 2 j , corresponding to the rotors at low UAV speed, otherwise at higher speed we need to add the body drag effect since it depends on square of velocity. However the constant K D , (as in Table 1 and by [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF]) has been identified for the forward velocity less than 5 m/s taking into account the interactions between the rear and the front rotors, and considering the whole UAV body and rotors.

Rotors dynamics

The rotors are driven by DC-motors, which binds electrical and mechanical quantities. As in [START_REF] Bouabdallah | Backstepping and Sliding-mode Techniques Applied to an Indoor Micro Quadrotor[END_REF] [START_REF] Bouabdallah | Full control of a quadrotor[END_REF], rotors can be represented as in Fig. 2, and described by the dynamic system

( L di dt = u Ri K e ! m J d!m dt = ⌧ m ⌧ d ,
where u is the motor input, R is the motor inertial resistance, L is the motor inductance, K e is the motor constant corresponding to the back EMF constant, ⌧ m is the motor torque, ! m is the motor angular speed, ⌧ d is the motor load, J is the motor moment of inertia, i is the current. Considering a small motor with a very low inductance the given system can be simplified and

G(s) = 1 bs + 1 , (9)
where the time constant b = 0.1 can be found considering Parrot Drones. Time constant slightly varies based on the considered motor and thus on the considered UAV mass. Precisely, in the Parrot Drones, the time constant is less than 0.1 for the increasing rotor angular velocities, and it is almost equal to 0.1 for decreasing rotor angular velocities. In this section, two UAV models with different complexity have been presented (a complete one, and a reduced model for low/moderate velocity). The model with simplifications, which are accepted following [START_REF] Planckaert | Quadrotor UAV aerodynamic model identification using indoor flight experiment and feasibility of UAV as wind gust sensor[END_REF], will be used to design a SMC, the other one will be used to validate the control in numerical experiments, since it considers nonlinear coefficients and a more realistic UAV behavior.

Control system equations

The model presented in previous section can be rewritten in the state-space form

Ẋ = f (X, U, d),
where f is expressed in ( 1), ( 4), ( 5) and the state vector X is chosen as

X = [x y z ẋ ẏ ż ✓ p q r] T , U = [U z U ✓ U U ]
T is the control input, and disturbances d are described in the next subsection. The relations between the control inputs and the rotor

6 6 4 U z U ✓ U U 3 7 7 5 ⇡ 2 6 6 4
K f K f K f K f K f lc j K f lc j K f lc j K f lc j K f ls j K f ls j K f ls j K f ls j K m K m K m K m 3 7 7 5 2 6 6 4 ! 2 1 ! 2 2 ! 2 3 ! 2 4 3 7 7 5 , with ! min  ! j  ! max , where K f = ⇢AR 2 C T stat and K m = ⇢AR 3 C D0 8 + stat a( ✓ 0 6 stat 4
) . In this work, the control inputs are selected to be proportional to the terms with ! 2 j . Thus, expanding ( 1) and ( 4), the other terms dependent linearly on ! j and wind velocities are considered as disturbances. Since we do not know in advance the wind perturbations, then we cannot use these terms in controls. Such a decomposition of thrust (which is proportional to ! j and ! 2 j ) and selection of disturbances are almost exact in the hover flight, where we have (p, q, r) ⇡ ( ˙ , ✓, ˙ ).

In the following, linear rotors speed in body frame are computed as a function of the state in body frame 2 4 u j v j w j

Definition of disturbances

In the literature, it is often assumed that each component of the disturbance input vector d admits a fixed upper bound, which means |d|  D for some known D 0. Unfortunately, it is a rather conservative hypothesis, and that is why the varying state-dependent bounds will be considered in our case for d. However, we will assume boundedness of the wind velocities:

|u w |  D x , |v w |  D y , |w w |  D z , for some known D x 0, D y 0, D z 0, which is a reasonable restriction.
In this work we will call signals, whose influence is compensated by robust abilities of the control, as "disturbances", which include external inputs (the wind velocities), state components and controls. Substituting ( 8), ( 6) and ( 7) in ( 1), ( 4

d x = 4 X j=1 ⇢ARK D (u j u w ) |! j |, (10)
d y = 4 X j=1 ⇢ARK D (v j v w ) |! j |, (11)
d z = 4 X j=1 ⇢ARK z (w j w w ) |! j |, (12)
d = 4 X j=1 ⇣ ! j ⇢AR 2 (u j u w ) a 2 ( ✓ 0 3 stat 4 ) + |! j |⇢AR hK D (v j v w ) lK z (w j w w ) s j + sign (! j ) 1 2 ⇢ARK z (u j u w ) (w j w w ) ⌘ , (13)
d ✓ = 4 X j=1 ⇣ ! j ⇢AR 2 (v j v w ) a 2 ( ✓ 0 3 stat 4 ) + |! j |⇢AR hK D (u j u w ) + lK z (w j w w ) c j + sign (! j ) 1 2 ⇢ARK z (v j v w ) (w j w w ) ⌘ , (14)
d = 4 X j=1 ✓ ! j ⇢AR 2 K z (w j w w ) ( 2✓ 0 3 2 stat ) |! j |⇢lARK D (v j v w ) c j (u j u w ) s j sign (! j ) ⇢AR ⇣ C D0 8 (u j u w ) 2 + (v j v w ) 2 4 a K 2 z (w j w w ) 2 ⌘ ◆ . (15)

Disturbance upper bounds

To design a control, which is able to compensate the disturbances, we have to evaluate the upper bounds for them.

For x dynamics, the upper bound in body frame is computed from equation ( 10)

|d x | = | KD 4 X j=1 (u j u w ) ! j |  KD 4 X j=1 |u j u w ||! j  KD 4 X j=1 (|u j | + |u w |) |! j |  KD (max j |u j | + |u w |) 4 X j=1 |! j |,
where KD = ⇢ARK D . Using the control U z = K f P 4 j=1 ! 2 j and applying the Jensen's inequality, an upper estimate can be obtained

4 X j=1 |! j |  K p |U z |, K = 2 p K f . ( 16
)
An upper bound of the disturbance becomes

|d x |  KD (|X| + D x ) p |U z | = d xx , (17)
where KD = K KD . In the earth frame using the rotation matrix (3) and the upper bound for the disturbances ( 17), ( 18), ( 19) it becomes

|d xe | = c c ✓ d x + (c c ✓ c s c )d y + (s s + c s ✓ c )d z  d xx + 2d yy + 2d zz ,
where d y , d z , d yy and d zz are derived below. For the disturbances d ye and d ze the computations are similar to the previous ones, therefore only final expressions are given next. For y dynamics, the upper bounds in body and earth frames respectively are computed from the equation ( 11)

|d y |  KD (|X| + D y ) p |U z | = d yy , ( 18
)
|d ye |  d xx + 2d yy + 2d zz .
For z dynamics, the upper bounds in body and earth frames respectively are computed from the equation ( 12)

|d z |  Kz (|X| + D z ) p |U z | = d zz , (19)
|d ze |  (f ze (|X|) + D ze ) p |U z |. ( 20
)
where

Kz = K⇢ARK z , f ze (|X|) = Kz max j |w j |+ KD max j |u j |+ KD max j |v j |, D ze = Kz D z + KD (D x + D y ).
For roll dynamics, the upper bound estimate is computed from the equations ( 13) and ( 16)

|d |  K (f 1 (X) + D 1 ) p |U z | + K (f 2 (X) + D 2 ) (21)
where

K 1 = ⇢AR 2 a 2 ( ✓ 0 K 2 D y +K 3 D z max j |s j |, K = 1 2 ⇢ARK z , f 2 (X) = max j |u j | 2 +max j |w j | 2 , D 2 = D 2
x + D 2 z . For pitch dynamics, the upper bound estimate is computed from the equations ( 14) and ( 16)

|d ✓ |  K✓ (f ✓1 (X) + D ✓1 ) p |U z | + K✓ (f ✓2 (X) + D ✓2 )
where

K ✓1 = ⇢AR 2 a 2 ( ✓ 0 3 stat 4 ), K ✓2 = ⇢ARhK D , K ✓3 = lK z ⇢AR, K✓ = K, f ✓1 (X) = K ✓1 max j |v j | + K ✓2 max j |u j | + K ✓3 max j |w j c j |, D ✓1 = K ✓1 D y + K ✓2 D x + K ✓3 D z max j |c j |, K✓ = 1 2 ⇢ARK z , f ✓2 (X) = max j |v j | 2 + max j |w j | 2 , D ✓2 = D 2 y + D 2 z .
For yaw dynamics, the upper bound estimate is computed from the equations ( 15) and ( 16)

|d |  K (f 1 (X) + D 1 ) p |U z | + K (f 2 (X) + D 2 )
where

K 1 = ⇢AR 2 K z ( 2✓ 0 3 2 stat ), K 2 = ⇢ARlK D , K 3 = ⇢ARlK D , K = K, f 1 (X) = K 1 max j |w j | + K 2 max j |v j c j | + K 3 max j |u j s j |, D 1 = K 1 D z + K 2 D y max j |c j | + K 3 D x max j |s j |, K = ⇢AR, f 2 (X) = C D0 8 ⇥ (max j |u j | 2 + max j |v j | 2 ) + 4 a max j |w j | 2 , D 2 = C D0 8 (D 2 x + D 2 y ) + 4 a D 2 z .

Controls design

For synthesis of control laws the SMC approach will be applied. First order and quasi-continuous SMCs (respectively 1-SMC, qc-SMC) will be designed in this section.

First order control design

This control methodology takes into account and compensates the matched disturbances. The big issue for the considered problem is that the disturbance d depends on wind signals, the control itself, and state of the system, as shown above. Thus, a mild development of SMC approach is needed. To this end, the sliding surfaces in this work are selected proportional to the regulation errors e i , in this way we are going to control the dynamics proportional to position and velocity

S i = ėi + ↵ i e i , ↵ i > 0, (22)
where i 2 (x, y, z, , ✓, ). The Lyapunov function is then chosen as

V i = 1 2 S 2 i . (23)
The dynamics of z can be rewritten in the earth frame

z = g (cos cos ✓) 1 m (U z + d ze ). (24)
In hover state we have cos ✓ cos ⇡ 1, thus with a rotation of the UAV, a reasonable assumption is that | cos cos ✓| > 0, where is our operating point limit. To build the altitude control, the regulation error has been chosen as

e z = z z des , (25)
where z des is the desired altitude for UAV. Thus, using [START_REF] Perruquetti | Finite-Time Observers: Application to Secure Communication[END_REF] and its derivative in [START_REF] Bouabdallah | Full control of a quadrotor[END_REF], the derivative of the sliding surface can be written as follows Ṡz = z zdes + ↵ z ėz and using [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF], it also equals to

Ṡz = g cos ✓ cos m (U z + d ze ) zdes + ↵ z ėz . ( 26
)
The following expression of the control can be selected

U z = m cos ✓ cos (g ũz zdes + ↵ z ėz ) , (27)
where ũz is an auxiliary control defined later. After substitution of the control (27) in (26) we obtain

Ṡz = ũz + d ze cos ✓ cos m . ( 28
)
Using the Lyapunov function ( 23) with (28), its derivative takes the form

V =S z ũz + S z d ze cos ✓ cos m  S z ũz + |S z ||d ze 1 m |.
Then, using [START_REF] Johnson | Helicopter Theory[END_REF], the Lyapunov derivative becomes

V  S z ũz + |S z | 1 m (f ze (|X|) + D ze ) p |U z |.
From the expression (27) we can derive

|U z |  m (|g zdes + ↵ z ėz | + |ũ z |), (29)
and then, substituting in Lyapunov derivative, we obtain

V  S z ũz + |S z | %(X, zdes , żdes ) + ⌫(X) p |ũ z | , (30)
where

%(X, zdes , żdes ) = 1 p m (f ze (|X|) + D ze ) ⇥ p |g + ↵ z ėz zdes |, ⌫(X) = 1 p m (f ze (|X|) + D ze ) .
For simplicity of notation, in the following part of the article we will denote %(X, zdes , żdes ) = %(X). Let us look for the control in the form

ũz = (X)sign(S z ),
where (X) is a function to synthesize. Substituting this control in the obtained inequality (30) we get

V  |S z | %(X) + ⌫(X) p (X) (X)
and it is necessary to ensure by a choice of (X)

%(X) + ⌫(X) p (X) (X) = .
Solving this quadratic inequality with respect to (X) we found

(X) = 1 2 ⌫(X) 2 + 2%(X) + ⌫(X) p ⌫ 2 (X) + 4%(X) + , (31)
where > 0 is a tuning parameter. Substituting (31) in the inequality for

V we get V < p 2 p V ,
consequently, the system state trajectory reaches and stays on the sliding surface S z , which means asymptotic convergence of e z to the origin. Moreover, finite-time stability with respect to S z can be proved according to [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF].

The dynamics of x can be rewritten in the earth frame as follows ẍ = (sin sin + cos sin ✓ cos ) U z m d xe .

Since x dynamics cannot be directly controlled, a virtual control U x = sin sin + cos sin ✓ cos is introduced, which will be used to find the desired Euler angles des , ✓ des that will be the inputs of attitude controller next. Then the dynamics of x takes the form

ẍ = U x U z m d xe .
In order to guarantee negativity of the time derivative of the Lyapunov function ( 23) the following expression of control in x is selected

U x = m U z (ũ x + ẍdes ↵ x ėx ) , (32)
where ũx = (d xx + 2d yy + 2d zz ) signS x is an auxiliary SMC for x dynamics. Such a design is admissible since the control U z is always separated with the zero and the discontinuous term in ũx , which leads to the chattering effect, is replaced with the continuous saturation approximation.

The y dynamics, by introducing the virtual control U y = sin sin ✓ cos cos sin , can be rewritten in the earth frame as follows ÿ = (sin sin ✓ cos cos sin )

U z m d ye = U y U z m d ye .
Using similar arguments as for x, the auxiliary SMC for y dynamics takes the form ũy = (d xx + 2d yy + 2d zz ) signS y , with the following expression of control in y

U y = m U z (ũ y + ÿdes ↵ y ėy ) . ( 33
)
The dynamics of can be rewritten as

ṗ = qr I yy I zz I xx + 1 I xx (U + d ). (34)
To build the roll control U , the error has been chosen as

e = des . (35)
Using (35) and substituting its derivative in [START_REF] Bouabdallah | Full control of a quadrotor[END_REF], the derivative of the sliding surface can be calculated

Ṡ = ṗ ¨ des + ↵ ė ,
which due to (34) also equals to

Ṡ = qr I yy I zz I xx + 1 I xx (U + d ) + ↵ ė ¨ des . (36)
The following expression of the control can be selected

U = I xx ✓ qr I yy I zz I xx + ũ ↵ ė + ¨ des ◆ ,
where ũ is an auxiliary SMC. After substitution of this control in (36), we obtain

Ṡ = ũ + d I xx . (37)
Using the Lyapunov function ( 23) with (37) and ( 21), its derivative becomes

V =S ũ + S d  S ũ + |S ||d 1 I xx |  S ũ + |S | I xx ⇥ ⇣ K f 1 (X) + D 1 p |U z | + K f 2 (X) + D 2 ⌘ . ( 38
)
To have (38) negative, the auxiliary control for dynamics must be

ũ = 1 I xx signS ⇣ K f 1 (X) + D 1 p |U z | + K f 2 (X) + D 2 ⌘ .
Pitch and yaw controls can be designed following computations similar to the roll one, so only final expressions are given below.

The dynamics of ✓ can be rewritten as

q = pr I zz I xx I yy + 1 I yy (U ✓ + d ✓ ).
The auxiliary control for ✓ dynamics

ũ✓ = 1 I yy signS ✓ ⇣ K✓ f ✓1 (X) + D ✓1 p |U z | + K✓ f ✓2 (X) + D ✓2 ⌘ ,
with the following expression of control

U ✓ = I yy ✓ pr I zz I xx I yy + ũ✓ ↵ ✓ ė✓ + ✓des ◆ .
The dynamics of can be rewritten as

ṙ = pq I xx I yy I zz + 1 I zz (U + d ).
The auxiliary control for dynamics

ũ = 1 I zz signS ⇣ K f 1 (X) + D 1 p |U z | + K f 2 (X) + D 2 ⌘ ,
with the following expression of control

U = I zz ✓ pq I xx I yy I zz + ũ ↵ ė + ¨ des ◆ .
In the control design above it has been assumed that the rotors possess an immediate response on the desired values ! j assigned to them by the control law. In reality they admit some dynamics, and for validation and comparison of the proposed control strategy, the transfer functions for the rotors have to be introduced in the realistic form as in eq. ( 40). Since it is a dynamics not considered during design, then it adds an undesired delay to the control. Recall also that the SMC is built for a simplified model with identified coefficients to diminish the design problem complexity, however it will be tested with the complete nonlinear UAV model validated through indoor experiments at ONERA lab.

Since a big shortage of SMC is the chattering (high frequency oscillations of a discontinuous control signal in the steady-state mode caused by imprecision of values of model parameters, digital and measurement noises), which can ruin the rotors, then below several solutions for chattering reduction are compared. To this end, several modifications are introduced into the control algorithm for a comparison.

In the literature the problem of chattering reduction is a well-known issue discussed in many articles, see for an example [START_REF] González | Chattering-Free Sliding Mode Altitude Control for a Quad-Rotor Aircraft: Real-Time Application[END_REF]. Saturation functions are popular solutions used for chattering reduction in SMC that leads to a practical stability in the closed-loop system. Arc-tangent function replaces the function sign with a smooth function:

sat ⇠ (x) = ( sign(x) if |x| ⇠ 4 ⇡ arctan ⇣ x ⇠ ⌘ if |x| < ⇠ , ( 39
)
with the tunable gain ⇠ > 0. According to [START_REF] Khalil | Nonlinear Systems[END_REF], for a sufficiently small ⇠ if for a sign function all trajectories converge in a finite time to an equilibrium, then with the saturation all trajectories will converge to a compact set around that equilibrium with the size proportional to the value of ⇠. In practice ⇠ should be chosen small enough to find a trade-off between chattering reduction and minimal acceptable steady-state error.

Quasi-continuous control design

Another way to reduce the chattering is to use a high order sliding mode (HOSM) control [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF]. In this subsection for synthesis of a control law, a kind of HOSM control called the quasi-continuous SMC [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF] will be applied, which can be considered as an approximation of the sign on the plane (see the Appendix for the explanation of the idea of qc-SMC). For this design the transfer functions for the rotors are taken into account in a generic and realistic form as previously in [START_REF] Zheng | Second order sliding mode control for a quadrotor UAV[END_REF], which leads to additional dynamics

b Lj = U j L j j 2 (z, , ✓, ), (40)
where L j are the controls subjected by rotor dynamics.

Recall that for z dynamics, the second derivative of the position error, e z = z z des , has the form

ëz = g cos ✓ cos m (L z + d ze ) zdes ,
where the control U z is substituted by rotor dynamic output L z from (40).

The second derivative of the position error can be rewritten

ëz = L z z + z ,
where z = (cos ✓ cos )/m, z = g zdes (cos ✓ cos d ze )/m. The first derivative of the sliding surface is obtained

Ṡz = ëz + ↵ z ėz = L z z + z + ↵ z ėz .
Using (40), its second derivative is computed as

Sz = U z L z b z + ˙ z L z + ˙ z + ↵ z ( z L z + z ) = U z z b + dz , where dz = ˙ z + ↵ z ( z L z + z ) + (L z z )/b + ˙ z L z .
In Theorem 1 (see Appendix) it is stated that if the second order control is selected as In [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF] it is stated that if

U z = b z Dz (t, X) d Ṡz c 2 + S z | Ṡz | 2 + |S z | ,
ëx = Dx (t, X) d ėx c 2 + ↵ x e x | ėx | 2 + ↵ x |e y | , ëy = Dy (t, X) d ėy c 2 + ↵ y e y | ėy | 2 + ↵ y |e y | , then e i = ėi = 0, i = x, y is reached in a finite time provided that Dx (t, X) > |d xe + ẍdes |, Dy (t, X) > |d ye + ÿdes |.
x and y dynamics are not influenced directly by the rotors, hence their stability does not need the introduction of an auxiliary sliding surface. The respective controls for x, y positions are

U x = m Dx (t, X) U z d ėx c 2 + ↵ x e x | ėx | 2 + ↵ x |e x | , U y = m Dy (t, X) U z d ėy c 2 + ↵ y e y | ėy | 2 + ↵ y |e y | .
Such a design is admissible since the control U z is always separated with the zero and U x , U y are always continuous by definition of the quasi-continuous control.

The controls for other dynamics can be designed following similar computations as for z, so only final expressions are given for roll, pitch and yaw, respectively. The second derivative of the position error, e = des , taking into account (40), has the form

ë = qr I yy I zz I xx + 1 I xx (L + d ) ¨ des ,
The roll control is The pitch control is

U = b D (t, X) d Ṡ c 2 + S | Ṡ | 2 + |S | , where d = ˙ + ↵ ( L + ) (L )/b ˙ L , = 1
U ✓ = b ✓ D✓ (t, X) d Ṡ✓ c 2 + S ✓ | Ṡ✓ | 2 + |S ✓ | , where d✓ = ˙ ✓ + ↵ ✓ ( ✓ L ✓ + ✓ ) (L ✓ ✓ )/b, ✓ = 1/I yy , ✓ = pr(I zz I xx )/I yy + d /I yy ✓des , with D✓ (t, X) > | d✓ + ✓des |.
Recall that for dynamics, the second derivative of the position error, e = des , has the form

ë = pq I xx I yy I zz + 1 I zz (L + d ) ¨ des ,
The yaw control is

U = b D (t, X) d Ṡ c 2 + S | Ṡ | 2 + |S | , where d = ˙ + ↵ ( L + ) (L )/b, = 1/I zz , = pq(I xx I yy )/I zz + d /I zz ¨ des , with D (t, X) > | d + ¨ des |.
The proposed qc-SMC can be modified to counteract the chattering avoiding the saturation functions, and using the quasi-continuous function itself as an approximation of the sign on the plane, with a mild modification by adding a small constant % i > 0 in the denominator:

d Ṡi c 2 + S i % i + | Ṡi | 2 + |S i | , d ėi c 2 + ↵ i e i % i + | ėi | 2 + ↵ i |e i | .
where % i is strictly related with accuracy. The smaller is % i the higher is the effort on the rotors, which results in a more accentuated oscillation of the control, but with a smaller convergence error, and vice versa. According to [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF], a finite-time convergence of the system can be achieved in the ideal case, when % i = 0 and there is no measurement or implementation (digital) noises. In our case, since these restrictions are not satisfied, the convergence is assured with respect to a compact set around the desired trajectory. Then % i are tuned accordingly to achieve a trade-off between control oscillations and convergence error.

Concluding the previous preliminary study, the quasi-continuous control is suitable to be applied in quadrotor regulation taking into account rotor dynamics without modification of the sliding surfaces, which have been designed for the conventional SMC without [START_REF] Zheng | Second order sliding mode control for a quadrotor UAV[END_REF]. It means that the function sign can be approximated efficacely with quasi-continuous function in SMC for quadrotors, and its finite time stability is proven in [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF], considering % i = 0. Then, considering the imposed quadrotor dynamics as in 1-SMC design ( ėi + ↵ i e i = 0 for i = x, y, z, , ✓, ), the final expressions of qc-SMC chosen for implementation and comparison are given below.

The control for z is given by

U z = m cos ✓ cos ✓ g zdes + ↵ z ėz + ⇣ 1 2 ⌫(X) 2 + 2%(X) + ⌫(X) p ⌫ 2 (X) + 4%(X) + L z b L z ↵ z + $ z ⌘ d Ṡz c 2 + S z % z + | Ṡz | 2 + |S z | ◆ ,
where

$ z > | ˙ z + ˙ z L z |
is a tuning parameter. The controls for i = x, y are given by

U i = m U z ✓ (d xx + 2d yy + 2d zz + $ i ) d ėi c 2 + e i % i + | ėi | 2 + |e i | + ïdes ↵ i ėi ◆ , (41)
where $ i > 0 is a tuning parameter with Di (t, X) > |d ie + ïdes |.

The controls for roll, pitch and yaw respectively are given by

U =I xx ✓ qr I yy I zz I xx + ũ ↵ ė + ¨ des ◆ , U ✓ =I yy ✓ pr I zz I xx I yy + ũ✓ ↵ ✓ ė✓ + ✓des ◆ , U =I zz ✓ pq I xx I yy I zz + ũ ↵ ė + ¨ des ◆ ,
where the auxiliary controls are defined as

ũ = 1 I xx d Ṡ c 2 + S % + | Ṡ | 2 + |S | ⇣ K f 1 (X) + D 1 p |U z | + K f 2 (X) + D 2 + $ + ↵ L b L ⌘ , ũ✓ = 1 I yy d Ṡ✓ c 2 + S ✓ % ✓ + | Ṡ✓ | 2 + |S ✓ | ⇣ K✓ f ✓1 (X) + D ✓1 p |U z | + K✓ f ✓2 (X) + D ✓2 + $ ✓ + ↵ ✓ ✓ L ✓ ✓ b L ✓ ⌘ , ũ = 1 I zz d Ṡ c 2 + S % + | Ṡ | 2 + |S | ⇣ K f 1 (X) + D 1 p |U z | + K f 2 (X) + D 2 + $ + ↵ L b L ⌘ , with $ i > | ˙ i |, i = , ✓,
tuning parameters. Homogeneous differentiator [START_REF] Perruquetti | Finite-Time Observers: Application to Secure Communication[END_REF] (see the expression (42) below) is used to estimate the value of the first derivative of the sliding surfaces Ṡi .

Physical constraints and tools

Since UAVs are physical objects, then they are subjected to physical constraints, which correspond to their maximal linear and angular velocities and accelerations. Consequently, a smooth position trajectory is desired.

In literature many articles are presented to generate desired trajectories admissible to the system dynamics. Since the focus of the presented article is the control of a mini drone, then simple method is used which can take into account UAV physical constraints imposed during indoor experiments. The desired step references for coordinates x, y, z are filtered (using a third order filter) to obtain realistic velocity and acceleration trajectories

⇠ f ⇠ wp = 1 (1 + ⌧ pos s) 3 ,
where ⇠ wp is the way-point step reference, ⇠ f = (x des , y des , z des ) is the filtered signal, ⌧ pos is the weight to create the desired smooth trajectory. The physical constraints are chosen where T max is 90% of the maximal thrust achievable by the UAV to ensure that we still have some available thrust to proper control the drone in case of safe maneuvers; T min is the minimal thrust. The reference angles des , ✓ des are saturated in an admissible range of min , ✓ min minimal values and max , ✓ max maximal values, to avoid the singularities.

v min =
The desired angles des , ✓ des are derived from the expressions for U x and U y , using the desired value of

U ,✓, min = 2 4 l⇢AR 2 C T stat sin(✏)(! 2 max ! 2 min ) l⇢AR 2 C T stat cos(✏)(! 2 max ! 2 min ) 2⇢AR 3 C Q (! 2 max ! 2 min ) 3 5 , U ,✓, max = 2 4 l⇢AR 2 C T stat sin(✏)(! 2 max ! 2 min ) l⇢AR 2 C T stat cos(✏)(! 2 max ! 2 min ) 2⇢AR 3 C Q (! 2 max ! 2 min ) 3 5
.

In order to implement the designed controls the first and the second derivatives of the desired angles have to be calculated, which are obtained using a homogeneous finite time differentiator, as in [START_REF] Perruquetti | Finite-Time Observers: Application to Secure Communication[END_REF]

ż1 (t) = 1 dz 1 (t) f (t)c 1 + z 2 (t), ż2 (t) = 2 dz 1 (t) f (t)c 2 + z 3 (t), ż3 (t) = 3 dz 1 (t) f (t)c 3 , (42)
where z 2 R 3 is the differentiator state and f (t) is the signal to be differentiated, i = 1 j⌧ for j = 1, 2, 3 and any ⌧ 2 ( 1 3 , 0), while the parameters 1 , 2 , 3 are selected to ensure the Hurwitz property of the matrix

A = 2 4 1 1 0 2 0 1 3 0 0 3 5 .
Then the system (42) performs estimation of the first and second derivative of f (t) in a finite time: f (t) = z 2 (t), f (t) = z 3 (t). Increasing the smallest eigenvalue of A improves the rate of convergence. In this work ⌧ < 0 is sufficiently big, the eigenvalues of A are chosen by trial and error accordingly to the desired dynamics and in such a way to avoid undesired response delay.

Control design summary

The generic scheme of hierarchical control algorithm presented above is given in Fig. 3. Homogeneous differentiator is used to estimate the first and second derivatives of angles and the first derivative of the sliding surfaces. The stability of the system can always be proved respecting the maximal value of , ✓, D x , D y , D z . In the proposed algorithms, the control gains are functions of the state and wind velocity, then using a wind estimator the wind velocity estimates can be substituted on-line in the control law, as in Fig. 4. In this way, by adapting the control amplitude, the regulator effort on the rotors can be reduced when it is possible. A good wind estimator is desired to reduce the uncertainty that will affect directly the performance of the proposed control. However, in this paper maximal admissible wind values are considered as input to the control.

Numerical Validation

Numerical experiments are performed in the simulator using UAV Simulink library, built with nonlinear aerodynamic coefficient, and using the rotor model in eq. ( 4). The quadrotor model was developed using blade momentum theory in aerodynamic science and next validated through indoor experiments with the commercial Parrot Ar Drone 2.0 at low/moderate velocity at the ONERA Aerodynamics Aeroelasticity Acoustics Department. Wind perturbations, affecting the UAV nonlinear model, are simulated as sinusoids, which is a reasonable assumption since the gust generator at ONERA lab allows to create a sinusoidal wind gust profile, hence the sinusoid wind effect over time is correlated to the motion of the drone through the turbulent area. For simplicity of demonstration the wind velocities are selected as constant input in qc-SMC and in 1-SMC the wind velocities are replaced with their constant upper bounds, see Fig. 5 (wind x, wind y, wind z). This restriction is imposed in the article to better present the new design strategy. Conventional SMC (conv-SMC), used for comparisons, is built as ũi = C i signS i for constant C i kd i k 1 , i = x, y, z, , ✓, and suitable under a popular hypothesis that perturbations in the basic quadrotor model are bounded functions of time. In conv-SMC, as usual for implementation, the function sign is replaced by the approximation sat ⇠ i in eq. (39). Constant gains C i are tuned to compensate the disturbances (dist x, dist y, dist z, dist , dist ✓, dist ) in Fig. 5 and taking their maximal values. Since three nonlinear controllers are studied, a rule is imposed to make the comparisons: the same performance of the controllers are imposed by tuning and keeping the same ⇠ i for conv-SMC and 1-SMC, % i are adjusted to have a similar performance for the qc-SMC, and same values of ↵ i are also considered. Unfortunately, 1-SMC and conv-SMC cannot achieve the same accuracy as the qc-SMC since they become unstable due to the high control oscillations if the convergence error is very small. For this reason, a second try is done keeping the same ⇠ i for 1-SMC and conv-SMC and "relaxing" them (changing ⇠ i to have bigger position convergence error), while a more accuracy is achieved for qc-SMC. Hence, we consider the values listed in Tables 1,3, 2, 4 where żmax , z max are with negative signs because z axes is taken with the positive direction towards down as in Fig. 1. For simulations, a generic flight case is illustrated which couples together forward, lateral and vertical flights, moving the UAV in x, y, z positions. The corresponding external disturbances, influenced by the proposed varying gain functions, and wind are illustrated in Fig. 5, where wind in z is chosen small because of the small UAV mass and limited thrust. Figures 6,7, 8, 9, 10, 11 illustrate the contribution of the controls together with state variables and reference states. After the presentation of the working hypothesis, the following conclusions are identified after these numeric experiments:

First, the simulation clearly shows that using a varying-gain function instead of a constant gain improves the transients, because the SMCs using this function are subjected to less effort on the rotors. Under the same ⇠ i , conv-SMC is slightly more robust than the 1-SMC due to the fact that in the tuned constant gain of conv-SMC some inner dynamics are automatically considered (derivative of the control U z , derivative of the wind velocity, etc. . . ) but at the price of much higher control oscillation. This fact is due to the simplifications adopted in the disturbances design in sec. 3.1. However, more accuracy can be achieved also for the 1-SMC, changing the parameter ⇠ i . While the tuned gain for the conv-SMC has to be done "a posteriori", the varying-gain of 1-SMC and of qc-SMC can be tuned "a priori" making an hypothesis on maximal wind velocities. This first point shows the effectiveness of the first contribution: tuning the varying-gains influenced by the state, the control itself, and the wind velocity.

Second, as already discussed, 1-SMC and conv-SMC tend to become unstable if higher accuracy (changing ⇠ i ) is imposed because the rotors dynamics is not considered in their designs. qc-SMC can achieve better accuracy among all the controls ensuring an acceptable control oscillation which is similar to the control in 1-SMC, but the resulted gain tuning is far to be optimal. Then, the trajectory can be reached changing freely % i and enforcing the control to be more robust as in Figures 12,[START_REF] Besnard | Quadrotor vehicle control via sliding mode controller driven by sliding mode disturbance observer[END_REF] without compromising the stability. This second point shows the effectiveness of the second contribution: the importance of including rotors dynamics inside the control design.

Last, regarding the qc-SMC, the error convergence of the system is affected by the introduction of % i in the denominator of quasi-continuous function. Reducing the convergence error, increases the control oscillations, hence a good trade-off between effectiveness of the control and the effort on the rotors has to be found. However, thanks to the varying-gains built in subsection 3.2, a wind estimator can be coupled with the control, then wind velocity can be substituted on-line. In this way, by adapting the control amplitude, the regulator effort on the rotors can be reduced when it is possible, allowing to achieve a better performance with less control oscillations. This last point highlights the importance of the first contribution: tuning the varying-gains influenced by the wind velocity.

As a result, the presented simulations demonstrate the stability of the system affected by wind velocity and realistic external disturbances, respecting the limitations of the mini drone mass and thrust.

Conclusion

In this work, a robust sliding mode control approach is used to stabilize a small quadrotor UAV under wind perturbations. It is considered that the disturbance bound for a UAV quadrotor at a low/moderate velocity depends on the control itself, the wind speed and the state of the UAV. To this end, In this work constant maximal admissible wind values are considered as direct input to the control. However, since this approach is strictly related to the wind estimation, the coupling with the estimator, as in Fig. 4 is desired. Further works lead to a deeper study of their interaction, which must ensure the stability of the system in a compact set around the equilibrium depending by the estimation incertitude. functions are considered in eq. ( 17) which leads to the validity of the Theorem 4. where ↵ > 0 and > 0 are tuning parameters. In our case = 1 and g, h are respectively the computed lower and the upper-bounds. Note that the control (A.3) is continuous everywhere outside of the origin.

Theorem 1. Consider the system in (A.1) and assume that the restrictions (A.2) are satisfied, then there exist ↵ > 0 sufficiently big and > 0 such that the control (A.3) makes the system globally finite-time convergent.

Figure 1 :

 1 Figure 1: UAV Parrot Drone 2.0 and its configuration

Figure 2 :

 2 Figure 2: Rotor model

  ), and performing calculations with the indicated simplifications, we have the disturbance d = [d x , d y , d z , d , d ✓ , d ] represented by the following terms

where d•c 2 =

 2 | • | 2 sign(•), then the point S z = Ṡz = 0 is reached in a finite time provided that Dz (t, X) > | dz |. The position errors, e x = x x des , e y = y y des , have the dynamics ëx = U x U z m d xe ẍdes , ëy = U y U z m d ye ÿdes .

  /I xx , = qr(I yy I zz )/I xx + d /I xx ¨ des , with D (t, X) > | d + ¨ des |. The second derivative of the position error, e ✓ = ✓ ✓ des , taking into account (40), has the form ë✓ = pr I zz I xx I yy + 1 I yy (L ✓ + d ✓ ) ✓des ,

  min m sin ✓ min ÿmin = T min m sin min zmin = T min m

  des = arcsin (U x sin des U y cos des ) ,✓ des = arcsin ✓ U x cos des + U y sin des cos des ◆ .U x and U y are illustrated in equations (32), (33) with the saturation (39) for 1-SMC, and in the eq. (41) for qc-SMC. Physical constraints for angles and controls are chosen to be

Figure 3 :

 3 Figure 3: Hierachical control scheme

Figure 4 :

 4 Figure 4: Control-wind estimator coupling. The red arrow indicates how varying wind estimation can be merged with the control.

Figure 6 :

 6 Figure 6: x, y, z position and respective controls (-: desired trajectory; -: 1-SMC; -: conv-SMC; -: qc-SMC).

Figure 7 :

 7 Figure 7: , ✓, angles and respective controls (-: 1-SMC; -: conv-SMC; -: qc-SMC).

Figure 8 :

 8 Figure 8: Linear and angular velocities (-: 1-SMC; -: conv-SMC; -: qc-SMC).

Figure 9 :

 9 Figure 9: Rotors angular velocities for conventional SMC.

Figure 10 :

 10 Figure 10: Rotors angular velocities for 1-SMC.

Figure 11 :

 11 Figure 11: Rotors angular velocities for qc-SMC.

Figure 12 :

 12 Figure 12: x, y, z position and respective controls (-: qc-SMC changing % i ; -: desired trajectory).

Figure 13 :

 13 Figure 13: , ✓, angles and respective controls for qc-SMC changing % i .

Table 1 :

 1 Parrot drone parameters

	Par.	R	l	h	g	✓ 0	m	I xx
	Val.	0.1	0.185	0.0125 9.81	23.9	0.472 0.00356
	Unit	m	m	m	m/s 2	deg	Kg	Kg.m 2
	Par.	I yy	I zz	⇢		a	C D0	stat
	Val. 0.00402 0.00712	1.25	0.111 4.6542 2.15	0.1056
	Unit Kg.m 2 Kg.m 2 Kg/m 3
	Par. C T stat	K D	K z
	Val. 0.0223	0.06	0.09
	Unit

Table 2 :

 2 Controls saturationU z (N ) U (Nm) U ✓ (Nm) U (Nm) min.

		1.4	0.14	0.14	0.30
	max.	5.6	0.14	0.14	0.30

Table 3 :

 3 Saturation parametersPar.

			max	✓ max	max		min	✓ min	min	˙	max	✓max
	Val.	45	45	180		45	45	180		40	40
	Unit	deg	deg	deg	deg	deg	deg deg /s deg /s
	Par.	˙	max	˙	min	✓min	˙	min	x max y max	z max	x min
	Val.	10		40	40		10	100	100		100	0
	Unit deg /s deg /s deg /s deg /s	m	m		m	m
	Par.	y min	z min	ẋmax	ẏmax	żmax ẋmin	ẏmin	żmin
	Val.		0		0	5		5	1	5		5	1
	Unit		m		m	m/s	m/s	m/s m/s	m/s	m/s
	Par. ! max	! min
	Val.	400	200
	Unit rad/s rad/s

stat

), K 2 = ⇢ARhK D , K 3 = lK z ⇢AR, K = K, f 1 (X) = K 1 max j |u j |+K 2 max j |v j |+K 3 max j |w j s j |, D 1 = K 1 D x +

Acknowledgment

The authors want to acknowledge Hauts-de-France region and ONERA to have financially supported this work.

the constant gain, which is proportional to the sign functions in conventional sliding mode controls, is replaced with a time-varying function. In this way we aim to reduce as much as possible the control effort on the rotors, adapting in real-time the weighted function used in the control design.

The issue originated by introduction of the rotor transfer functions is studied and several solutions for decreasing chattering are applied. Then, the contribution with respect to conventional first order sliding mode control is highlighted, and a recent quasi-continuous high order sliding mode control is also tested and discussed in simulations.

The presented results of numerical experiments confirm the effectiveness of the designed control to stabilize the drone under varying wind, and its stability is proved for a known maximal admissible wind value.

Appendix A. 1

Following [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF] consider a double integrator system:

where X(t) = [x(t), ẋ(t)] T 2 R 2 is the state vector, u(t) 2 R is the control input, two functions g : R 3 ! R and h : R 3 ! R ensure forward existence and uniqueness of the system solutions at least locally. In addition, there are two known functions g : R 3 ! R and h : R 3 ! R such that for all X 2 R 2 and t 0 g(t, X) g(t, X) > 0, |h(t, X)|  h(t, X).

(A.

2)

The following control for quasi-continuous SMC can be proposed for (A.1) (a more generic case is studied in [START_REF] Ding | Simple homogeneous sliding-mode controller[END_REF] where bounded time-varying positive


François Dubois

Equivalent partial differential equations of a lattice Boltzmann scheme

Keywords:

published or not. The documents may come

Equivalent partial differential equations of a lattice

Boltzmann scheme

1) Discrete geometry

• We denote by d the dimension of space and by L a regular d-dimensional lattice. Such a lattice is composed by a set L 0 of nodes or vertices and a set L 1 of links or edges between two vertices. From a practical point of view, given a vertex x, there exists a set V (x) of neighbouring nodes, including the node x itself. We consider here that the lattice L is parametrized by a space step ∆x > 0. For the fundamental example called D2Q9 (see e.g. [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF], the set V (x) is given with the help of the family of vectors (e j ) 0≤j≤J defined by J = 8, (1.1)

e j = 0 0 , 1 0 , 0 1 , -1 0 , 0 -1 , 1 1 , -1 1 , -1 -1 , 1
-1 and the vicinity (1.2)

V (x) = { x + ∆x e j , 0 ≤ j ≤ J } .

• In the general case, we still suppose that the equation (1.2) holds but we do not make any precise definition concerning the integer J and the nondimensionalized vectors (e j ) 0≤j≤J .

Nevertheless if x is a node of the lattice (x ∈ L 0 ), then y j = x + ∆x e j is an other node of the lattice, i.e. y j ∈ L 0 .

2) Lattice Boltzmann framework

• We introduce a time step ∆t > 0 and we suppose that the celerity λ defined according to

(2.1) λ = ∆x ∆t remains fixed. Then we introduce a local velocity v j in such a way that (2.2) ∆t v j = ∆x e j , 0 ≤ j ≤ J .

In this d-dimensional framework we will denote by v α j (1 ≤ α ≤ d) the Cartesian components of velocities v j . Recall that if x is a node of the lattice, the point x + ∆t v j is also a node of the lattice:

(2.3)

x ∈ L 0 =⇒ x + ∆t v j ∈ L 0 , ∀j = 0, . . . J.

• According to D' Humières (1992), the lattice Boltzmann scheme describes the dynamics of the density f j (x, t) of particles of velocity v j at the node x and for the discrete time t. We introduce the d + 1 scalar "conservative variables" W (x, t) composed by the density ρ and the momentum q. Note that it is also possible to take into account the conservation of the total energy (see D'Humières's article for example). We have

(2.4) ρ(x, t) = J j=0 f j (x, t) ≡ W 0 (x, t) (2.5) q α (x, t) = J j=0 v α j f j (x, t) ≡ W α (x, t) , 1 ≤ α ≤ d ,

and

(2.6)

W (x, t) = ρ(x, t), q 1 (x, t), • • • , q d (x, t) .
When a state W is given in space R d+1 , a Gaussian (or any other choice) equilibrium distribution of particles is defined according to (2.7)

f j eq = G j (W ) , 0 ≤ j ≤ J in such a way that (2.8) J j=0 G j (W ) ≡ W 0 , J j=0 v α j G j (W ) ≡ W α , 1 ≤ α ≤ d .
• Following D 'Humières (1992), we introduce the "moment vector" m according to (2.9)

m k = J j=0 M k j f j , 0 ≤ k ≤ J .
For 0 ≤ i ≤ d, the moments m i are identical to the conservative variables:



Equivalent partial differential equations of a lattice Boltzmann scheme

(2.10)

m 0 ≡ ρ , m α ≡ q α , 1 ≤ α ≤ d .
In other words, the matrix M satisfies (2.11)

M 0 j ≡ 1 , M α j ≡ v α j , 0 ≤ j ≤ J , 1 ≤ α ≤ d .
We assume that vectors (e j ) 0≤j≤J are chosen such that the (d + 1) × (J + 1) matrix (M k j ) 0≤k≤d, 0≤j≤J is of full rank. With this hypothesis, the conservative moments W introduced in relation (2.6) are independent variables.

• When a particle distribution f is given, the moments are evaluated according to (2.9). The matrix M is supposed to be invertible and the inverse relation takes the form:

(2.12)

f j = J k=0 (M -1 ) j k m k , 0 ≤ j ≤ J .
When f j eq is determined according to the relation (2.7), the associated equilibrium moments m k eq are given simply according to (2.9), i.e. in this case

(2.13)

m k eq = J j=0 M k j f j eq , 0 ≤ k ≤ J .
We remark also that by construction (relation (2.8)), we have

(2.14) m i eq = m i = W i , 0 ≤ i ≤ d .

3) Collision step

• The collision step is local in space and is naturally defined in the space of moments. If m k (x, t) denotes the value of the k th component of the moment vector m at position x and time t, the same component m k * (x, t) of the moment after the collision is trivial by construction for the conservative variables:

(3.1) m i * (x, t) = m i (x, t) , 0 ≤ i ≤ d .
For the non-conservative components of the moment vector, we fix the ratio s k (k ≥ d + 1) between the time step ∆t and the relaxation time τ k of an underlying process:

(3.2) s k = ∆t τ k , d + 1 ≤ k ≤ J .
• Then m k * (x, t) after the collision is defined according to

(3.3) m k * (x, t) = (1 -s k ) m k (x, t) + s k m k eq , d + 1 ≤ k ≤ J .
Proposition 1. Explicit Euler scheme. The numerical scheme (3.3) is exactly the explicit Euler scheme relative to the continuous in time relaxation equation

(3.4) d dt m k -m k eq + 1 τ k m k -m k eq = 0 , d + 1 ≤ k ≤ J .

François Dubois

Proof of Proposition 1. Following e.g. [START_REF] Strang | An introduction to applied mathematics[END_REF], we know that the explicit Euler scheme for the evolution (3.4) takes the form

(3.5) 1 ∆t m k -m k eq (t + ∆t) -m k -m k eq (t) + 1 τ k m k -m k eq (t) = 0 .
We have by construction the relation (3.1), that is m i (t + ∆t) = m i (t) for 0 ≤ i ≤ d with these notations. Then W (t+∆t) = W (t) and, due to the relation (2.7), f j eq (t+∆t) = f j eq (t) after the collision step for all the components j of the particle distribution. Due to (2.13), we deduce that m k eq (t + ∆t) = m k eq (t) for all k ≤ J. Thus the expression (3.5) takes the simpler form

(3.6) 1 ∆t m k (t + ∆t) -m k (t) + 1 τ k m k -m k eq (t) = 0 , which is exactly (3.3), except the change of notations: m k (t + ∆t) is replaced by m k * .
• We remark also that the classical stability condition for the explicit Euler scheme (see again e.g. the book of Strang) takes the form

(3.7) 0 ≤ ∆t ≤ 2 τ k .
We will suppose in the following that

(3.8) 0 < s k ≤ 2 , d + 1 ≤ k ≤ J .
to put in evidence that the moments m k are not conserved for index k greater than d + 1.

We remark also that for the physically relevant Boltzmann equation, the relaxation times τ k have a physical sense. With the lattice Boltzmann scheme itself, these physical constants are no longer correctly approximated whereas the ratios s k = ∆t τ k are supposed to be fixed in all what follows. Despite the usual "LBE" denomination, a lattice Boltzmann scheme is not a numerical method to approach the Boltzmann equation !

• The particle distribution f j * after the collision step follows the relation (2.12). We have precisely after the collision step (3.9)

f j * = J k=0 (M -1 ) j k m k * , 0 ≤ j ≤ J .

4) Advection step

• The avection step of the lattice Boltzmann scheme claims that after the collision step, the particles having velocity v j at position x go in one time step ∆t to the j th neighbouring vertex. Thus the particle density f j (x + v j ∆t, t + ∆t) at the new time step in the neighbouring vertex is equal to the previous particle density f j * (x, t) at the position x after the collision: (4.1) f j (x + v j ∆t, t + ∆t) = f j * (x, t) . We re-write this relation in term of the "arrival" node x + v j ∆t. We set x = x + v j ∆t, then we have x = x -v j ∆t and going back to the notation x, we write the relation (4.1) in the equivalent manner  Equivalent partial differential equations of a lattice Boltzmann scheme

(4.2) f j (x, t + ∆t) = f j * (x -v j ∆t, t) , 0 ≤ j ≤ J , x ∈ L 0 .
Proposition 2. Upwind scheme for the advection equation. The scheme (4.2) for the advection step of the lattice Boltzmann method is nothing else that the explicit upwind scheme for the advection equation

(4.3) ∂f j ∂t + v j • ∇f j = 0 , 0 ≤ j ≤ J ,
with a so-called Courant-Friedrichs-Lewy number σ j in the j th direction of the lattice defined by (4.4)

σ j ≡ | v j | ∆t ∆x | e j |
equal, due to the definition (2.2), to unity: σ j = 1.

Proof of Proposition 2.

When the Courant-Friedrichs-Lewy number σ j is equal to unity, it is classical (see e.g. [START_REF] Strang | An introduction to applied mathematics[END_REF]) that the upwind scheme is exact for the advection equation.

5) Equivalent equation at zero order

• The lattice Boltzmann scheme is defined by the relations (2.4) to (2.9), (3.3) and (4.2). It is parametrized by the lattice step ∆x, the matrix M linking the particle distribution f and the moment vector m, the choice of the conservative moments, the nonlinear equilibrium function G(•), the time step ∆t and the ratios s k between the time step and the collision time constants for nonequilibrium moments. In what follows, we fix the geometrical and topological structure of the lattice L, we fix the matrix M and the equilibrium function G(•), we fix also the ratio λ defined in (2.1) and last but not least, we suppose that the parameters s k for k ≥ d + 1 have a fixed value. Then the whole lattice Boltzmann scheme depends on the single parameter ∆t.

• We explore now formally what are the partial differential equations associated with the Boltzmann numerical scheme, following the so-called "equivalent equation method" introduced and developed by [START_REF] Lerat | Noncentered Schemes and Shock Propagation Problems[END_REF] and [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF]. This approach is based on the assumption, that a sufficiently smooth function exists which satisfies the difference equation at the grid points. This assumption gives formal responses to put in evidence partial differential equations that minimimize the truncation errors of the numerical scheme. Nevertheless, we note here that this method of analysis fails to predict initial layers and boundary effects properly, as discussed by [START_REF] Griffiths | On the scope of the method of modified equations[END_REF] or [START_REF] Chang | A critical analysis of the modified equation technique of Warming and Hyett[END_REF]. The idea of the calculus is to suppose that all the data are sufficiently regular and to expand all the variables with the Taylor formula.



François Dubois

Proposition 3. Taylor expansion at zero order. With the lattice Boltzmann defined previously, we have (5.1) f j (x, t) = f j eq (x, t) + O(∆t) , 0 ≤ j ≤ J , (5.2) f j * (x, t) = f j eq (x, t) + O(∆t) , 0 ≤ j ≤ J , with f j eq defined from the conservative variables W according to the relation (2.7).

Proof of Proposition 3.

The key point is to expand the relation (4.2) relative to the infinitesimal ∆t. We have on one hand (5.3)

f j (x, t + ∆t) = f j (x, t) + O(∆t)
and on the other hand (5.4)

f j * (x -v j ∆t, t) = f j * (x, t) + O(∆t) Then m k * (x, t) = J j=0 M k j f j * (x, t) = m k (x, t) + O(∆t) and
(5.5)

m k * (x, t) -m k (x, t) = O(∆t) . But, due to (3.3), we have (5.6) m k * (x, t) -m k (x, t) = -s k m k (x, t) -m k eq (x, t
) . From (5.5) and (5.6) we deduce, due to the fact that s k = 0 when k ≥ d + 1 :

(5.7) m k (x, t) = m k eq (x, t) + O(∆t) , k ≥ d + 1 . We insert (5.7) into (5.5) and we deduce (5.8) m k * (x, t) = m k eq (x, t) + O(∆t) , k ≥ d + 1 . Taking into account the relations (2.14) and (3.1) on one hand and (2.12) and (3.9) on the other hand, we deduce (5.1) and (5.2) from (5.7) and (5.8).

6) Taylor expansion at first order

• We expand now the relation (4.2) one step further with respect to the time step ∆t. We introduce the second order moment (6.1)

F α β ≡ J j=0 v α j v β j f j eq , 1 ≤ α, β ≤ d .
We denote in the following ∂ t instead of ∂ ∂t and ∂ β in place of ∂ ∂x β . Then we have the following result at the first order. Proposition 4. Euler equations of gas dynamics. With the lattice Boltzmann scheme previously defined, we have the conservation of mass and momentum at the first order: (6.2)

∂ t ρ + d β=1 ∂ β q β = O(∆t) 
Equivalent partial differential equations of a lattice Boltzmann scheme

(6.3) ∂ t q α + d β=1 ∂ β F α β = O(∆t) .
Proof of Proposition 4.

We expand both sides of relation (4.2) up to first order:

f j (x , t + ∆t) = f j (x , t) + ∆t ∂ t f j + O(∆t 2 ) f j * (x -v j ∆t , t) = f j * (x , t) -∆t v β j ∂ β f j * + O(∆t 2
) . We take the moment of order k of this identity:

m k (x , t) + ∆t ∂ t m k + O(∆t 2 ) = m k * (x , t) -∆t J j=0 M k j v β j ∂ β f j * + O(∆t 2 )
and we use the previous Taylor expansions (5.1) (5.2) at the order zero:

(6.4) m k (x , t) + ∆t ∂ t m k eq = m k * (x t) -∆t J j=0 M k j v β j ∂ β f j eq + O(∆t 2 ) .
We take k = 0 inside the relation (6.4). We get (6.2) since m 0 (x , t) ≡ m 0 * (x , t) ≡ ρ(x , t). Considering now the particular case k = α with 1 ≤ α ≤ d, we have also m α (x , t) ≡ m α * (x , t) ≡ q α (x , t) and the relation ( 6.3) is a direct consequence of the definition (6.1) and the property (2.11).

Proposition 5. Technical lemma. We introduce the "conservation defect" θ k according to the relation

(6.5) θ k (x , t) = ∂ t m k eq + J j=0 M k j v β j ∂ β f j eq ≡ J j=0 M k j ( ∂ t f j eq + v β j ∂ β f j eq ) .
Then we have the following properties:

(6.6) m k (x , t) = m k eq (x , t) - ∆t s k θ k + O(∆t 2 ) , k ≥ d + 1 , (6.7) m k * (x , t) = m k eq (x , t) - 1 s k -1 ∆t θ k + O(∆t 2 ) , k ≥ d + 1 , (6.8) ∂ β f j * = ∂ β f j eq -∆t J k=d+1 1 s k -1 (M -1 ) j k ∂ β θ k + O(∆t 2 ) .
Proof of Proposition 5.

We start from the relation (6.4) and we have observed at the previous proposition that (6.9)

θ i = O(∆t) , 0 ≤ i ≤ d .
We remark also that from the relation (5.6), we have

m k (x , t) -m k eq (x , t) = 1 s k m k (x , t) -m k * (x , t) if k ≥ d + 1.
Then the relation (6.6) is a direct consequence of (6.4) and the definition (6.5). In consequence, the relation (6.7) follows from (6.6) and (6.4). Due to (6.7), (6.9) and (3.9), we have (6.10)

f j * (x , t) = f j eq (x , t) -∆t k≥d+1 1 s k -1 (M -1 ) j k θ k + O(∆t 2 )  François Dubois
and the relation (6.8) follows from derivating (6.10) in the direction x β .

7) Equivalent equation at second order

• We introduce the tensor Λ α β k according to

(7.1) Λ α β k ≡ J j=0 v α j v β j (M -1 ) j k , 1 ≤ α, β ≤ d , 0 ≤ k ≤ J .
We can now establish the major result of our contribution.

Proposition 6. Navier-Stokes equations of gas dynamics.

With the lattice Boltzmann method defined in previous sections and the conservation defect θ k defined in (6.5), we have the following expansions up to second order accuracy:

(7.2) ∂ t ρ + d β=1 ∂ β q β = O(∆t 2 ) (7.3) ∂ t q α + d β=1 ∂ β F α β -∆t k≥d+1 1 s k - 1 2 Λ α β k θ k = O(∆t 2 ) .
• A consequence of relation ( 7.3) is the fact that a lattice Boltzmann scheme approximates at second order of accuracy a Navier-Stokes type equation with viscosities µ k of the form (7.4)

µ k = ∆t 1 s k - 1 2 .
We refer for the details to D. D 'Humières (1992), [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] or to our recent survey (2007). The relations (7.4) are known as the "D'Humières relations". We observe that in practice, the scalar µ k is imposed by the physics and by the parameter ∆t is constrained by the space discretization ∆x and the relation (2.1). Then the parameter s k must be chosen in order to satisfy the D 'Humières relations (7.4).

Proof of Proposition 6. We start again from the identity (4.2). We expand both terms up to second order accuracy:

f j (x , t + ∆t) = f j (x , t) + ∆t ∂ t f j + 1 2 ∆t 2 ∂ 2 tt f j + O(∆t 3 ) f j * (x -v j ∆t , t) = f j * (x , t) -∆t v β j ∂ β f j * + 1 2 ∆t 2 v β j v γ j ∂ 2 βγ f j * + O(∆t 3
) . We take the moment of order i (0 ≤ i ≤ d) of this identity. We obtain:

(7.5)          m i (x , t) + ∆t ∂ t m i + 1 2 ∆t 2 ∂ 2 tt m i + O(∆t 3 ) = m i * (x , t) + -∆t J j=0 M i j v β j ∂ β f j * + 1 2 ∆t 2 J j=0 M i j v β j v γ j ∂ 2 βγ f j * + O(∆t 3 ) .
We use the microscopic conservation m i * (x , t) ≡ m i (x , t) in (7.5) and the previous Taylor expansion at order one, in particular the relation (6.8). We divide by ∆t and we deduce:  Equivalent partial differential equations of a lattice Boltzmann scheme

∂ t m i + 1 2 ∆t ∂ 2 tt m i = - J j=0 M i j v β j ∂ β f j eq + + ∆t J j=0 k≥d+1 M i j v β j 1 s k -1 (M -1 ) j k ∂ β θ k + + 1 2 ∆t J j=0 M i j v β j v γ j ∂ 2 βγ f j eq + O(∆t 2 ) . Then (7.6)                        ∂ t m i + d β=1 J j=0 M i j v β j ∂ β f j eq = = ∆t d β=1 J j=0 k≥d+1 M i j v β j 1 s k -1 (M -1 ) j k ∂ β θ k + + ∆t 2 -∂ 2 tt m i + d β=1 J j=0 M i j v β j v γ j ∂ 2 βγ f j eq + O(∆t 2 ) .
• We set i = 0 in the relation (7.6) and we look for the conservation of mass. Due to the property M 0 j ≡ 1, the sum over j in the second line of (7.6) is null since

J j=0 v β j (M -1 ) j k
is equal to zero. We have also the following algebraic calculus:

∂ 2 tt m 0 = ∂ 2 tt ρ = - d β=1 ∂ 2 tβ q β + O(∆t) = - d β=1 ∂ β ∂ t q β + O(∆t) = = d β=1 d γ=1 ∂ 2 βγ F β γ + O(∆t) = d β=1 d γ=1 J j=0 v β j v γ j ∂ 2 βγ f j eq + O(∆t)
and the third line of (7.6) is null up to second order accuracy. Thus the conservation of mass (7.2) up to second order accuracy is established.

• We set i = α with 1 ≤ α ≤ d and we look for the conservation of momentum. In this particular case, the relation (7.6) takes the form:

(7.7)                        ∂ t q α + d β=1 J j=0 v α j v β j ∂ β f j eq = = ∆t k≥d+1 1 s k -1 d β=1 J j=0 v α j v β j (M -1 ) j k ∂ β θ k + + ∆t 2 -∂ 2 tt q α + d β=1 J j=0 v α j v β j v γ j ∂ 2 βγ f j eq + O(∆t 2 ) .
We have now to play with some algebra:

 François Dubois -∂ 2 tt q α + d β=1 J j=0 v α j v β j v γ j ∂ 2 βγ f j eq = = d β=1 ∂ t ∂ β F α β + J j=0 v α j v β j v γ j ∂ 2 βγ f j eq + O(∆t) = d β=1 ∂ β J j=0 v α j v β j ∂ t f j eq + v γ j ∂ γ f j eq + O(∆t) = d β=1 ∂ β J j=0 v α j v β j J k=0 (M -1 ) j k θ k + O(∆t) = d β=1 ∂ β k≥d+1 J j=0 v α j v β j (M -1 ) j k θ k + O(∆t) = d β=1 ∂ β k≥d+1 Λ α β k θ k + O(∆t)
due to the definition (7.1). We deduce from (6.1), (7.7) and the above calculus: and the relation (7.3) is established.

∂ t q α + d β=1 ∂ β F α β = ∆t

8) Equivalent equation at second order

• The previous propositions establish that the equivalent partial differential equations of a Boltzmann scheme are given up to second order accuracy by the same result as the formal Chapman-Enskog expansion. We find Euler type equation at the first order (Proposition 4) and Navier-Stokes type equation at the second order (Proposition 6). Note that with the above framework no a priori formal two-time multiple scaling is necessary to establish the Navier-Stokes equations from a lattice Boltzmann scheme, as done previously in the contribution of D'Humières. We remark also that a so-called diffusive scaling like ∆t ∆x 2 = constant, instead of our condition (2.1) ∆t ∆x = constant, leads to the incompressible Navier-Stokes equations, as proposed by [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF]. In both cases, we have just to use the Taylor formula for a single infinitesimal parameter.

∂

  β θ k + O(∆t 2 ) .

9) Acknowledgments

The author thanks Li-Shi Luo to his kind invitation to present the scientific work of Orsay's team at ICMMES Conference in July 2005. The author thanks also with a great emphasis Pierre Lallemand for very helpfull discussions all along the elaboration of this contribution. Last but not least, the referees transmitted to the author very good remarks that have been incorporated inside the present edition of the article.