Van Cuong Nguyen

Gilney Damm
email: gilney.damm@univ-eiffel.fr

MIMO Conditional Integrator Control for Unmanned Airlaunch

Keywords: Conditional Integrator, Airlaunch, F-16 Aircraft, Nonlinear System

A Nonlinear Multiple Input Multiple Output (MIMO) controller based on the conditional integrator technique is designed for the robust stabilization of a new satellite launching strategy called (unmanned) airlaunch. This control technique performs as a nonlinear robust controller when the system is outside a boundary layer defined by the controller (large errors), and as a nonlinear controller that provides an integral term when the system enters into the boundary layer (small errors). While the airlaunch strategy consists in using a two-stages launching system. The first stage is composed of an airplane (manned or unmanned) that carries a rocket launcher which constitute the subsequent stages. The control objective is to stabilize the aircraft in the launch phase. It is developed separately for two nonlinear MIMO motion modes of the model, the longitudinal mode and lateral mode, and is then applied to the full model of the aircraft. The considered model is highly nonlinear, mostly as a consequence of possible large angle of attack, sideslip and roll angle. Finally, the present work illustrates through simulations the good performance of the proposed control algorithm.

Introduction

Satellites launching is a strategical activity today. Launchers are able to carry from micro-satellites of some tens of kilograms up to 10 tons in the case of French Ariane 5 launcher. Recently, new applications have called upon very small satellites mostly used in groups (see [START_REF] Burlacu | An analysis of the nanosatellites launches between 2004 and[END_REF]). These small satellites need a new class of launchers since launching implies in many fixed costs that are independent of the size and weight of the launched device. For this reason, the ratio price-per-kilogram launched in space becomes too high. A quite logical solution in this case would be to pack many small devices to be launched together. Unfortunately this implies many additional risks in the split phase and is not envisaged. A more efficient solution in this case is to use the procedure of airlaunch (see [START_REF] Sarigul-Klijn | Selection of a carrier aircraft and a launch method for air launching space vehicles[END_REF], [START_REF] Gary | Quickreach responsive launch system[END_REF]). It consists of using a two stages launching system (see Fig. 1). The first stage is composed of an air vehicle (manned or unmanned) that carries (inside, beneath or above) a launcher which constitutes the second stage. There are many advantages in airlaunch, mainly because there is no need for specific large non populated launching areas. The aerial vehicle takes-off from a standard runway and fly to open ocean, avoiding populated areas or ship and airplane paths. For this reason there is also a minimization of weather constraints, since the vehicle can fly to open sky, and as consequence the launch delay can be significantly shortened. Similarly, instead of waiting for specific launch windows (to attain desired orbits), the vehicle may be flown to a better suited launch point, with a better alignment with the desired orbit. The fact that the first stage is a reusable aerial vehicle allows a much smaller launching delay. In the same way, launching reuse time may be very short (one or two days). These characteristics provide great flexibility, and allow to deploy small satellites designed for specific tasks of communication or data gathering in real time for urgent situations.

Airlaunch provides the advantages of two stage launchers. The second (dropped) stage may use specific nozzles and propellants for the low outside atmospheric pressure at altitude (20 000 meters or 60000 feet). This is obtained without the complex, expensive and relatively dangerous high pressure ground-launched first stage that is replaced by the aerial vehicle. Most current airlaunch projects use standard or lightly modified airplanes as first stage. For example, there has been tests using F15, C17, B52, L-1011 in Rascal, Quick-Reach, Proteus and Pegasus projects.

It is important to remark that airplanes use the wing's lift force to fly. For this reason, higher (low altitude) air density benefit the flight while the aircraft uses standard fuel to keep flying. A first stage rocket would use a much more complex, dangerous and expensive fuel while in this higher air density. From a certain altitude, air density is too low to be useful for an airplane, while not representing anymore a drawback for rockets.

Unlike the before-mentioned projects, other ones aim in developing an airlaunch system that uses an Unmanned Aerial Vehicle (UAV) instead of a standard aircraft with a human pilot inboard. The current paper addresses their launching phase, it intends to introduce modeling and a robust nonlinear controller for this delicate procedure. In fact, airlaunch may be very challenging since the rocket may be almost as heavy as the UAV. This means that the aircraft will instantaneously lose almost half of its mass. As consequence, the stabilization task is much more complex during and after the launching phase with a much more adverse mass ratio. To our best knowledge, it does not exist an equivalent research line, and then there is no results in the literature considering this problem.

The paper is organized as follows: in section 2, we describe the nonlinear mathematical system model. A conditional integrator control design and its application to the full system model is discussed in section 3. The paper is completed by computer simulations and conclusions. Modeling an air launch system and its separation phase is a difficult task which requires many data and informations about the real system. We have based our work on the model of an F-16 aircraft (that already has been tested for manned airlaunch). We have also used it to verify in simulations the reactions of the system when modeling the air launch phase.

Modeling

We suppose that the studied system is a set of a reusable launch vehicle and a down stage, whose mass is equivalent to the reusable launch vehicle's own mass. For the sake of simplicity, this set is considered as a complete aircraft before air launching. In this paper the airlaunch system is considered as a hybrid system composed by two continuous models that are switched in three phases, representing the system before, during and after the separation phase (see [START_REF] Van Cuong Nguyen | Modeling and lqr robust control of an unmanned aerial vehicle for airlaunch[END_REF]).

1. before the separation ⇒ the first aircraft model (representing the UAV and the rocket) is in a stable operating condition 2. during the separation ⇒ a second aircraft model representing only the UAV, starting on the previous operating condition, is disturbed by impulses on forces and moments representing a not perfect separation inside a time interval T int . Furthermore the initial conditions, inherited from the first phase, are not an equilibrium point for the second aircraft model.

3. after the separation ⇒ the disturbances stop (UAV and rocket are not in physical contact anymore).

It can be shown that the effect of launching the rocket from the UAV impacts most the lift force, and the roll and pitch moments. We suppose that these perturbing force and moments are constant during interval T int , and we represent then F zp , L p and M p for the perturbations on the lift force, on the roll moment and pitch moment respectively. In the present work we will study a worst case of disturbance. We consider that the separation phase is not simultaneous in all links that attach the rocket and the UAV. For this reason, the rocket remains attached to one end of the UAV during T int . We have then studied how long the disturbance could last and that the control algorithm could still stabilize the aircraft back.

We suppose that:

the perturbation on lift force during T int is equal to the air launch vehicle's mass, that means F wp = m r g cos θ 0 .

the perturbation on drag force is F up = -F wp sin θ 0 = -m r g sin θ 0 where θ 0 is the initial pitch angle of the first model at the launching phase and m r is the rocket's weight.

the perturbation on pitch moment during T int is an worst case that is represented by the rocket that remains attached to the aircraft by only one end during T int , applying a rotational movement to the aircraft, so a moment with value M p = m r gl r cos θ 0 /2 where l r is the rocket length.

the perturbation on roll moment during T int is small because of the rocket shape (long and thin).

the model following the launch phase is the F-16 model. Its initial condition is the state at an equilibrium point of the model previous the launch phase that is the F-16 model but with twice its standard mass. That means the rocket's mass is equal to the aircraft carrier. It is a hypothesis in the worst case because in practice the aircraft carrier's weight is greater than the rocket's weight.

The model of the dynamic airlaunch after the split phase is described by the Newton-Euler's law in the aerodynamic axes (OX W Y W Z W in Fig. 2) i.e. the reference frame attached to the airspeed vector V (see [START_REF] Lee | Nonlinear adaptive flight control using backstepping and neural networks controller[END_REF], [START_REF] Sonneveldt | Nonlinear F-16 Model Description[END_REF] and [START_REF] Peter | Modeling and Simulation of Aerospace Vehicle Dynamics[END_REF]), instead of using variables in the body fixed axes because of the measurability of these state variables:

                                           α = -cos α tan βp + q -sin α tan βr -sin α mV cos β (T + F x ) + cos α mV cos β F z + g V cos β [sin α cos θ + cos α cos φ cos θ] β = sin αp -cos αr -cos α sin β mV [T + F x ] + cos β mV F y -sin α sin β mV F z + g V [cos α sin β sin θ + cos β cos θ sin φ -sin α sin β cos φ cos θ] V = cos α cos β m [T + F x ] + sin β m F y + sin α cos β m F z + g[cos α cos β sin θ + sin β sin φ cos θ + sin α cos β cos φ cos θ] ṗ = 1 IxxIzz-I 2 xz [(I yy I zz -I 2 zz -I 2 xz )rq -I xz (I xx + I zz -I yy )pq + I zz L -I xz N ] q = 1 Iyy [(I zz -I xx )pr + I xz (p 2 -r 2 ) + M ] ṙ = 1 IxxIzz-I 2 xz [(-I xx I yy + I 2 zz + I 2 xz )pq + I xz (I xx + I zz -I yy )rq + I xx N -I xz L] φ = p + tan θ(q sin φ + r cos φ) θ = q cos φ -r sin φ ψ = q sin φ+r cos φ cos θ (1)
in which I xx , I yy , I zz , I xz are the moments of inertia, m is the mass of the system (kg) and g the gravity constant. α, β, V, p, q, r, φ, θ, ψ are the state variables of the airlaunch aircraft model, they are the angle of attack, sideslip, airspeed, roll rate, pitch rate, yaw rate, roll angle, pitch angle and yaw angle respectively. α, β, φ, θ, ψ are expressed in rad, p, q, r in rad/s and V in m/s. T is the thrust force, F x , F y , F z and L, M, N are aerodynamic forces and moments respectively. All forces and moments are expressed in N and Nm.

These aerodynamic forces and moments are function of all considered states. In this model, these aerodynamic forces and moments are under look-up table from wind tunnel data measurements as may be found in [START_REF] Nguyen | Simulator study of fighter airplane with relaxed longitudinal static stability[END_REF]. Finally, the control inputs are respectively the aileron (δ a ), rudder (δ r ) and elevator (δ e ) angles. all actuators are modeled as a first order model (τ = 0.0495s) with limit rates 60 • /s for aileron and elevator, and 120 • /s for rudder.

In particular, we use the low quality mode of the F-16 model, and the aerodynamic data is interpolated and extrapolated linearly in simulation from tables found in [START_REF] Nguyen | Simulator study of fighter airplane with relaxed longitudinal static stability[END_REF].

Control design 3.1 Conditional integrator control design

The MIMO conditional integrator (CI) controller design for the output regulation of a class of minimum-phase nonlinear systems in case of asymptotically constant references is studied in [START_REF] Mahmoud | Robust control for a nonlinear servomechanism problem[END_REF], [START_REF] Hassan | On the design of robust servomechanisms for minimum phase nonlinear systems[END_REF] and extended in [START_REF] Damm | Mimo conditional integrator control for a class of nonlinear systems[END_REF]. Our present work is dedicated to use a nonlinear extension of these results developed in [START_REF] Damm | Mimo conditional servocompensator control for a class of nonlinear systems[END_REF] and [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF], for stabilizing an unmanned aircraft during and immediately after the airlaunch phase.

This modified CI is composed of two main nonlinear terms, one saturated and one that may grow unbounded. The condition is responsible to choose where each term dominates the controller. In this way it behaves differently following some conditions (how big the errors are) that are tunned by the design parameters.

We remind the theory of conditional integrator control design for a nonlinear system.

Consider the system:

ė1 = e 2 ė2 = f (e 1 , e 2 ) + g(e 1 , e 2 )u (2)
where e 1 (t) ∈ R n is an output error vector, e 2 = ė1 , u ∈ R n control input and f (e 1 , e 2 ) ∈ R n , g(e 1 , e 2 ) ∈ R n×n are continuous functions.

Let us define an intermediate variable:

s = k 0 σ + K 1 e 1 + e 2 (3)
where σ ∈ R n is the output of the conditional integrator

σ = -k 0 σ + µsat(s/µ) (4)
in which µ is the boundary layer, k 0 is a positive parameter, K 1 ∈ R n×n is a positive definite matrix. The saturation function is determined as:

sat(s/µ) = s/ s if s ≥ µ s/µ if s < µ (5)
We denote O µ as the region neighborhood of (e 1 , e 2 ) = (0, 0) with a radius R µ for s < µ

O µ = {e = (e 1 , e 2 ) ∈ R n × R n | e ≤ R µ } (6)
We state the following assumptions on the forcing terms f (e 1 , e 2 ) and g(e 1 , e 2 ) to design the control algorithm.

∆ 0 = f (0, 0) f (e 1 , e 2 ) ≤ γ( e 1 + e 2 ) + ∆ 0 (7)
for (e 1 , e 2 ) ∈ R n × R n and while the sliding surface does not enter the boundary layer, i.e. s ≥ µ.

Inside the boundary layer, the function f (e 1 , e 2 ) is required to be Lipschitz for (e 1 , e 2 ) ∈ O µ , as a consequence

f (e 1 , e 2 ) -f (0, 0) ≤ L 1 K 1 e 1 + L 2 e 2 (8)
γ( e 1 + e 2 ) is also required to be Lipschitz for (e 1 , e 2 ) ∈ O µ :

γ( e 1 + e 2 ) ≤ γ 1 K 1 e 1 + γ 2 e 2 (9)
where L 1 and L 2 , γ 1 and

γ 2 ∈ R + . Assumption 3.2 g(e 1 , e 2 ) is invertible for (e 1 , e 2 ) ∈ R n × R n .
Following these assumptions, the controller u defined below in [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF] can be applied to (2) to stabilize the system. u = -π(e 1 , e 2 )sat(s/µ) [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF] in which,

π(•) = (π 0 + γ(•) + k 0 µ + ∆ 0 )g -1 (•) (11)
π 0 and k 0 are positive constants, µ is the boundary layer as defined above.

Remark 1 It is important to remark that the control can even grow unbounded since the term [START_REF] Van Cuong Nguyen | Modeling and lqr robust control of an unmanned aerial vehicle for airlaunch[END_REF] is not necessarily bounded. Functions γ(•) and g -1 (•) can grow continuously. For this reason we call this controller a modified Conditional Integrator, composed of two terms (see [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF][START_REF] Van Cuong Nguyen | Modeling and lqr robust control of an unmanned aerial vehicle for airlaunch[END_REF]) one saturated and one not. The later will dominate for small errors, and as a consequence the controller will behave as an integrator. In the case of large errors, it is the first that dominates, and the controller acts as a robust controller.

The stability in the general case of the control law [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF] for system (2) is demonstrated in [START_REF] Damm | Mimo conditional servocompensator control for a class of nonlinear systems[END_REF].

Lateral control design

In the lateral control design, we assume that all longitudinal state variables are null or constant, only lateral states are time varying. Moreover it is assumed that the airspeed's response is much slower than other states, and that the control surface deflection has no effects on the aerodynamic force components (lift and drag) but only on moments. Aerodynamic force F y and moments L, N are calculated by their aerodynamic coefficients (see more in [START_REF] Lee | Nonlinear adaptive flight control using backstepping and neural networks controller[END_REF]).

F y = (C y (β) + (C yp (α)p + +C yr (α)r) b/(2V ))qS L = (C l (β) + C lp (α, β)p b/(2V ) + C lr (α, β)r b/(2V ) + C l δa (α)δ a + C l δr (α)δ r )qS b N = (C n (β)+C np (α, β)p b/(2V )+C nr (α, β)r b/(2V )+C n δa (α)δ a +C n δr (α)δ r )qS b
By replacing F y , moments L, N and α = α 0 , θ = θ 0 in (1), the lateral nonlinear dynamic model used for the control design procedure is consequently reduced as: in which S is the wing area, q dynamic pressure, b is reference wing span,

                                 β= 1 mV (-cos (
I 3 = Izz (IxxIzz-I 2 xz ) , I 4 = Ixz (IxxIzz-I 2 xz ) , I 9 = Ixx (IxxIzz-I 2 xz ) . C y (α, δ e ), C yp (α 0 ), C yr (α 0 ), C l (α 0 , β), C n (α 0 , β), C lp (α 0 ), C np (α 0 ), C lr (α 0 ), C nr (α 0 ), C l δa (α 0 ), C n δa (α 0 ), C l δr (α 0 ), C n δr (α 0 )
are lateral aerodynamic coefficients taken from [START_REF] Morelli | Global nonlinear parametric modeling with application to f-16 aerodynamics[END_REF].

The previous equation can be rearranged as:

       β φ =f β 11 (β, φ) + f β 12 (β, φ) p r ṗ ṙ =f β 21 (β, φ) + f β 22 (β, φ) p r + g β 2 (β, φ) δ a δ r (13)
where

f β 11 (•), f β 12 (•), f β 13 (•), f β 21 (•), f β 22 (•)
, and g β 2 (•) represent the terms of ( 12) respectively (see Appendix A).

Let us define

x β 1 = [β, φ] T , x β 2 = ẋβ 1 = [ β, φ]
T and u β = (δ a , δ r ) T , that allow us to rewrite equation ( 13) to:

ẋβ 1 = x β 2 ẋβ 2 = F β (x β 1 , x β 2 ) + G β (x β 1 , x β 2 )u β (14)
where

             F β (•) = ( ∂f β 11 (•) ∂x β 1 + (f β (•) + f β 12 (•)f β 22 (•))(f β 12 (•)) -1 )x β 2 -(f β (•) + f β 12 (•)f β 22 (•))(f β 12 (•)) -1 f β 11 (•) + f β 12 (•)f β 21 (•) G β (•) = f β 12 (•)g β 2 (•) f β (•) [p, r] T = ∂(f β 12 (•)[p,r] T ) ∂x β 1 ( 15
)
We define an output error vector e β 1 = x β 1 -x β 1ref and e β 2 = ėβ 1 where x β 1ref = (β ref , φ ref ) T is the output reference considered as constant. Equation ( 15) can be transformed into (16) with two new state variables e β 1 and e β 2 .

ėβ 1 = e β 2 ėβ 2 = F β (e β 1 , e β 2 ) + G β (e β 1 , e β 2 )u β (16a) (16b)
G β (x β 1 , x β 2 ) is invertible in the considered domain of x β 1 = [β, φ] T , ẋβ 1 = x β 2 with β ∈ (-30 • , 30 • ) and φ ∈ (-180 • , 180 • ).
Application of control law [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF] for system (16) leads to the controller:

u β = -π β (e β 1 , e β 2 )sat(s β /µ β ) π β (•) = (π β 0 + γ β (•) + k β 0 µ β + ∆ β 0 )(G β (•)) -1 (17) with s β = k β 0 σ β + K β 1 e β 1 + e β 2 σβ = -k β 0 σ β + µ β sat(s β /µ β ) (18)
where π β 0 is a constant large enough, k β 0 is a positive parameter, µ β is the boundary layer and K β 1 is a positive definite matrix chosen such a way that K β 1 + sI 2 is Hurwitz. Theorem 3.1 System (16) with F β (•) which satisfies assumption 3.1, G β (•) satisfying assumption 3.2, and applying the control law (17-18), will globally reach an arbitrary error region in finite time, and there on will be exponentially stabilized towards its equilibrium point.

Proof: In order to demonstrate the exponential stability of designed controller ( 17) and ( 18) for the lateral mode in (16) which is a nonlinear MIMO system where the sideslip and roll angle are the outputs and aileron and rudder are the inputs, we will consider two regions: outside the boundary layer ( s β ≥ µ β ) and inside the boundary layer ( s β ≤ µ β ).

In the region s

β ≥ µ β , sat(s β /µ β ) = s β / s β .
The derivative of variable s can be expressed as:

ṡβ = k β 0 σβ + K β 1 ėβ 1 + ėβ 2 = -k β 0 s β + k β 0 µ β sat(s β /µ β ) + k β 0 (K β 1 e β 1 + e β 2 ) + K β 1 e β 2 + F β (•) + G β (•)u β
Now by letting

∆ β (•) = k β 0 (K β 1 e β 1 + e β 2 ) + K β 1 e β 2 + F β (•) (19)
Its derivative becomes:

ṡβ =-k β 0 s β + k β 0 µ β sat(s β /µ β ) + ∆ β (•) + G β (•)u β (20) Because of boundedness of F β (x β 1 , x β 2 , θ), ∆ β (•) is bounded by a function of γ β ( e β
1 + e β 2 ) (where γ β (•) is a class K function) and a positive constant ∆ β 0 (assumption 3.1):

∆ β (e β 1 , e β 2 ) ≤ γ β ( e β 1 + e β 2 ) + ∆ β 0 (21)
and as a consequence,

∆ β (e β 1 = 0, e β 2 = 0) = F β (0, 0) ≤ ∆ β 0 (22) for (e β 1 , e β 2 ) ∈ R n × R n . Let's consider the product (s β ) T ṡβ (s β ) T ṡβ = -(s β ) T k β 0 s β + k β 0 µ β (s β ) T sat(s β /µ β ) + (s β ) T ∆ β (e β 1 , e β 2 ) + (s β ) T G β (e β 1 , e β 2 )u β (23)
This product (s β ) T ṡβ can be developed with the definition of saturation function (5):

(s β ) T ṡβ =-(s β ) T k β 0 s β + µ β (s β ) T k β 0 s β / s β + (s β ) T ∆ β (•) -(s β ) T G β (•)π β (•)s β / s β ≤-(s β ) T k β 0 s β + µ β (s β ) T k β 0 s β / s β + ∆ β (•) s β -(s β ) T G β (•)π β (•)s β / s β ≤-(s β ) T k β 0 s β + µ β (s β ) T k β 0 s β / s β + (γ β (•) + ∆ β 0 ) s β -(s β ) T G β (•)π β (•)s β / s β ≤-(s β ) T k β 0 s β -(s β ) T (G β (•)π β (•) -(µ β k β 0 + γ β (•) + ∆ β 0 )In)s β / s β (24)
Replacing the control law in (17) and ( 18), the term (s β ) T ṡβ can be expressed as:

(s β ) T ṡβ ≤-(s β ) T k β 0 s β -(s β ) T (G β (•)π β (•) -(µ β k β 0 + γ β (•) + ∆ β 0 )In)s β / s β ≤-(s β ) T k β 0 s β -(s β ) T π β 0 s β / s β ≤-k β 0 s β 2 -π β 0 s β (25)
The product (s β ) T ṡβ is then not positive and we have also

d( s β 2 ) dt =2 s β d( s β ) dt = 2 (s β ) T ṡβ dt ≤ 2(-π β 0 s β -k β 0 s β 2 ) ∴ d( s β ) dt ≤-π β 0 -k β 0 s β ∴ s β (t) ≤ s β (0) -π β 0 t -s β (0) (e β ) -k β 0 t -1) (26)
Then s β (t) reaches the boundary layer µ β in finite time.

In the region

s ≤ µ β , sat(s/µ β ) = s/µ β .
Consider again ( 18) and (20), which inside the boundary layer may be rewritten as:

       σβ = -k β 0 σ β + s β ėβ 1 = -K β 1 e β 1 + s β -k β 0 σ β ṡβ = ∆ β (•) -G β (•)π β (•)s β /µ β (27a) (27b) (27c)
It can be shown that this system has an equilibrium point:

ēβ 1 = ēβ 2 = 0, s β = sβ , σ β = σβ with sβ = k β 0 σβ = µ β (π β (0, 0)) -1 (G β (0, 0)) -1 F β (0, 0) = µ β F β (0, 0)/(π β 0 + k β 0 µ β + ∆ β 0 )
. System (27) may be rewritten with respect to sβ and σβ :

       σβ = -k β 0 σβ + sβ ėβ 1 = -K β 1 e β 1 + sβ -k β 0 σβ ṡβ = ∆ β (•) -G β (•)π β (•)s β /µ β -G β (•)π β (•)s β /µ β (28a) (28b) (28c) where σβ = σ β -σβ , sβ = s β -sβ . F β (x β 1 , x β 2
) is a Lipschitz function inside the boundary region e.g. s β ≤ µ β , such that:

F β (e β 1 , e β 2 ) -F β (0, 0) ≤ l β 1 e β 1 + l β 2 e β 2 (29)
where l β 1 and l β 2 ∈ R + .

Following assumption 3.1, γ(•) is also a Lipschitz function:, such that:

γ β (•) F β (0, 0) ≤ (π β 0 + k β 0 µ β + ∆ β 0 )(γ β 1 e β 1 + γ β 2 e β 2 ) (30)
where γ β 1 and γ β 2 ∈ R + .

We would like to demonstrate that every trajectory of system (28) starting inside the boundary layer, will approach the equilibrium point as time tends to infinity when the control law (10) is applied. Toward that end, we take:

W β = λ β 1 2 (σ β ) T σβ + λ β 2 2 (e β 1 ) T e β 1 + (s β ) T sβ 2 (31)
as a Lyapunov candidate, where λ β 1 and λ β 2 are positive constants. Its derivative can be easily developed as:

Ẇ =λ 1 k β 0 (σ β ) T σβ + λ 2 (e β 1 ) T K β 1 ėβ 1 + (s β ) T ṡβ =λ 1 k β 0 (σ β ) T (-k β 0 σβ + sβ ) + λ 2 (e β 1 ) T K β 1 (-K β 1 e β 1 + sβ -k β 0 σβ ) +(s β ) T (∆ β (•) -G β (•)π β (•)s β /µ β -G β (•)π β (•)s/µ β ) (32) Since s β ≤ µ β , ∆ β (•) can be expressed as: ∆ β (•)=k β 0 sβ -(k β 0 ) 2 σβ -(K β 1 ) 2 e β 1 + K β 1 sβ -k β 0 K β 1 σβ + F β (•) (33)
Replacing system (28) and ∆ β (•) into the derivative of the Lyapunov function, we have then (reminding that sβ = µ β (π β (0, 0)

) -1 (G β (0, 0)) -1 F β (0, 0) and π β (0, 0) = (π β 0 + k β 0 µ β + ∆ β 0 )(G β (0, 0)) -1 ): Ẇ =λ 1 k β 0 (σ β ) T (-k β 0 σβ + sβ ) + λ 2 (e β 1 ) T K β 1 (-K β 1 e β 1 + sβ -k β 0 σβ ) +(s β ) T (k β 0 sβ -(k β 0 ) 2 σβ -(K β 1 ) 2 e β 1 + K β 1 sβ -k β 0 K β 1 σβ -G β (•)π β (•)s β /µ β ) +(s β ) T (F β (•) -G β (•)π β (•)s β /µ β ) =-λ 1 (k β 0 ) 2 (σ β ) T σβ + λ 1 k β 0 (σ β ) T sβ -λ 2 (e β 1 ) T (K β 1 ) 2 e β 1 + λ 2 (e β 1 ) T K β 1 (s β -k β 0 σβ ) +((s β ) T (k β 0 I n + K β 1 )s β + (s β ) T (k β 0 + K β 1 )k β 0 σβ -(s β ) T (K β 1 ) 2 e β 1 -(s β ) T G β (•)π β (•)s β /µ β ) +(s β ) T (F β (•) -F β (0, 0) - γ β (•) π β 0 +k β 0 µ β +∆ β 0 F β (0, 0))
Using equations (16a), (28b) and Assumption 29, we can have:

(s β ) T (F β (•) -F β (0, 0) - γ β (•) π β 0 +k β 0 µ β +∆ β 0 F β (0, 0)) ≤ s (l β 1 K β 1 e β 1 + l β 2 e 2 ) + γ β (•) π β 0 +k β 0 µ β +∆ β 0 s F β (0, 0) ≤ s (l β 1 K β 1 e β 1 + l β 2 e 2 ) + ∆ β 0 π β 0 +k β 0 µ β +∆ β 0 s (γ β 1 K β 1 e β 1 + γ β 2 e 2 ) ≤ s (l β 1 K β 1 e β 1 + l β 2 e 2 ) + s (γ β 1 K β 1 e β 1 + γ β 2 e 2 ) ≤ (l β 1 + γ β 1 ) s K β 1 e β 1 + (l β 2 + γ β 2 ) s e 2 ≤ (l β 1 +γ β 1 ) 2 ((s β ) T s + (e β 1 ) T (K β 1 ) 2 e β 1 ) + (l β 2 +γ β 2 ) 2 ((s β ) T s + e T 2 e 2 ) ≤ (l β 1 +γ β 1 ) 2 ((s β ) T s + (e β 1 ) T (K β 1 ) 2 e β 1 ) + (l β 2 +γ β 2 ) 2 ((s β ) T s +(s -k β 0 σβ -K β 1 e β 1 ) T (s -k β 0 σβ -K β 1 e β 1 )) ≤ (l β 1 +γ β 1 ) 2 ((s β ) T s + (e β 1 ) T (K β 1 ) 2 e β 1 ) + (l β 2 +γ β 2 ) 2 ((s β ) T s +3((s β ) T s + (k β 0 ) 2 (σ β ) T σβ + e β 1 (K β 1 ) 2 e β 1 )) ≤ 3(l β 2 +γ β 2 ) 2 (k β 0 ) 2 (σ β ) T σβ + (l β 1 +γ β 1 )+3(l β 2 +γ β 2 ) 2 ((e β 1 ) T (K β 1 ) 2 e β 1 ) + (l β 1 +γ β 1 )+4(l β 2 +γ β 2 ) 2 ((s β ) T s) ≤ c 1 (k β 0 ) 2 (σ β ) T σβ + c 2 (s β ) T s + c 3 (e β 1 ) T (K β 1 ) 2 e β 1 ( 34
)
where

c 1 = 3(l β 2 +γ β 2 ) 2 and c 2 = (l β 1 +γ β 1 )+3(l β 2 +γ β 2 ) 2 and c 3 = (l β 1 +γ β 1 )+4(l β 2 +γ β 2 ) 2
. Using (28) and (34), the derivative of the Lyapunov function is developed:

T = -k P (V -V ref ) -k I ( V -Vref )
where V ref is the airspeed reference, k P = 1242 and k I = 955.

Simulation Results

In section 3, the design methodology of the conditional integrator controller to stabilize the angle of attack, sideslip and roll angle is proposed when full knowledge of the aerodynamic characteristics is available. This section presents numerical simulation results for the controller to demonstrate the performance of the proposed conditional integrator control laws in the drop phase.

As mentioned in section 2, we have considered the launch phase as impulses on aerodynamic force and moments during a time interval T int , and that the model following the launch phase is that of an F-16. This model is used since it has already been applied for (manned) airlaunch, and because its nonlinear model, wind tunnel informations and data are widely known and used for control design. It is important to remark that the model used in the following simulations is even more complete than that used in the control design, for example it includes actuator dynamics and their limitations. As a consequence, simulations also illustrate some properties of robustness to unmodeled dynamics.

In the following simulations, we have simultaneously applied the SISO longitudinal controller for angle of attack, and the MIMO lateral one for the sideslip and roll angles in the full nonlinear F-16 aircraft model. We may note that the control inputs are limited by their physical bounds introduced in section 2.

The control law in [START_REF] Van Cuong Nguyen | Modeling and lqr robust control of an unmanned aerial vehicle for airlaunch[END_REF] whose π(•) can be written simpler as:

π(•) = (π 0 + γ(•))G(•) -1 (70) in which, γ(•) = γ 1 e 1 2 + γ 2 e 2 2
, γ 1 and γ 2 are positive constant. G(•) -1 is defined in (43) or in (16) depending on longitudinal mode or lateral mode.

Application of this control law to two motion modes presented in subsections (3.3) and (3.2) is done by determining the set of parameters π i 0 , γ i 1 , γ i 2 , µ i , K i 1 and k i 0 with i = α, β corresponding to longitudinal mode and lateral mode respectively. We use the design parameters in the Table 1 for the longitudinal controller, and the design parameters in Table 2 for the lateral controller. We will stabilize the second model following the launch phase to its equilibrium point (V, h) = (154m/s, 6500m) corresponding to angle of attack α r to 2.7 • , sideslip β r to 0 • , and roll angle φ r to 0 • . Its initial condition is the final state of the first model (α = 12.5 • , β = 0 • and φ = 0 • ) as in Section 2. Moreover, we add on its initial condition a small disturbance on system output. That means the initial condition of second model is (α 0 = 17.5 • , β 0 = 4 • and φ 0 = 10 • ) for all numerical simulations.

The second model is disturbed on its aerodynamic force and moments during an interval T int as in Section 2. 1. T int = 0.2s (corresponding to solid lines in Fig. 4 to Fig. 6) produces damped oscillations.

2. T int = 0.227s, the Conditional Integrator is able to stabilize the system as can be seen in the dashed lines in Fig. 4 to Fig. 6, but the bounds on actuator dynamics are reached.

3. T int = 0.3s (corresponding to dash dotted lines in Fig. 4 to Fig. 6) one can see that the controller is still able to stabilize the system.

Fig. 4 illustrates the convergence of the system output to the operating point of the aircraft at the end of 5s without static steady error for the three cases of T int = (0.2s, 0.227s, 0.3s). All system outputs are still under their physical limitation. Figs. 6 shows that angular rates converge to zero in all cases. In Fig. 5, it can be seen that the control variables are saturated by their physical limitations due to a high perturbation on aerodynamic forces and moments.

Finally we show in Fig. 7 and Fig. 8 that the system will be unstable for an interval T int = 0.39s.

Conclusion

We have developed a MIMO nonlinear Conditional Integrator Controller for the staging phase of an airlaunch system. Because our airlaunch system has a down stage mass close to the launch vehicle, the separation phase produces large changes in aerodynamic force and moment, as shown in section 2, which may turn the system unstable. These changes may, specially in the case of a non perfect separation phase, produce large impulses on forces and moments, making the system unstable. These impulses are considered to last a time interval, that is then evaluated in simulations. To stabilize the airlaunch system after this stage separation phase, a conditional integrator control is considered. This controller is designed using an F-16 MIMO nonlinear model representing the aircraft during and after dropping the second stage, and as a consequence disturbed by large impulses.

For a perturbation on aerodynamic forces and moments during an interval T int , the stability of the system after the drop stage may be assured for small time intervals. When T int becomes large, the system becomes unstable even with the proposed controller, that is bounded by actuator's limitations.

Figure 1 :

 1 Figure 1: Stratolaunch system credited © Stratolaunch Systems, Inc.

Figure 2 :

 2 Figure 2: Frames: Body fixed axes OX B Y B Z B , Stability axes OX S Y S Z S , Aerodynamic axes OX W Y W Z W

Figure 3 :

 3 Figure 3: The launcher attached to the aircraft carrier in the worst case

  This model is based on wind tunnel data from NASA, considering the following conditions: angle of attack is in the range of [-10 • , 45 • ] and sideslip of [-30 • , 30 • ] flap deflection is ignored physical constraints for aileron (|δ a | ≤ 21.5 • ), rudders (|δ e | ≤ 25 • ) and elevator (|δ r | ≤ 30 • )

  Assumption 3.1 f (e 1 , e 2 ) is bounded by a function γ( e 1 + e 2 ) of class K and a positive constant

  α0) sin(β)(T + Cx(α0)qS) + cos(β)Cy(β)qS -sin(α0) sin(β)Cz(α0, β)qS) + sin(α0)p -cos(α0)r + ρS 4m (cos(β)Cy p (α0) bp + cos(β)Cy r (α0) br) + g V (cos(α0) sin(β) sin(θ0) + cos(β) cos(θ0) sin(φ) -sin(α0) sin(β) cos(φ)) φ=p + cos(φ) tan(θ0)r ṗ =I3C l (α0, β)qS b + I4Cn(α0, β)qS b + ρV S b 4 [(I3C lp (α0) +I4Cn p (α0))p + (I3C lr (α0) + I4Cn r (α0))r] +qS[(I3C l δa (α0) + I4Cn δa (α0))δa + (I3C l δr (α0) + I4Cn δr (α0))δr] ṙ =I4C l (α0, β)qS b + I9Cn(α0, β)qS b + ρV S b 4 [(I4C lp (α0) +I9Cn p (α0))p + (I4C lr (α0) + I9Cn r (α0))r] +qS[(I4C l δa (α0) + I9Cn δa (α0))δa + (I4C l δr (α0) + I9Cn δr (α0))δr] (12)

Figure 4 :

 4 Figure 4: System Output: Angle of attack, Sideslip Angle and Roll angle stabilized

Figure 5 :

 5 Figure 5: Control surface: Aileron, Elevator and Rudder

Figure 6 :

 6 Figure 6: State variables: Angular rates of system

Figure 7 :

 7 Figure 7: Angle of attack, Sideslip Angle and Roll angle unstable

Figure 8 :

 8 Figure 8: Saturation of Aileron, Elevator and Rudder

Table 2 :

 2 Parameters for the lateral mode

It can be verified that by taking λ β 1 , λ β 2 and π β (•) large enough and µ β small enough, the derivative of Lyapunov function is negative. We establish the additional design parameter's condition:

Inequality (36) implies that the design condition of parameter k β 0 and matrix K β 1 must satisfy:

In this way, W β (t) satisfies W β (t) > 0 and Ẇ β < -w 0 W β (where w 0 is a positive constant) for all σ β = σβ , e β 1 = 0 and s β = sβ . Then W β (t) reaches exponentially zero when time tends to infinite. As consequence, the output error e β 1 (t) tends to zero, σ β and s β tend to their equilibrium values as time tends to infinite. We may assure the exponential stability of the system in the region of s β ≤ µ β .

Longitudinal control design

As in the case of lateral control design, in the longitudinal case it is considered that only longitudinal state variables are time varying. It is a single input single output system where angle of attack is the output and elevator is the input. Aerodynamic forces F x , F z and moment M can be calculated by its aerodynamic coefficients (see more in [START_REF] Lee | Nonlinear adaptive flight control using backstepping and neural networks controller[END_REF]).

By replacing F x , F z , moment M and β = 0, φ = 0, p = 0, r = 0 in (1), the model for longitudinal dynamic can be written as:

are aerodynamic coefficients taken from [START_REF] Morelli | Global nonlinear parametric modeling with application to f-16 aerodynamics[END_REF].

Equation (38) can be rearranged as:

where f α 11 (α), f α 12 (α), f α 13 (α, θ), f α 21 (α), f α 22 (α) and g α 2 (α) represent the terms of (38) respectively (see Appendix B).

Let us define x α 1 = α, x α 2 = ẋα 1 = α and u α = δ e , which allow us to rewrite (39) into:

where

2 and θ, it can be expressed as:

Remark 2 We have some remarks:

22 and g α 2 are function of aerodynamic coefficients under analytical forms by interpolation from wind tunnel test data. F α and h , formed from these functions, can be then bounded by a class K function and be a Lipschitz function in a flighting envelop.

As a consequence, they fulfill Assumptions 3.1 and 3.2.

We define now the reference for the angle of attack α ref considered as constant in this study, and the error vector of angle of attack e α 1 = x α 1 -x α 1ref = α -α ref and the variable e α 2 = ėα 1 . Equation (40b) can be transformed into:

Since function h(•) depends on θ under cosines and sinus functions, it is easy to show that h(•) is bounded by a function of class K γ α 1 (•) and a constant

Here we remark that G α (x α 1 , x α 2 ) is invertible for the entire domain of attraction of the system. It is important to remark that pitch angle θ is not included in the considered system, and is left free in the analysis that follows. The reason for this is to acknowledge the possibility of airlaunching under different possible angles θ, and even under a looping like trajectory that naturally brings the aircraft away of the rocket.

Applying the control algorithm presented in [START_REF] Van | Modeling and conditional integrator control of an unmanned aerial vehicle for airlaunch[END_REF] for system (43) which in this case is a nonlinear single input single output system, gives the controller:

where π α 0 , µ α , K α 1 and k α 0 are positive constants.

Theorem 3.2 System (43) with F α (•) satisfying assumption 3.1, G α (•) which satisfies assumption 3.2, and the control law (45-46), will globally reach an arbitrary error region in finite time, and there on will be exponentially stabilized towards its equilibrium point.

Proof: In order to demonstrate the exponential stability of designed controller (45) and ( 46) for the longitudinal mode in (43), we will consider in two regions outside the boundary layer (|s α | ≥ µ α ) and inside the boundary layer (|s α | ≤ µ α ). The longitudinal mode in this study is a single input single output system.

In the region |s

The derivative of the integral error measurement surface can be then expressed as:

From ( 43) and ( 46), the previous equation may be written again :

Now by letting

The derivative of the integral error measurement surface becomes:

The term

) outside the boundary region e.g. |s α | ≥ µ α and a positive constant f α 0 , where

Function h(•) is also bounded in the same way (see (44)). Then δ α (•) is bounded by a γ α (•) class K function and a positive constant δ α 0 :

where

and as a consequence,

Let's consider the product s α ṡα (since s is a scalar)

1 , e α 2 ) + s α g α (e α 1 , e α 2 )u α (53) This product (s α ) T ṡα can be developed with the definition of saturation function (5):

Replacing the control law in (45) and ( 46), the term s α ṡα can be expressed as:

The product s α ṡα is then not positive and we have also

Then the integral error measurement surface s α (t) reaches the boundary layer µ α in finite time.

In the region |s|

Consider again (46) and (49), which inside the boundary layer may be rewritten as:

It can be shown that this system has an equilibrium point: ēα 1 = ēα 2 = 0, s α = sα , σ α = σα with s = k α 0 σα = µ α (π α (0, 0)) -1 (g α (0, 0)) -1 (f α (0, 0) + h(0, 0, θ)) = (f α (0, 0) + h(0, 0, θ))/(π α 0 + k α 0 µ α + δ α 0 ). System (57) may be rewritten with respect to sα and σα :

where σα = σ α -σα , sα = s α -sα .

Function

is Lipschitz inside the boundary region e.g. |s α | ≤ µ α (see remark 2), such that:

where l α 1 and l α 2 ∈ R + .

γ α 2 (•) and γ α 1 (•) is are chosen so that are differentiable, then Lipschitz. As a consequence, function γ α (•) is then Lipschitz.

where

We would like to demonstrate that every trajectory of system (58) starting inside the boundary layer, will approach the equilibrium point as time tends to infinity when the control law (45) is applied. Toward that end, we take:

as a Lyapunov candidate, where λ α 1 and λ α 2 are positive constants. Since |s α | ≤ µ α , δ α (•) can be expressed as:

From ( 58) and reminding that s = (f α (0, 0) + h(0, 0, θ))/(π

, the derivative of the Lyapunov function is developed:

) Using equations (43a) and (58b), the relation in (59) and in (60) can be expressed as:

where

For any value of θ, it is very important to remark that this variable was not included in this analysis, even if its derivative is. We will discuss this point again later, but the main reason is that the airlaunch may be performed under any pitch angle θ, and even under a looping-like trajectory. From (65), the derivative of the Lyapunov function can be further developed:

The derivative of Lyapunov function negative implies that λ α 1 , λ α 2 and π α 0 large enough and µ α small enough, to satisfy:

In this way, W α (t) satisfies W α (t) > 0 and Ẇ α < -w 0 W α + Cµ (where w 0 is a positive constant) for all σ α = σα , e α 1 = 0 and s α = sα . Then W α (t) is ultimately bounded towards a neighborhood of zero when time tends to infinite, independent of the pitch angle θ. As consequence, e α 1 (t), σ α and s α tend to a region around their equilibrium point as time tends to infinite and for every pitch angle θ. The system is then said to be ultimately bounded to a small region which is function of pitch angle θ. Finally, it is interesting to remark that, as said before, variable θ is left free in order to allow situations as a looping, where θ is continuously varying. Its derivative on the other hand is bounded, and also goes to a residual set. In fact, the best trajectory for airlaunch is still an open problem. There are proposals of launching in horizontal, constant angle (climbing), concave or convex (zero gravity instant) trajectories.

Since the airspeed control is only a secondary objective, we design a simple PI controller for the thrust to regulate airspeed. Its form is:

A Functions in the lateral mode


Khadidja Benzemrane

Gilney Damm

Giovanni L Santosuosso

Nonlinear Observability of Aerial Vehicles without GPS Measurements

The present work addresses the problem of estimating the translation speed vector of an unmanned aerial vehicle only based on measurements of angles, angular speeds and translation accelerations. This is a delicate problem since this kind of system is commonly known as unobservable. In order to build this result we have first demonstrated that this problem, even if the linearized model is unobservable, the considered nonlinear system is in fact unobservable only inside a domain, but observable in another state space region. This result has been established by the study of the nonlinear observability of the considered system. In addition we briefly introduce a nonlinear adaptive observer that allow the estimation of translation speeds even in the case of outputs disturbed by measurement noise, wind or model errors, with rigorous theoretical proofs of the proposed results. The designed observers, based on adaptive control techniques, need the well known condition of persistency of excitation. Finally, we have shown that this condition in fact brings the system to an observable state space manifold. The considered manifold is easily reachable for a wide class of unmanned aerial vehicles including helicopters and quadrotors since these manifolds depend on the measurable and fully controllable states of this class of system, and as consequence, independent of the observer's results.

I. INTRODUCTION

The objective of this work is to study the possibility of estimating the translation speed vector of an unmanned aerial vehicle only based on measurements of angles, angular speeds and translation accelerations. This is a well known problem in robotics and mechanics in general, and was classically known as non-observable.

This estimation may be constructed based on other methods using different sensors like cameras and lasers, relying on probabilistic and statistic techniques like sensor fusion. These techniques on the other hand can not be applied to small drones where size, payload, energy consumption and computer resources are extremely limited. Furthermore, these techniques do not provide rigorous convergence proofs and leave unanswered questions about the fundamental concepts of system's observability. For all these reasons, the proposed work does not consider these techniques. In the same way, in many situations (indoors, urban or close to large structures for example) it is not possible to use GPS related devices, and those are also not considered in the present work.

Many recent articles were dedicated to the (nonlinear) estimation of state variables and external or internal parameters Khadidja Benzemrane and Gilney Damm are with IBISC Laboratory -Université d'Evry Val d'Essonne, Evry, France [khadidja.benzemrane,gilney.damm]@ibisc.fr Giovanni L. Santosuosso is with the Dipartimento di Ingegneria Elettronica, Universitá di Roma Tor Vergata, Rome, Italy santosuosso@ing.uniroma2.it of unmanned aerial vehicles (see [START_REF] Bonnabel | A non-linear symmetrypreserving observer for velocity-aided inertial navigation[END_REF], [START_REF] Cheviron | Robust nonlinear fusion of inertial and visual data for position, velocity and attitude estimation of UAV"journal=In practice[END_REF], [START_REF] Pflimlin | Nonlinear attitude and gyroscope's bias estimation for a VTOL UAV[END_REF], [START_REF] Vissiere | Using distributed magnetometers to increase IMU-based velocity estimation into perturbed area[END_REF], [START_REF] Zhao | Discrete nonlinear observers for inertial navigation[END_REF]). But in the best knowledge of authors, none has addressed the question of system observability, and in particular the translational speed. This part of the state vector was in general considered as non-observable, and most examples rely on the direct integration of the acceleration coupled to the measurements of other sensors by a sensor fusion algorithm. This scheme is very sensitive to measurement noises and uncertainties and in practice presents large errors.

In previous works, ( [START_REF] Benzemrane | Développement d'observateurs et application au drone XSF[END_REF], [START_REF] Benzemrane | Adaptive Observer and Kalman Filtering[END_REF], [START_REF] Benzemrane | Adaptive Observer and Kalman Filtering[END_REF], [START_REF] Benzemrane | Unmanned aerial vehicle speed estimation via nonlinear adaptive observers[END_REF]) we had proposed an estimation scheme based on nonlinear adaptive techniques. Under some well known, and easily fulfilled, constraints of persistency of excitation this estimator may in fact obtain exponentially stable estimations of the state space for the undisturbed system. This is a quite unexpected result for a classically known unobservable system. This lead us to a deeper study of the system's observability. In fact, the tangent linearization of the considered system is non-observable, but a study of observability for the nonlinear system shows that this is only valid in part of the state space. This means that there exist observable and unobservable manifolds. Furthermore, the persistency of excitation condition brings the system to an observable manifold, and this explains why it is possible to design an observer to the considered system. This analysis was only made using the definition of local observability, what implies in analysing the observability close to some trajectories of states. In those results it was shown that a sinusoidal trajectory was inside the observability space of a flying object.

In the present work, we will extend this result using the nonlinear observability definition, but will extend our study on a family of manifolds created by sinusoidal inputs of different frequencies and amplitudes.

II. OBSERVABILITY

The concept of observability of nonlinear systems is a quite delicate matter. In most cases the standard procedure is to linearize the complete system around an operating point, and then analyse its observability by the standard procedure for linear systems.

Considering a generic unmanned aerial vehicle (UAV), its dynamic model may be expressed as follows, where η 1 [x, y, z] T be the UAV position vector represented in the global reference frame, let η 2 [φ, θ, ψ] T be the Euler angles vector represented in the global reference frame (roll pitch and yaw respectively) ; we also assume that ν 1 [u, v, w] T is the translational speed vector represented in the local reference frame (surge, sway and heave respectively) and ν 2

[p, q, r] T is the angular speed vector represented in the local reference frame. The UAV model in state space form collects the set of the first order differential equations expressing η1 , η2 , ν1 , ν2 , as a function of η 1 , η 2 , ν 1 , ν 2 . To be more specific, the linear velocity in the fixed reference frame η1 = ( ẋ, ẏ, ż) T can be expressed as a function of η 2 , ν 1 , via the set of differential equalities

ẋ = cosθ cosψ u + (sinφ sinθ • cosψ -cosφ sinψ) v + (cosφ sinθ cosψ + sinφ sinψ) w ẏ = cosθ sinψ u + (sinφ sinθ sinψ + cosφ cosψ) v + (cosφ sinθ sinψ -sinφ cosψ) w ż = -sinθ u + sinφ cosθ v + cosφ cosθ w (1)
The time derivatives of the Euler angles η2 =( φ, θ, ψ) satisfy the following equations :

φ = p + (sinφ q + cosφ r) tanθ θ = cosφ q -sinφ r ψ = (sinφ q + cosφ r) cosθ -1 (2)
We assume that X, Y and Z are the aerodynamic forces applied by the motors to the UAV such that :

m u = m(-qw + rv) + X -m g sinθ m v = m(-ru + pw) + Y + m g sinφ cosθ m ẇ = m(-pv + qu) + Z + m g cosφ cosθ (3)
where g is the gravity acceleration constant, m is the drone mass. In the angular setting we obtain

I xx ṗ = L -(I zz -I yy ) q r I yy q = M -(I xx -I zz ) r p I zz ṙ = N -(I yy -I xx ) p q (4)
L, M and N being the aerodynamic torques provided by the motors. This model consider, without less of generality, a simplification of the inertia matrix in respect to symmetry plans.

The present problem considers as available for measurements the orientation (Euler Angles), angular velocities and acceleration of the vehicle. In this context 1 , defining the system state vector

   S = [φ, θ, ψ, p, q, r, u, v, w] T y = [φ, θ, ψ, p, q, r, u, v, ẇ] T Ω = [X, Y, Z, L, M, N ] T (5)
equations ( 2)-( 4) can be written in a compact notation as

Ṡ = f (S, Ω), S ∈ R 9 y m = h(S, Ω), y ∈ R 6 (6)
where f : R 9 × R 6 → R 9 is a smooth function depending on the state and the control input Ω, and h(S, Ω) is the second 1. We will not consider position estimation. For this reason, the state vector does not include it. equality in [START_REF] Benzemrane | Unmanned aerial vehicle speed estimation via nonlinear adaptive observers[END_REF] where u, v, ẇ are expressed as functions of S and Ω via (3).

The paper objective may then be stated as : given system (6) study if it is possible to design observation strategies yielding exponentially converging estimates of the state S from the measurable output y, to be eventually included in any control strategy of system [START_REF] Bonnabel | A non-linear symmetrypreserving observer for velocity-aided inertial navigation[END_REF].

In practical terms, from [START_REF] Benzemrane | Unmanned aerial vehicle speed estimation via nonlinear adaptive observers[END_REF] we deduce that the observation algorithm focuses on the estimation of the unmeasurable state variables ν 1 , based on the measurable variables η 2 , ν1 and ν 2 given by the standard sensors embedded in the drone.

Observability of linearized system

The standard procedure to study the observability of a nonlinear system is to linearize the model around an equilibrium point and to use the linear techniques to study this linearized model.

In the case of the estimation of the translational speed of an aerial vehicle, if one does the linearisation of the system around an equilibrium point S =   η2 ν1 ν2   given by :

η 2 ⇒ η2 =   φ θ ψ   ν 2 ⇒ ν2 =   0 0 0   ν 1 ⇒ ν1 =   ū v w  
will obtain an standard observability matrix given by :

O =      C C A . . . C A n-1     
that in the present case has a rank equal to 6. As a consequence the system is considered as not observable.

Nonlinear Observability

The design of state observers for linear systems has been extensively addressed in the literature. In that case, it is easy to check if a system is observable, and this property is valid for all state space. For a nonlinear system, the fact that the system is observable in a region does not means that this property is kept for the whole space. This means that two distict states x 1 and x 2 produce the same output y(x 1 ) = y(x 2 ) = ȳ for a given input Ω. This follows from the standard definition of distinguishability (Isidori and other references).

In fact, the observability of a nonlinear system is also dependent of the system's inputs as well as the states. One simplified procedure to verify this observability is to check the local observability of the system for constant inputs. To this purpose, we recall the following definition (see [?]).

Definition 1: Let G denote the set of all finite linear combinations of the Lie derivatives of h 1 , . . . , h p with respect to f for various values of Ω = constant. Let dG denote the set of all their gradients. If we can find n linearly independent vectors within dG, then the system is locally observable.

The observability matrix O = dG is given by :

O =   L 0 f h ... L p-1 f h  
The system is locally observable if O has full rank.

This method allows to find observable and unobservable points and trajectories, but is very limited. We have applied this method in ( [START_REF] Benzemrane | Adaptive Observer and Kalman Filtering[END_REF]) to establish the existance of observable trajectories for the nonlinear in the case of pure sinusoidal inputs in the three aerodynamic torques provided by the motors. Unfortunatelly, nothing could be said in respect of the effects of sinusoidal amplitudes and frequencies. In the present work we extend this result in order to better characterise the effect of these parameters.

The first step was to augment the system states by an exciting exosytem. We were inspired by recent results from A. Isidori like [START_REF] Byrnes | Nonlinear internal models for output regulation[END_REF] and [START_REF] Byrnes | Further results on output regulation by pure error feedback[END_REF] for example. We then define the exosystem whose states are ζ L , ζ M and ζ N , with dynamics described by :

ζL = 0 f 1 -f 1 0 ζ L = A L ζ L ζM = 0 f 2 -f 2 0 ζ M = A M ζ M ζN = 0 f 3 -f 3 0 ζ N = A N ζ N (7)
Variables ζ L , ζ M and ζ N are considered available for measurement, f 1 , f 2 and f 3 stand for the sinusoid's frequencies. Amplitudes and phases on the other hand are a function of initial conditions L (0), ζ M (0) and ζ N (0).

These oscillators are then used to provide torques L, M and N applied to the aerial vehicle. These torques may be defined as :

L * = 1 0 ζ L 1 0 M * = 1 0 ζ M 1 0 N * = 1 0 ζ N 1 0 (8)
In this way, the extended system may be expressed as :

                    φ θ ψ ṗ q ṙ u v ẇ ζL ζM ζN                     =                   (sin(φ) q + cos(φ) r) tan(θ) cos(φ) q -sin(φ) r (sin(φ) q + cos(φ) r) cos(θ) -1 -(Izz -Iyy)rq -(Ixx -Izz)pr -(Iyy -Ixx)pq (-qw + rv -g sinθ) (-ru + pw + g sinφ cosθ) (-pv + qu + g cosφ cosθ) 0 0 0                   + B                   φ θ ψ p q r u v w ζ L ζ M ζ N                   +                   0 0 0 L * M * N * X Y Z 0 0 0                   y = C                   φ θ ψ p q r u v w ζ L ζ M ζ N                   +                   0 0 0 0 0 0 (-qw + rv -g sinθ) + X (-ru + pw + g sinφ cosθ) + Y (-pv + qu + g cosφ cosθ) + Z 0 0 0                   B =                  
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A L 0 0 0 0 0 0 0 0 0 0 0 0 A M 0 0 0 0 0 0 0 0 0 0 0 0

A N                   C =                  
1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1

                 
In this new extended system, the inputs where replaced by states, and the observability study, even if restricted to pure sinusoidal inputs, is much easier. This study has shown that the system is observable if vectors ζ L , ζ M and ζ N are not identically equal to zero. In the beginning, we have studied the case where the three sinusoids were of diferent frequencies and amplitudes. In that case, system's observability was assured, confirming our previous works. A more interesting, and in some sense unexpected, result is that even in the case where frequencies and/or amplitudes are not different, the observability is still assured. Since the observers developed in ( [START_REF] Benzemrane | Adaptive Observer and Kalman Filtering[END_REF], [START_REF] Benzemrane | Adaptive Observer and Kalman Filtering[END_REF] and [START_REF] Benzemrane | Unmanned aerial vehicle speed estimation via nonlinear adaptive observers[END_REF]) were based on adaptive techniques, we expected to find the standard "rule of the tumb" of distinct sinusoids in order to attain persistance of excitation. But the present work shows that even only one frequency brings the system to an observable manifold. We can then conclude that each sinusoidal input makes system's dynamics evolve on an observable manifold, and that there exists a family of manifolds parametrized on input's amplitudes and frequencies.

III. OBSERVER DESIGN

Notice that the measured output can be expressed as

y = [y 1 , y 2 ] T = [φ, θ, ψ, p, q, r y1 , u, v, ẇ y2 ] T
where the vector y 2 = ν1 = [ u, v, ẇ] T is available for measurement, so that equation ( 3) is re-written as

ν1 (t) = Φ(t)ν 1 (t) + Λ(t) y 2 (t) = ν1 (t) (9)
with the matrix Φ as :

Φ(t) =   0 r -q -r 0 p q -p 0   ( 10
)
and the vector Λ defined as

Λ(t) =   -g sinθ + X g sinφcosθ g cosφ cosθ + Z   (11)
In the following we design a reduced-order observer which estimates the translational velocity of the UAV recalling that, in order to guarantee that p(t), q(t) and r(t) are different from zero, the control inputs have to be chosen such that the system oscillates around the equilibrium point. We introduce an estimate of the linear velocity :

ν1 =   û v ŵ  
The observer described by the following equations :

ν1 = ν1 + KΦ T ζ ζ = -cζ + Φν 1 (12)
whith a constant c > 0, a diagonal constant matrix K ∈ R 3 × R 3 , K > 0, a vector ζ that is an internal variable, for all initial conditions ν ∈ R 3 and ζ ∈ R 3 garantees, under some conditions on matrix Φ that will be introduced during the proof, the convergence to zero of estimation error ν1 = ν 1 -ν1 that satisfies the following equation :

ν1 = -KΦ T ζ (13)
Proof:

The stability proof is standard and follows Lyapunov techniques. Consider the Lyapunov function :

V = 1 2 ζ T ζ + 1 2 νT 1 K -1 ν1 (14)
By computing the time derivative of V , we have :

V = -c ζ 2 + ζ T Φν 1 + νT 1 K -1 (-KΦ T ζ) = -c ζ 2 + νT 1 Φ T ζ + νT 1 Φ T ζ V = -c ζ 2 (15)
Then, V ≤ 0. The previous equations imply that ζ and ν1 are bounded and that 0 ≤ V (t) ≤ V (0). Here we introduce the Hypothesis (1) that Φ(t) is bounded, and then the second equation of system [START_REF] Nijmeijer | Nonlinear dynamical control systems[END_REF] assures that ζ is also bounded, and then

V = -cζ ζ is bounded.
As a consequence V is uniformly continuous, and according to Barbalat's Lemma, since V ≥ 0, V tends to zero. This in turns implies, from [START_REF] Vissiere | Using distributed magnetometers to increase IMU-based velocity estimation into perturbed area[END_REF], that ζ also converges to zero. Then, according to [START_REF] Pflimlin | Nonlinear attitude and gyroscope's bias estimation for a VTOL UAV[END_REF] and the convergence of ζ to zero, results that ν1 also converges to zero. Now, stating the Hypothesis (2) that Φ is also bounded, it is shown using the equation :

ν1 = Λ(Φ T ζ + ΦT ζ)
that ν1 is also bounded and then ν1 is uniformly continuous. As a consequence it can now be stated that ν1 converges to a constant, since ν1 converges to zero.

From the second equation of system [START_REF] Nijmeijer | Nonlinear dynamical control systems[END_REF], since ζ converges to zero, and its forcing term Φν 1 is uniformly continuous :

Φν 1 = 0 (16)
Multiplying equation ( 16) by Φ(t), and integrating it by a time period T , one arrives to :

t+τ t Φ(s) Φ T (s)ds ν1 = 0 (17)
It was shown that ν1 converges to a constant, and as a consequence this equation has a unique zero solution if and only if the matrix Φ(t) satisfies the persistency of excitation condition, which means that there exist strictly positive constants α and T such that for any t > t 0 :

t+T t Φ(τ )Φ T (τ )dτ ≥ kI (18)
where k > 0 (Hypothesis (3)).

Then the observer [START_REF] Nijmeijer | Nonlinear dynamical control systems[END_REF] assures global assymptotical estimation of the unmeasurable states representing the linear velocities of the considered vehicles.

Remark 1: Hypothesis (1-3) are easily satisfied for a large class of aerial vehicles like helicopters, quadrotors ans PVTOLs, since motors assure the boundness of the angular speeds [p, q, r] T with bounded control magnitude.

Remark 2: Hypothesis (3) reflects the observability study presented in section II. The persistance of excitation condition comes from the adaptive control field, and is not necessarily related with the concept of observability. It is interesting to remark that in the present case, both concepts (observability and persistance of excitation) coincide.

IV. SIMULATION RESULTS

In this section we illustrate the observer designed above to estimate the linear velocity of an quadrotor UAV (see [START_REF] Azouz | Modeling and development of a 4rotors helicopter UAV[END_REF]) based on measurements of angles, angular velocities and linear accelerations.

In all simulations the observer estimations have initial conditions set to zero, while the desired states are time varying (and different from zero at t = 0). In all simulations, we have used the following parameters values : The orientations and angular speeds of the UAV are plotted on Fig. 1 and2. Small periodic orbits were considered in order to satisfy the persistency of excitation condition (equation ( 18)). Some constraints must also be satisfied to garantee the stability of the physical system : the angles have to be less than 0,1 rad (about 6 degrees) and the frequency of the needed oscilations must be less than 1Hz.

We have considered the case of additive measurement noise on the accelerations. The measured acceleration is presented in Fig. 3 where the noise µ is 10 percent of the measured acceleration.

The time history of the three speeds to be observed is shown in Fig. 4(a), and observer estimates in Fig. 4(b). Finally, Fig. 4c illustrates the observation error going to a residual set given by the noise amplitude. Pursuiting previous results, this work has shown that in the case of simple sinusoidal inputs in torques will produce a family of observable manifolds, parametrised in frequency and amplitude of the considered sinusoids. Furthermore, in oposition as what would be expected, the observability property does not need different frequencies and amplitudes in the inputs. Finally, we have shown that the condition of persistance of excitation needed for a proposed observer indeed brings the system to an observable region.

The characterization of system's observability allows the developpement of other observer algorithms since the system is brougth to an observable region. This is easily achievable for a large class of unmmaned aerial vehicles like helicopters, quadrotors and PVTOLs. It is important to notice that in these vehicles, angular dynamics are fully measurable and controlable. As a consequence, it is possible to bring these systems to an observable region independetly of observer convergence.

m = 2

 2 .500Kg l b = 23cm I R = 100.10 -7 Kgm 2 I xx = 224931.10 -7 Kgm 2 I yy = 222611.10 -7 Kgm 2 I zz = 325130.10 -7 Kgm 2 k T = 10 -5 N.s 2 k M = 9.10 -5 N s 2 m I R is the inertia moment of the rotating elements in the rotors, while I xx , I yy and I zz are the diagonal terms of the body's inertia tensor. Finally k T represents a coefficient related to the lift torque (F k T ω 2 ) and k M is the coefficient related to the drag torque (M D k M ω 2 ).

  Fig. 1. Orientation of the UAV

Fig. 4 .

 4 Fig. 4. Linear speed (a), Speed estimation (b) and Estimation error (c)

TABLE I SIMULATION PARAMETERS

 IPARAMETERS