Daniele Pucci
email: daniele.pucci@iit.it

Tarek Hamel
email: thamel@i3s.unice.fr

Pascal Morin
email: morin@isir.upmc.fr

Claude Samson
email: csamson@i3s.unice.fr

Nonlinear feedback control of axisymmetric aerial vehicles

Keywords: aerodynamic forces, aircraft control, axisymmetric body, flight control, stabilizing feedback

come

Introduction

Flight control techniques have been developed primarily for full-size commercial airplanes, designed to fly along very specific trajectories (trim trajectories with a very narrow range of angles of attack). Control design is then typically achieved from the linearized equations of motion along desired trajectories. However, some aerial vehicles are required to fly in conditions that involve large and rapid variations of the angle of attack. Examples are given by fighter aircraft, convertible aircraft, or small Unmanned Aerial vehicles (UAVs) operating in windy environments. In fact, some Vertical Take-Off and Landing (VTOL) vehicles, like e.g. ducted fans, are subjected to large variations of the angle of attack when transitioning from hover to horizontal flight. It then matters to ensure large stability domains that are achievable via the use of nonlinear feedback designs. Nonlinear feedback control of aircraft can be traced back to the early eighties. Following [START_REF] Singh | Output feedback nonlinear decoupled control synthesis and observer design for maneuvering aircraft[END_REF], control laws based on the dynamic inversion technique have been proposed to extend the flight envelope of military aircraft (see, e.g., [START_REF] Wang | Robust nonlinear flight control of highperformance aircraft[END_REF] and the references therein). The control design relies on tabulated models of aerodynamic forces and moments, like the High-Incidence Research Model (HIRM) [START_REF] Muir | Robust flight control design challenge problem formulation and manual: The high incidence research model (hirm)[END_REF]. Compared to linear techniques, this type of approach yields extended flight domains without involving gain scheduling strategies. The angle of attack is still assumed to remain in a limited range, however, away from the stall zone. Nonlinear feedback control of VTOL vehicles has been addressed with a larger variety of techniques, like dynamic inversion [START_REF] Hauser | Nonlinear control design for slightly non-minimum phase systems: Application to V/STOL[END_REF], Lyapunov-based design [START_REF] Isidori | Robust autonomous guidance: an internal-model based approach[END_REF], backstepping [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrotor[END_REF], sliding modes [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrotor[END_REF][START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF], or predictive control [START_REF] Kim | Nonlinear model predictive tracking control for rotorcraft-based unmanned aerial vehicles[END_REF] (See e.g. [START_REF] Hua | Introduction to feedback control of underactuated vtol vehicles[END_REF] for more references). Since most of these studies address the stabilization of hover flight or lowvelocity trajectories, aerodynamic effects are typically ignored, or modeled as a simple additive perturbation. In highly dynamic flight or harsh wind conditions, however, aerodynamic effects become important. This raises several questions such as, e.g., which aerodynamic models should be considered for the control design? and which feedback control solutions can be inferred from these models to ensure large stability domains and robustness? Methods used in aerodynamics to precisely describe aerodynamic forces, e.g. computational fluid dynamics or wind tunnel measurements, do not provide analytical expressions of aerodynamic characteristics (AC). From a control viewpoint they are useful to finely tune a controller around a given flight velocity, but exploiting them in the case of large flight envelopes (i.e., that involve strong variations of either the flight velocity or the angle of attack) is difficult. In this paper we advocate the use of simple analytical models of AC that can account for important structural properties of the system in a large flight envelope. The idea is then to exploit these properties at the control design level and rely on the robustness of feedback controllers to cope with discrepancies between the model and the true AC. More precisely, for (heavier than air) underactuated vehicles with a body-shape symmetric with respect to (w.r.t.) the thrust axis, we provide conditions on the AC under which the vehicle's equilibrium orientation associated with a desired flight velocity is explicitly (and uniquely) defined. We also show that such conditions are satisfied by simple models that approximate at the first order AC of real systems reported in the literature. The control design then essentially consists in aligning the thrust direction with the desired equilibrium orientation and monitoring the thrust intensity to compensate for the intensity of external forces. This thrust direction control paradigm has been exploited for underactuated VTOL vehicles either by neglecting aerodynamic effects [START_REF] Guenard | Dynamic modeling and intuitive control strategy for an "x4-flyer[END_REF], or by considering systems submitted to drag forces only [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF]. Determination of the vehicle's equilibrium orientation is straightforward in these cases but it is a major issue for more general vehicles [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF]. This paper shows that the thrust direction control paradigm can be extended to aerial vehicles submitted to significant lift forces.

The paper is organized as follows. Section 2 provides the notation and background. In Section 3, we show that for a class of symmetric bodies the dynamical equations of motion can be transformed into a simpler form. This transformation is used in Section 4 to propose a control design method applicable to several vehicles of interest.

2 Notation and background E 3 denotes the 3D Euclidean vector space and vectors in E 3 are denoted with bold letters. Inner and cross products in E 3 are denoted by the symbols • and × respectively. I = {O; i 0 , j 0 , k 0 } denotes a fixed inertial frame w.r.t. which the vehicle's absolute pose is measured (see Figure 1). This frame is chosen as the NED frame (North-East-Down). B = {G; i, j, k} denotes a frame attached to the body, with G the body's center of mass. The linear and angular velocities v and ω of the body frame B are then defined by

v := d dt OG , d dt (i, j, k) := ω × (i, j, k) , (1)
Throughout the paper, the time-derivative is taken w.r.t. the inertial frame I and, for any couple of points A, B, AB ∈ E 3 denotes the vector from A to B.

j 0 i 0 O k 0 j i G k T ω v Figure 1. Notation.

Equations of motion

Aircraft are modeled as rigid bodies of constant mass controlled via four control inputs: the thrust intensity T ∈ R of a body-fixed thrust force T = -T k and the three components of a control torque vector Γ G . This class of systems covers a large variety of vehicles like multi-copters, helicopters, convertibles UAVs, or even conventional airplanes. The vehicle's equations of motion are given by

m v = mg + F a -T k, (2a) d dt (i, j, k) = ω × (i, j, k), (2b)
d dt (J ω) = GP × F a + Γ G . (3)
with m the vehicle's mass, g = gk 0 the gravitational acceleration vector, (F a , P ) the resultant of the aerodynamic forces and its point of application, G the vehicle's center of mass, and J the inertia operator at G. Here we implicitly assume that so-called gyroscopic torques and body-forces are negligible or that they have already been compensated for via a preliminary torque control action (See, e.g., [START_REF] Pounds | Modelling and control of a large quadrotor robot[END_REF], [START_REF] Bristeau | The role of propeller aerodynamics in the model of a quadrotor UAV[END_REF], [START_REF] Pflimlin | Commande d'un minidrone à hélice carénée: de la stabilisation dans le vent à la navigation autonome[END_REF]Ch. 3] for more details).

Aerodynamic forces

The modelling of aerodynamic forces and torques F a and M a := GP ×F a acting on the vehicle is of particular importance, see, e.g. [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF] [17, Ch. 2] for fixed-wing aircraft, [START_REF] Huang | Aerodynamics and Control of Autonomous Quadrotor Helicopters in Aggressive Maneuvering[END_REF] [14] for quadrotors, [START_REF] Johnson | Modeling, control, and flight testing of a small ducted fan aircraft[END_REF] [20] [15, Ch. 3] for ductedfan tail-sitters, and [START_REF] Prouty | Helicopter Performance, Stability, and Control[END_REF] for helicopters. Denote by v a the air velocity, i.e. v a = vv w with v w the wind's velocity. The aerodynamic force is typically decomposed as F a = F L + F D where the lift force F L is the component perpendicular to v a , and the drag force F D is the component along v a 's direction. Consider any pair of angles (α, β) characterizing the orientation of v a with respect to the body frame (e.g. Figure 2). Knowing that the intensity of the steady aerodynamic force varies approximately like |v a | 2 , there exist two dimensionless positive functions C L and C D depending on the Reynolds number R e , the Mach number M, and (α, β), and such that

F L = k a |v a |C L (R e , M, α, β)r(α, β, |v a |) × v a , F D = -k a |v a |C D (R e , M, α, β)v a , r • v a = 0, |r| = 1, k a := ρΣ/2, (4)
with ρ the air density, Σ an area germane to the body shape, r(•) a unit vector-valued function. C D and C L are the aerodynamic characteristics of the body, i.e. the socalled drag and lift coefficients. Here we neglect the effects of the vehicle's rotational and unsteady motions on its surrounding airflow (see [17, p. 199] for more details).

Reduced control model

To develop control principles that apply to many aerial vehicles, one must get free of actuation specificities. In agreement with a large number of works on VTOL control and in view of Eq. ( 3), which points out how ω can be modified via the choice of the control torque Γ G , one can consider ω as an intermediate control input (see [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF] for details). The corresponding physical assumption is that "almost" any desired angular velocity can be obtained after a short transient time. The control model then reduces to Eqs. ( 2) with T and ω as control inputs.

3 Symmetric bodies and spherical equivalence From Eq. ( 2), for the body to move with a constant velocity the controlled thrust vector T k must be equal to the resultant external force F ext := mg + F a . When F a does not depend on the vehicle's orientation, as in the case of spherical bodies, F ext does not depend on this orientation either. The thrust direction at the equilibrium is then unique and it is explicitly given by the direction of F ext . The control strategy then basically consists in aligning the thrust direction k with the direction of F ext (using ω as control input) and in opposing the thrust magnitude to the intensity of F ext (using T as control input). This is the basic principle of the thrust direction control paradigm [START_REF] Guenard | Dynamic modeling and intuitive control strategy for an "x4-flyer[END_REF][START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF]. For most aerial vehicles, however, aerodynamic forces depend on the vehicle's orientation, and thus on the direction of k. The equilibrium relation T k = F ext then becomes an implicit equation with both sides of this equality depending on k. In this case, existence, uniqueness, and explicit determination of the equilibrium thrust direction(s) become fundamental questions for the control design [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF]. In this section, we provide answers to these questions for a class of axisymmetric vehicles. More precisely, let us consider vehicles whose external surface S is characterized by the existence of an orthonormal body frame B c = {G c ; i, j, k} that satisfies either one of the following assumptions: Here, rot Ov (ψ)(P ) denotes the rotation about the axis Ov by the angle ψ of the point P . Examples of "symmetric" and "bisymmetric" shapes satisfying these assumptions are represented in Figure 3 (with G = G c ). Various human-made aerial devices (rockets, missiles, airplanes with annular wings, etc.) satisfy Assumption 1 in the first approximation. When this assumption holds true, one can define α ∈ [0, π] as the angle of attack between -k and v a , and β ∈ (-π, π] as the angle between i and the projection of v a on the plane {G c ; i, j} (see Figure 2). Assumption 1 also implies the following properties: P1 : the aerodynamic force F a does not change when the body rotates about its axis of symmetry G c k; P2 : F a ∈ span{k, v a }. Property P1) in turn implies that the AC do not depend on β, whereas Property P2) implies that: i) the unit vector r in ( 4) is orthogonal to k and independent of α and, ii) the lift coefficient is equal to zero when α = {0, π}. Then, from Eq. ( 4),

F L = k a |v a |C L (R e , M, α)r(β) × v a , F D = -k a |v a |C D (R e , M, α)v a , r(β) = -sin(β)i + cos(β)j. (5)
The above choice of (α, β) also yields:

v a = |v a |(sin(α) cos(β)i + sin(α) sin(β)j -cos(α)k) so that, from (5) r(β) × v a = -cot(α)v a -|va| sin(α) k, F a = -k a |v a | C D (•) + C L (•) cot(α) v a + CL(•) sin(α) |v a |k . (6)
For constant Reynolds and Mach numbers the AC depend only on α and the following result follows from [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrotor[END_REF].

Proposition 1 Consider an axisymmetric thrustpropelled vehicle subjected to aerodynamic forces given by [START_REF] Isidori | Robust autonomous guidance: an internal-model based approach[END_REF]. Assume that

C D (α) + C L (α) cot(α) = C D0 , (7)
with C D0 constant. Then, Eq. (2) can be written as

m v = mg + F p -T p k, (8) 3 with
T p = T + k a |v a | 2 C L (α) sin(α) , F p (v a ) = -k a C D0 |v a |v a . ( 9
)
This proposition is the core result of this paper. It shows that an axisymmetric body subjected to both drag and lift forces can be seen as a sphere subjected to the orientation independent drag force F p and powered by the thrust force T p = -T p k. As shown further on, it then becomes possible to rely on existing techniques to design feedback stabilizers. In particular, it follows from (8) that given a desired reference velocity v r , there exists a unique (up to sign) equilibrium thrust direction k ref as long as |mg + F pm vr | = 0 along this reference velocity. More precisely,

k ref = mg + F p (v r,a ) -m vr |mg + F p (v r,a ) -m vr |
where v r,a = v rv w . This holds true provided that the AC satisfy [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF]. An example of simple functions C L , C D that satisfy [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF] and the π-periodicity property w.r.t. α associated with bisymmetric bodies is given next.

Proposition 2

The functions C D and C L defined by

C D (α) = c 0 + 2c 1 sin 2 (α), C L (α) = c 1 sin(2α), (10)
with c 0 and c 1 two real numbers, satisfy Condition [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF] with C D0 = c 0 + 2c 1 . Eq. ( 9) then yields:

T p = T + 2c 1 k a |v a | 2 cos(α), F p (v a ) = -k a C D0 |v a |v a .
The proof is straightforward.

A particular bisymmetric body is the sphere whose AC (zero lift coefficient and constant drag coefficient) are obtained by setting c 1 = 0 in [START_REF] Guenard | Dynamic modeling and intuitive control strategy for an "x4-flyer[END_REF]. Elliptic-shaped bodies are also bisymmetric but in contrast with the sphere they do generate lift. Approximation of measured AC with functions given by ( 10) is illustrated by the Figure 4 (left) with experimental data borrowed from [22, p.19] for an elliptic-shaped body (R e = 7.96 • 10 6 , M = 6).

The identified coefficients are c 0 = 0.43 and c 1 = 0.462. Since missile-like devices are "almost" bisymmetric, approximating their AC with such functions can also be attempted. For instance, the approximation shown in In both cases the approximation should be good enough for feedback control purposes.

Control design

The results of Section 3 are now exploited for the control design of axisymmetric vehicles. We start by revisiting the thrust direction control problem.

Measurements

Eq. ( 10)

Figure 4. AC of elliptic and missile-like bodies.

Thrust direction control

Consider a time-varying reference thrust (unitary) direction k r . It is assumed that k r varies smoothly with time so that kr (t) is well defined for all t. The following result provides control expressions for the angular velocity control input ω yielding a large stability domain.

Proposition 3

The feedback law

ω = k 1 (k, t) + γ(t) γ(t) k × k r + ω r + λ(k, t)k (11)
with ω r = k r × kr the instantaneous angular velocity of k r , λ(•) any real-valued continuous function, γ(•) any smooth positive real-valued function such that inf t γ(t) > 0, and k 1 (•) any continuous positive real-valued function such that inf k,t k 1 (k, t) > 0, ensures exponential stability of the equilibrium k = k r with domain of attraction {k(0

) : k(0) • k r (0) = -1}.
The proof is given in the appendix. The above expression of ω is a coordinate-free generalization of the solution proposed in [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF], for which the control gain γ was not present and a specific choice of k 1 was imposed. The first term in the right-hand side of ( 11) is a feedback term depending on the error between k and k r . The second term is a feedforward term. The last term is associated with the rotation about the axis k. It does not affect the thrust direction dynamics since k = ω × k. Let us comment on the choice of the control gains. Concerning λ(•), the simplest choice is λ(t) ≡ 0.

Another possibility is λ(t) = -ω r (t) • k(t). This yields ω(t)•k(t) = 0 so that the control does not induce instantaneous rotation about k. Concerning γ and k 1 , constant gains can be chosen but other choices can be preferred. For instance, k 1 (k, t) = k 1,0 /(1 + k • k r (t) + ǫ 1 ), with k 1,0 > 0 and ǫ 1 a small positive number, grows large when k gets close to -k r , thus making this undesired equilibrium more repulsive.

Velocity and position control

Let v r (•) denote a reference velocity time-function (at least three times differentiable everywhere). Velocity control consists in the asymptotic stabilization of the velocity error ṽ := vv r at zero. This control objective may be complemented by the convergence to zero of a position error p := pp r , with p r (•) a reference position time-function. In this case, v r = ṗr , and the error state vector to be stabilized at zero contains the six-dimensional vector ( p, ṽ). The error vector may further include an integral of the position error p. In order to address various control problems, we consider a "generalized" control objective consisting in the asymptotic stabilization at zero of an error vector denoted as ( ρ, ṽ), with ρ ∈ R p and such that ρ = f ( ρ, ṽ), with f (•, •) denoting a smooth vector-valued function. For instance, if ρ = p with p denoting either a position error or an integral of the velocity error ṽ, then f ( ρ, ṽ) = ṽ. If ρ = (I p , p), with I p denoting a saturated integral of the position tracking error such that d dt I p = h( ρ), then f ( ρ, ṽ) = (h( ρ), ṽ). The simplest case is pure velocity control without integral correction, i.e. ρ = ∅. Consider now an axisymmetric vehicle. From [START_REF] Kim | Nonlinear model predictive tracking control for rotorcraft-based unmanned aerial vehicles[END_REF],

m v = F p + m(g -a r ) -T p k = mξ + Fp -T p k, (12)
with a r := vr the reference acceleration and

Fp (v a , a r , ξ) := F p (v a ) + m(g -a r -ξ). ( 13
)
The idea is to end up working with the simple control system v = ξ where ξ plays the role of a virtual feedback term. To this aim Eq. ( 12) suggests to make Fp -T p k converge to zero. With T p preferred positive, this implies that the thrust direction k should tend to

k r := Fp /| Fp | . ( 14
)
Recall from (9) that F p does not depend on k. Thus, provided that ξ does not depend on k, Fp does not depend on k either, and k r is well defined as long as Fp does not vanish. This is precisely what makes Proposition 1 important for the control design. Convergence of Fp -T p k to zero also implies that T p must tend to Fp • k.

From ( 2) and ( 8), this is equivalent to the convergence of the thrust intensity T to Fa • k with

Fa := F a + m(g -a r -ξ) (15)
Once the reference thrust direction k r is well defined, a stabilizing feedback law can be designed:

Proposition 4 Consider an axisymmetric vehicle with aerodynamic characteristics satisfying relation [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF], and a smooth feedback controller ξ( ρ, ṽ) for the control system ρ = f ( ρ, ṽ)

(16a) v = ξ (16b)
Assume that A1 : ξ( ρ, ṽ) makes ( ρ, ṽ) = (0, 0) a locally exponentially stable equilibrium point of System (16); A2 : Fp does not vanish along the velocity reference trajectory v r , i.e., ∃δ > 0 :

δ ≤ | Fp (v r,a (t), a r (t), 0)|, ∀t, with v r,a := v r -v w .
Then, T = Fa • k and ω given by [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF], with k r defined by [START_REF] Bristeau | The role of propeller aerodynamics in the model of a quadrotor UAV[END_REF], γ = c 2 + | Fp | 2 , and c 2 any strictly positive constant, ensure local exponential stability of the equilibrium point ( ρ, v, k) = (0, v r , k r ) for the system (16a)-( 2).

The proof is given in the appendix. Let us comment on the above result. 1) Proposition 4 shows how to derive an exponential stabilizer for the underactuated System (2) from an exponential stabilizer for the fully-actuated system v = ξ, e.g., a linear feedback controller.

2) Given an exponential stabilizer ξ for System [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF], local exponential stability of zero tracking errors for an axisymmetric vehicle for which the AC satisfy relation [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF] relies on Assumption A2, which imposes that the reference thrust direction k ref , i.e. along the reference trajectory, is well defined at all times. This condition may be violated for very specific and aggressive reference trajectories. Note, however, that its satisfaction can be checked from the knowledge of the reference velocity only (provided that an accurate model of AC is available).

3) Feedback control calculation calls for a few remarks. A difficulty is that both k r and γ depend on Fp . Since γ and kr are involved in the calculation of ω, the timederivative of Fp has to be calculated also. In this respect, in view of ( 8) and the way T p is calculated, note that v and subsequently Ḟp do not depend on ω. This guarantees the well-posedness of the control calculation provided only that k r is well defined. In practice, it is also possible to estimate Ḟp using the calculation of Fp and a high-gain observer. Another possibility is to use the reference velocity instead of the vehicle's actual velocity to calculate an approximation of this term. Proposition 4 guarantees local asymptotic stability only. The difficulty to ensure a large domain of convergence comes from the risk of Fp vanishing at some point, which would in turn make k r , as specified by [START_REF] Bristeau | The role of propeller aerodynamics in the model of a quadrotor UAV[END_REF], ill-defined. This risk, although small, cannot be ruled out, especially because the term F p in Fp can be very large due to aerodynamic forces. In practice, having a control always well defined requires to modify the term Fp used in the control expression so that it cannot vanish. This is a subject of future studies. Assuming that Fp remains different from zero, then convergence of the tracking errors can be guaranteed, as specified by the following proposition the proof of which follows from that of Proposition 4.

Proposition 5 Consider the feedback law of Proposition 4. Assume further that ξ( ρ, ṽ) globally asymptotically stabilizes the origin ( ρ, ṽ) = (0, 0) of the system

ρ = f ( ρ, ṽ) v = ξ( ρ, ṽ) + ε(t) (17)
when ε(•) is identically zero, and still ensures the convergence to zero of the solutions to this system when ε(•) converges to zero exponentially. Then, any solution to the closed-loop system (16a)-( 2) along which Fp does not vanish (in the sense that ∃δ > 0 : δ ≤ | Fp (v a (t), a r (t), ξ( ρ, ṽ))|, ∀t) converges to the equilibrium point (0, v r , k r ).

Simulation results

We consider a missile-like body with AC given by the measurements of Figure 4 (right). The control objective is the asymptotic stabilization at zero of the velocity error ṽ. A saturated integral ρ = I v of this error is used in the control law in order to compensate for static modelling errors and additive perturbations. I v is obtained as the (numerical) solution to the following equation [START_REF] Hua | Time sub-optimal nonlinear pi and pid controllers applied to longitudinal headway car control[END_REF]:

İv = f ( ρ, ṽ) = -k I I v + k I sat δ I v + ṽ k I ; I v (0) = 0,
with k I a positive gain, δ > 0 the upperbound of |I v |, and sat δ a differentiable approximation of the classical saturation function defined by sat δ = min 1, δ |x| x. The feedback law of Proposition 4 is applied with ξ( ρ, ṽ)

= -k v ṽ -k i I v , k 1 (k, t) = k 1,0 /(1 + k • k r (t) + ǫ 1 ) 2 , and λ(k, t) = -ω r (t)•k(t), where k v = 5, k i = k 2 v /4
, k I = 50, k 1,0 = 10, and ǫ 1 = 0.01. feedforward term ω r is evaluated using the reference acceleration vr than the vehicle's acceleration v calculated from Newton's equation [START_REF] Wang | Robust nonlinear flight control of highperformance aircraft[END_REF] and the model of aerodynamic forces F a . The simulated vehicle's equations of motion are given by ( 2)-( 5), with the AC C L and C D obtained by interpolating the measurements reported in [23, p.54] (see Figure 4 (right)). These coefficients thus differ from the approximating functions [START_REF] Guenard | Dynamic modeling and intuitive control strategy for an "x4-flyer[END_REF] used in the control calculation, here computed by using the identified coefficients c 0 = 0.1 and c 1 = 11.55 reported previously. The missile's physical parameters are m = 100 [Kg] and (ρ, Σ) = (1.292, 0.5) ( Kg/m 3 , [m] 2 ), so that k a ≈ 0.3 [Kg/m]. These values are replaced by estimated ones, namely ka = 0.24 and m = 80 [Kg], in the control calculation in order to test the control robustness w.r.t. parametric uncertainties. In particular, the vector Fp in ( 13) is calculated with F p (v a ) = -ka (c 0 + 2c 1 )|v a |v a . The reference velocity expressed in Mach numbers is:

v r (t) =        0.7i 0 0 ≤t < 10, -0.7j 0 10 ≤t < 20, -0.7k 0 20 ≤t < 30, -0.7i 0 30 ≤t < 40, (18)
and v r (t) = -0.5 sin(tπ/5)i 0 + 0.6 sin(tπ/10)j 0 + 0.6 cos(tπ/10)k 0 when 40 ≤ t < 60. The applied thrust force and angular velocity ω = (i, j, k)ω are saturated as follows: 0 < T < 10 mg, |ω i | < 2π, i = {1, 2, 3}. The initial velocity and attitude are: v(0) = 0.5i 0 [M ach], φ 0 = ψ 0 = 0 • , θ 0 = -40 • where (φ, θ, ψ) denote standard roll, pitch, and yaw angles as defined in [17, p. 47].

From top to bottom, Figure 5 shows the evolution of the reference velocity v r = (i 0 , j 0 , k 0 ) ẋr , the vehicle's velocity v = (i 0 , j 0 , k 0 ) ẋ, angle of attack α, the angular velocity ω = (ı, , k)ω, the applied thrust-to-weight ratio, the norm of Fp (which must be different from zero to ensure the well-posedness of the control solution), and the angle θ between k and k r . There is no wind. The initial angle of attack at t = 0 is 50 • The attitude control makes this angle decrease rapidly. Sharp discontinuities of ẋr at t = 10, 20, 30, 40 are responsible for the observable transitions and temporarily large angles of attack. Thanks to the integral correction terms, ẋ converges to ẋr when ẋr is constant. For 40 ≤ t < 60, although ẋr varies rapidly velocity errors remain small, thanks to both feedforward and integral correction terms. Figure 6 illustrates the improvement brought by the control proposed in this paper w.r.t. a control that does not take the dependence of the aerodynamic forces upon the vehicle's orientation into account. To this aim, we consider the velocity control proposed in [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF] for spherical-like vehicles subjected to aerodynamic drag solely. This control is basically the same as the one considered in Proposition 4 with Fa used in place of Fp in the control law. Figure 6 shows the evolution of | Fa | and θ when applying this control, with the feedforward term ω r (whose calculation involves Ḟa ) set equal to zero for the sake of simplification. One can observe from this figure that i) relative variations of the norm of | Fa | are significantly more pronounced than those of | Fp | in Figure 5, ii) the amplitude of the orientation error θ after a discontinuous change of the reference velocity is much more important, and even more significantly iii) Fa crosses zero little after the reference velocity discontinuity occurring at t = 40. The reference direction k r and thus the control law, are not defined at this point, thus leading to an abrupt stop of the simulation.

Conclusion and perspectives

Extension of the thrust direction control paradigm to a class of vehicles with axisymmetric body shapes has been addressed. Application examples include, e.g., rockets and aerial vehicles with annular wings. Specific aerodynamic properties associated with these particular shapes allow for the design of nonlinear feedback controllers yielding asymptotic stability in a very large flight envelope. Simulation results for bisymmetric and quasibisymmetric bodies validate the proposed control approach. Further extension of the present approach to vehicles with non-symmetric body shapes (e.g. conventional airplanes) is currently investigated.

Appendix

We make use of the following classical vectorial relations:

∀x, y, z ∈ E 3 , x • (y × z) = y • (z × x) ∀x, y, z ∈ E 3 , x × (y × z) = (x • z)y -(x • y)z (19) Proof of Proposition 3: Consider the function V 0 := 1 -k • k r ,
which is non-negative and vanishes only when k = k r since |k| = |k r | = 1. Recall that ω r = k r × kr . From [START_REF] Johnson | Modeling, control, and flight testing of a small ducted fan aircraft[END_REF], kr = ω r × k r so that

V0 = -k • k r -k • kr = -(ω × k) • k r -k • (ω r × k r ) = (ω r -ω) • (k × k r ) (20)
where the last equality follows from [START_REF] Johnson | Modeling, control, and flight testing of a small ducted fan aircraft[END_REF]. Now, define

V 1 := γ 2 (t) 2 1 -k • k r 1 + k • k r = γ 2 (t) 2 1 -(k • k r ) 2 (1 + k • k r ) 2 = γ 2 (t) 2 |k × k r | 2 (1 + k • k r ) 2 (21)
where the last equality comes from that (k

• k r ) 2 + |k × k r | 2 = 1, since |k| = |k r | = 1. One verifies that V1 = γ(t) γ(t) 1 -k • k r 1 + k • k r + γ 2 (t) V0 (1 + k • k r ) 2
and it follows from ( 20) and ( 21) that

V1 = γ(t) γ(t) |k × k r | 2 (1 + k • k r ) 2 + γ 2 (t) (k × k r ) • (ω r -ω) (1 + k • k r ) 2 = γ(t) k × k r (1 + k • k r ) 2 γ(t) k × k r + γ(t)(ω r -ω)
Replacing ω by its expression [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF] yields [START_REF] Keyes | Aerodynamic characteristics of lenticular and elliptic shaped configurations at a mach number of 6[END_REF] Since inf k,t k 1 (k, t) > 0, V 1 converges exponentially to zero. Exponential stability of k = k r then follows from the definition of V 1 and the fact that inf t γ(t) > 0.

V1 = -k 1 (k, t)γ 2 (t) |k × k r | 2 (1 + k • k r ) 2 = -2k 1 (k, t)V 1
Proof of Proposition 4: In view of Assumption 2, k r is well defined in a neighborhood of the equilibrium point ( ρ, v, k) = (0, v r , k r ). Note also that γ(t) ≥ √ c 2 > 0.

Therefore, the feedback law is well defined in a neighborhood of the equilibrium point. From ( 2) and ( 8

where the last equality comes from [START_REF] Johnson | Modeling, control, and flight testing of a small ducted fan aircraft[END_REF]. Therefore, along the solutions to the controlled system, the variables ρ and ṽ satisfy Eq. ( 17) with ε := 1 m | Fp |(k × (k r × k)). Since k and k r are unit vectors, it follows from [START_REF] Prouty | Helicopter Performance, Stability, and Control[END_REF] and the definition of γ that

V 1 ≥ γ 2 (t) 8 |k × k r | 2 ≥ | Fp | 2 8 |k × k r | 2 ≥ | Fp | 2 8 |k × (k × k r )| 2 ≥ m 2 |ε| 2 8 (24)
From the proof of Proposition 3, V 1 converges to zero exponentially so that ε also converges to zero exponentially. From Assumption 1 and converse Lyapunov theorems (See, e.g., [START_REF] Khalil | Nonlinear systems[END_REF]Section 4.7]) there exists a quadratic Lyapunov function V 2 ( ρ, ṽ) for System [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF], i.e., such that in a neighborhood of ( ρ, ṽ) = (0, 0), V2 ( ρ, ṽ) ≤ -k 2 V 2 ( ρ, ṽ)

Using the triangular inequality, it follows from ( 22), ( 17), [START_REF] Hua | Time sub-optimal nonlinear pi and pid controllers applied to longitudinal headway car control[END_REF], and (25) that V = αV 1 +V 2 is a Lyapunov function for the controlled system for α > 0 large enough.

Assumption 1 (Figure 2 .Figure 3 .

 123 Figure 2. The (α, β) angles.

Fig- ure 4

 4 (right) uses experimental data taken from [23, p.54] for a missile moving at M = 0.7. In this case, c 0 = 0.1 and c 1 = 11.55.

Figure 5 .

 5 Figure 5. C-701 simulation with Fp in angular control.

Figure 6 .

 6 Figure 6. C-701 simulation with Fa in angular control.

  ), F p -T p k = F a -T k. Therefore Fp -T p k = Fa -T k and Fp • k = T p -T + Fa • k = T p .From this relation and (14), Eq. (12) can be written as:m v = | Fp |k r -T p k + mξ( ρ, ṽ) = | Fp |k r -( Fp • k)k + mξ( ρ, ṽ) = | Fp |k r -(| Fp |k r • k)k + mξ( ρ, ṽ) = | Fp |(k r -(k r • k)k) + mξ( ρ, ṽ) = | Fp |(k × (k r × k)) + mξ( ρ, ṽ)

Preprint submitted to AutomaticaDecember 1, 2014

This work was supported by the "Chaire d'excellence en Robotique RTE-UPMC" and by the French ANR Project SCAR.


Florent Le Bras

Tarek Hamel

Robert Mahony

Claude Samson

Observer design for position and velocity bias estimation from a single direction output

This paper considers the questions of observability, and design of an observer, for position estimate of an object moving in R n with direction measurements, typically a bearing direction in three dimensions derived from a vision system. We provide a comprehensive observability analysis and discuss stability of the observer error dynamics. We show that even with bias in the velocity measurement it is possible to recover uniform global exponential stability of error dynamics. Simulation results on demonstrate the performance of the proposed approach.

I. INTRODUCTION

There is is rich literature in vision based pose estimation driven by advances in the structure from motion in the field of computer vision [START_REF] Häming | The structure-from-motion reconstruction pipeline a survey with focus on short image sequences[END_REF]. Most of the recent structure from motion algorithms are formulated as an optimisation problem over a set of selected images [START_REF] Triggs | Bundle adjustment a modern synthesis[END_REF], however, recent work has emphasised the importance of considering motion models and filtering techniques [START_REF] Strasdat | Visual slam: Why filter[END_REF] for a class of important problems. Recursive filtering methods for vision based structure from motion and pose estimation themselves have a rich history primarily associated with stochastic filter design such as EKF, unscented filters and particle filters [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF], [START_REF] Soatto | Motion estimation via dynamic vision[END_REF], [START_REF] Armesto | Fast ego-motion estimation with multi-rate fusion of intertial and vision[END_REF], [START_REF] Civera | Inverse depth parametrization for monocular slam[END_REF]. A comparison of EKF and particle filter algorithms for vision based SLAM is available in [START_REF] Bekris | Evaluation of algorithms for bearing-only slam[END_REF]. Although nonlinear observer design does not provide a stochastic interpretation of the state estimate they hold the promise to handle the non-linearity of the vision pose estimation problem in a robust and natural manner [START_REF] Baldwin | A nonlinear observer for 6 DOF pose estimation from inertial and bearing measurements[END_REF]. Ghosh and subsequent authors consider nonlinear observers on the class of perspective systems [START_REF] Rehbinder | Pose estimation using line-based dynamic vision and inertial sensors[END_REF], [START_REF] Abdursul | Nonlinear observers for perspective time-invariant linear systems[END_REF], [START_REF] Abdursul | Nonlinear observers for perspective time-invariant linear systems[END_REF], [START_REF] Aguiar | Minimum-energy state estimation for systems with perspective outputs[END_REF], [START_REF] Dahl | Observer forms for perspective systems[END_REF], that is systems with output in a projective space obtained as a quotient of the state space. Perspective outputs y(x) are of the form

y P = ( x 1 x n , . . . x n-1 x n , 1)
and correspond to the nonlinear projection along rays through the origin onto an affine image plane perpendicular to the focal axis. The output representation is attractive in that it corresponds to the normal representation of vision data for perspective cameras. Indeed, there are a number of works that consider filtering for y P directly, rather than estimating the camera position [START_REF] Dixon | Range identification for perspective vision systems[END_REF], [START_REF] De Luca | On-line estimation of feature depth for image-based visual servoing schemes[END_REF], [START_REF] Dani | Globally exponentially stable observer for vision-based range estimation[END_REF], corresponding to image tracking. Although significant work has been based on this output representation, it tends to lead to complex observer and filter design and difficult analysis [START_REF] Rehbinder | Pose estimation using line-based dynamic vision and inertial sensors[END_REF], [START_REF] Abdursul | Nonlinear observers for perspective time-invariant linear systems[END_REF], [START_REF] Abdursul | Nonlinear observers for perspective time-invariant linear systems[END_REF], [START_REF] Aguiar | Minimum-energy state estimation for systems with perspective outputs[END_REF], [START_REF] Dahl | Observer forms for perspective systems[END_REF]. An additional question of importance concerns the rate of convergence of an observer and recent work has addressed this question in the context of controlling the camera motion to improve observability of the problem and increase the rate of convergence of the observer [START_REF] Spica | Active structure from motion: Application to point, sphere, and cylinder[END_REF].

The present paper contributes further to the field of nonlinear observer design for systems with direction outputs. The key contribution that we make is the development of an elegant and rigorous stability analysis for a simple filter design. The filter is designed for a single bearing measurement and relies on the motion of the camera to generate persistence of excitation of the innovation in order to guarantee global asymptotic convergence. Rather than using the perspective outputs favoured in previous papers we use direction outputs y = x/|x| = y P |y P | corresponding to projection onto a virtual spherical image plane and differing from perspective outputs only in the scaling. The two formulations are essentially equivalent from a systems perspective in the region where perspective outputs are defined, however, we believe that the direction output representation contributes to the simplicity of the observer proposed in the present paper. We characterise the rate of convergence of the filter in terms of the persistence of excitation property. We then consider the case when the measurement of velocity of the camera is perturbed by an unknown bias. To the authors knowledge, this problem has not been considered in the nonlinear observer literature. We provide a rigorous proof of the global asymptotic stability of the observer state for this case by exploiting a novel state transformation. The simulations provided demonstrate the performance of the filter.

The paper is organised along five sections. Following the present introduction, Section II introduces the system under consideration and points out observability properties attached to it. Section III develops the main results of the paper. Section IV present a few illustrative simulations. Concluding remarks are provided in Section V.

II. PROBLEM DESCRIPTION

The system considered is the kinematics of an object moving in R n ẋ = v + a

(1)

y = x |x| ∈ S n-1 (2)
where v ∈ R 3 is the velocity of the object and a ∈ R 3 represents any unknown bias. Let S n-1 denote the unit sphere, the space of measurements y ∈ R n such that 1 |y| = 1. An example of such a measurement is the bearing in S 2 obtained from a camera looking at a moving object.

In most applications the unknown velocity a ∈ R n (with n = 3) represents the velocity of the fluid in which evolves the moving object or/and any bias that affects the measurement of v. In this paper, for the sake of generality, we consider an arbitrary dimension n ≥ 2. We emphasize that the value of y ∈ S n-1 and v ∈ R n must be known at all times.

A. Observability analysis

We first give a general observability criterion. The following persistency of excitation condition will then yield an observability result for system (1-2).

Definition 2.1: The direction y ∈ S n-1 , is called persistently exciting if there exist δ > 0 and 0 < µ < δ such that for all t t+δ t π y(τ ) dτ ≥ µI, with π y = (I -yy )

For future use, note that (3) is equivalent to

∀b ∈ S n-1 : t+δ t |π y(τ ) b| 2 dτ ≥ µ, (4)
Another characterization of persistent of excitation, in terms of the property that the time-derivative of ẏ must satisfy, is pointed out in the following lemma Lemma 2.2: Assume that ẏ(t) is uniformly continuous, then relation (3) (respectively (4)) is equivalent to:

∃ε > 0, ∃τ ∈ [t, t + δ] such that | ẏ(τ )| ≥ ε (5)
Proof: The proof of this lemma is given in Appendix Recall that two different points

x 0 1 , x 0 2 ∈ R n are said distinguishable, if there exists an input v(t) ∈ R n and a time t 1 such that for solutions x 1 (t), x 2 (t) of (1) with x 1 (0) = x 0 1 , x 2 (0) = x 0 2 we have y(x 1 (t 1 )) = y(x 2 (t 1 )
). Equivalently, in this case one says that the admissible input distinguishes the two initial states, and also that two initial states of system (1-2) are indistinguishable if they are not distinguished by any admissible input.

Definition 2.3: A system is called strongly observable if all pairs of distinct initial states are distinguished by all admissible inputs. It is called weakly observable if every pair of distinct initial states is distinguished by at least one admissible input.

Reasons to differentiate between strong observability and weak observability are well explained in the non-linear control literature. For complementary details on this subject we refer the reader to a classical work by Sussmann [START_REF] Sussmann | Single-input observability of continuous-time systems[END_REF].

Lemma 2.4: The system (1-2) complemented with the equation ȧ = 0, with X = x a as the system state vector, v as the system input, and y as the system output, is weakly observable but not strongly observable.

Proof:

Choose, for instance, the input v(t) = (cos(t), sin(t), 0 . . . 0) T . The solutions to the system are then given by x(t) = x(0) + (sin(t) + at, -cos(t) + at, 0 . . . 0) T and one easily verifies that y 1 (t) = y 2 (t), ∀t, implies that x 1 (0) = x 2 (0). This establishes the weak observability property of the system. Note also that the chosen input renders both outputs y 1 (t) and y 2 (t) persistently exciting in the sense of the definition (2.1). On the other hand, one verifies that, if the input v is constant, then initial states x 1 (0) = k 1 v and x 2 (0) = k 2 v, with k 1 and k 2 denoting arbitrary positive numbers, can not be distinguished because y 1 and y 2 are constant and equal in this case. This proves that the system is not strongly observable. The weak observability property of the system justifies the introduction of the persistence condition evoked previously to characterize "good" outputs (produced by "good" inputs) yielding a property of "uniform" observability that renders the state-observation problem addressed in the next section well-posed.

III. OBSERVER DESIGN

The problem of state observation refers to the design of an algorithm that allows one to recover actual state values from the observation of previous outputs. We start by the observer design for the classical situation addressed in the literature where the unknown constant velocity bias a is equal to zero. The situation when this term is different from zero and unknown a priori is addressed subsequently.

Lemma 3.1: Consider the system (1-2) and the following observer:

ẋ1 = v -kπ y x1 , x1 (0) = x0 1 ∈ R n and k > 0 (6)
Assume that a ≡ 0, x is bounded and never crosses zero, so that the output y is always well defined. Let x1 = x -x1 ∈ R n denote the estimation error. If v(t) ∈ R n is bounded and such that the measured direction y(t) is persistently exciting, then the equilibrium x1 = 0 is Uniformly Globally Exponentially Stable (UGES).

Proof: Differentiating x and using (1) and ( 6) one gets the following linear time varying system:

ẋ1 = kπ y x = -kπ y x1 ,
Using the assumption of persistent excitation characterized by relation (3) a direct application of [20, Lemma 5] proves that x1 is UGES. More explicitly, one verifies that the transition matrix Φ associated with the above system satisfies

exp -k(t-τ ) ≤ ||Φ(t, τ )|| ≤ exp -γ(t-τ ) , (7)
where γ = µk δ(1+k 2 δ) 2 . The interest of this result lies in the extreme simplicity of the observer design and, more importantly, in the property of global stability and explicit bounds on the convergence rate of the observer.

Lemma 3.2: Consider the system (1-2) and the above filter [START_REF] Armesto | Fast ego-motion estimation with multi-rate fusion of intertial and vision[END_REF]. If

• v ∈ R n is bounded and such that the measured direction y is persistently exciting, • x is bounded and never crosses zero, and • a ∈ R n is constant then |x| (and hence |x 1 |) is uniformly bounded w.r.t. initial conditions and ultimately bounded by 1 γ |a|. Proof: It is straightforward to verify that, in this case, the error-system equation is:

ẋ1 = a -kπ y x1 , (8)
whose general solution is:

x1 (t, 0) = Φ(t, 0)x 0 1 + t 0 Φ(τ, t)adτ Using (7) it follows that |x 1 (t)| ≤ (|x 0 1 | + 1 γ |a|) and lim t→+∞ |x 1 (t)| ≤ 1 γ |a|.
Since x is bounded by definition, it follows that x1 is also bounded.

For the design of an exponentially stable observer in the case where a = 0 the following two technical lemmas are instrumental.

Lemma 3.3: Assume that y ∈ S n-1 is persistently exciting. The matrix-valued function M (t) solution to:

Ṁ = I -kπ y M, M (0) = M (0) T = M 0 > 0 (9)
is bounded and always invertible, and its condition number is bounded.

Proof: See appendix B. Lemma 3.4: Assume that y ∈ S n-1 is persistently exciting and ẏ is uniformly continous. The dual output y := M -1 y |M -1 y| is also persistently exciting.

Proof: See appendix C. The observer design presented hereafter is based on the association of the filter (6) that ensures, as we will show, that the variable z := x1 + M a converges to x, with a second filter that provides an estimate x of x := M -1 z. It then suffices to pre-multiply this second estimate by M to obtain an estimate of x. The following theorem specifies the observer design and its convergence properties in the case where the output y is persistently exciting.

Theorem 3.5: Consider the system (1-2) along with (6) and [START_REF] Baldwin | A nonlinear observer for 6 DOF pose estimation from inertial and bearing measurements[END_REF]. Define the virtual observer z as follows

z = x1 + M a (10)
and the dual observer x of M -1 z as follows

ẋ = v -k π y x , x (0) = x 0 (11)
with v := M -1 v -M -1 x1 a known term and k any positive gain. If y is persistently exciting in the sense of Lemma 2.2 then the virtual error xz := x -z, the dual position error x := x -x , the position error x = x-M x , and the adaptation error â -a (with â := x -M -1 x1 ), globally exponentially converge to zero. Proof: The proof proceeds step by step. Concerning the convergence of xz to zero, one easily verifies, using (6-9), that:

ż = v + a + kπ y z (12)
Differentiating xz , and using ( 1) and ( 12), one obtains:

ẋz = kπ y z = -kπ y xz
This equation being the same as the one for x1 in the case where a = 0, one concludes as in Lemma 3.1 that xz = 0 is uniformly globally exponentially stable, provided that y is persistently exciting. Concerning the convergence of x to zero, using the fact that the matrix M is invertible (from Lemma 3.3), we introduce the virtual estimate z ∈ R n , defined as follows:

z := M -1 z = M -1 x1 + a (13)
Since we already have established that (x -z) converges exponentially to zero, and since M -1 is bounded, we deduce that (x -z ) converges exponentially to zero. Differentiating (13), and using ( 6) and ( 9), one verifies that:

ż = v (14)
Now, differentiating z = z -x , and using ( 11) and ( 14), it comes that:

ż = k π y x = -k π y (z -z )
Using the fact that x ∈ ker(π y ), and defining x z = x -z one gets:

ż = -k π y z -k π y x z = -k π y z -k π y M -1 xz
Since y is a persistently exciting (from Lemma 3.4), the above equation is similar to the one of x1 in the case where a = 0, except for the additive "perturbation" term -k π y M -1 xz which converges exponentially to zero, due to the exponential convergence of xz to zero. It is immediate to show that this exponentially vanishing perturbation does not prevent z from globally converging to zero exponentially. Therefore, x = x -x , and subsequently x = x -M x , also globally converge exponentially to zero. Finally, since z converges to zero so that x -z also converges to zero, with z = M -1 x1 + a, it comes that a -(x -M -1 x1 ) = a -â converges exponentially to zero.

IV. SIMULATION

We consider the example of a moving target point observed by a camera. The point moves in the 3D space (n = 3) along a circular trajectory at a fixed altitude (z = 3m) above the ground. The frame associated to the camera is located at the origin of the inertial frame whose optical axis is aligned with the z-axis and looks up at the moving point. The measure y ∈ S 2 corresponds to the spherical projection of the point, given by the algebraic transformation y = y P |y P | , where y P is the projective measure provided by the camera. The measurement of the velocity v is biased by a = (0.33, 0.66, 0.99) , and v is chosen so that v + a = (-0.5 sin 0.5t, 0.5 cos 0.5t, 0) . The following values of the observer gains are used: k = 0.5 and k = 5. In Figures [START_REF] Häming | The structure-from-motion reconstruction pipeline a survey with focus on short image sequences[END_REF][START_REF] Triggs | Bundle adjustment a modern synthesis[END_REF][START_REF] Strasdat | Visual slam: Why filter[END_REF] the performance of the observer in the ideal noise-free case is shown. From these figures one can observe the exponential convergence of all estimation errors to zero. In figures [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF][START_REF] Soatto | Motion estimation via dynamic vision[END_REF][START_REF] Armesto | Fast ego-motion estimation with multi-rate fusion of intertial and vision[END_REF], the observer algorithm is simulated in the case where the 3D bearing measurement y is calculated from the position x to which a uniform noise w taking values in the interval [-0.5m, 0.5m] is added. Figures [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF][START_REF] Soatto | Motion estimation via dynamic vision[END_REF][START_REF] Armesto | Fast ego-motion estimation with multi-rate fusion of intertial and vision[END_REF] show that the high frequency part of the noise is filtered by the proposed algorithm so that the performance of the proposed observer is not much reduced. In this paper, we discussed the issue of observability of a moving object in R n from bearing measurement and proposed nonlinear observer and a detailed analysis of this observer in the case where the bearing measurement satisfies a condition of persistent excitation. There is an increasing number of emerging applications that can make use of such an observer. We think, in particular, of applications involving cameras for relative localization of mobile robot teams. We believe that extending the observer design methodology described in the paper to the estimation of the relative . . .

ẏn (s n )    (τ -t)| 2 ≥ |b T (y(t)| 2 -2nd(δ)(τ -t) -4n 2 d(δ) 2 (τ -t) 2
for some

s i ∈ [t, τ ] (i = 1, . . . , n) and d(δ) = sup τ ∈[t,t+δ] | ẏ(τ )|. Choose b = y(t) so that b T y(t) = 1, then |b T y(τ )| 2 ≥ 1 -2nd(δ)(τ -t) -4n 2 d(δ) 2 (τ -t) 2 and t+δ t |b T y(τ )| 2 dτ ≥ δ -nd(δ)δ 2 - 4 3 n 2 d(δ) 2 δ 3
Clearly there exists > 0 (independent of t) such that We now show that (5) implies (4) One has |b T y(t)| 2 = cos(θ(t)) 2 with θ(t) the angle between b and y(t). Using the (assumed) uniform continuity of ẏ, (5) implies the existence of an interval [t 1 , t 2 ] ⊂ [t, t + δ] such that (t 2 -t 1 ) = 1 ( ) > 0 and |y(t 2 ) -y(t 1 )| ≥ 2 ( ) := 2n 1 ( ) > 0. This in turn implies that the distance between θ(t 1 ) and θ(t 2 ) is larger than some 3 ( ) > 0. Since the uniform continuity of ẏ implies the uniform continuity of θ, and since cos(θ(t 1 )) and cos(θ(t 2 )) cannot be both equal to one, one deduces the existence of a positive number 4 ( ) such that

t2 t1 cos(θ(τ )) 2 dτ ≤ (t 2 -t 1 )(1 -4 ( )). Therefore, t+δ t |b T y(τ )| 2 dτ = t1 t |b T y(τ )| 2 dτ + t2 t1 |b T y(τ )| 2 dτ + t+δ t2 |b T y(τ )| 2 dτ ≤ (t 1 -t) + (t 2 -t 1 )(1 -4 ( )) +(t + δ -t 2 ) ≤ δ -µ( ) with µ( ) = 1 ( ) 4 ( ) > 0.

B. Proof of Lemma 3.3

To prove that M is bounded, is suffices to ensure that, for any constant vector b ∈ S n-1 , |M b| is bounded. Define u = M b, it follows that:

u = b -kπ y u
This equation is similar to the equation (8) of x1 . Therefore

|M (t)b| ≤ |M (0)b| + 1 γ , ∀b ∈ S n-1
This implies that ||M || is bounded. To show that M is an invertible matrix, define ∆ := det(M ). From Jacobi's formula, one has

∆ = ∆tr(M -1 Ṁ ) = ∆tr(M -1 -kM -1 π y M ) = ∆tr(M -1 ) -k∆tr(π y M M -1 ) = -k(n -1)∆ + ∆tr(M -1 ) (15)
Note that this equation holds even if M is not invertible. Indeed, using the fact that det(M ) = n i=1 λ i and tr(M

-1 ) = n i=1 1 λi , with λ i (i = 1 . . . n) the eigenvalues of M , one verifies that ∆ = -k(n -1)∆ + n i=1 n j=1,j =i λ j (16)
Since M (0) is symmetric positive definite by assumption, all eigenvalues of M (0) are positive and ∆(0) > 0. Assume that ∆ is equal to zero for the first time at the time instant t 0 > 0. Then, tr(M (t)) > 0 on [0, t 0 [ and tr(M (t)) ≥ 0.

In view of ( 16), ∆(t) ≥ r(t) with r(t) the solution to the equation ṙ = -k(n -1)r, with r(0) = ∆(0). Therefore ∆(t) ≥ y(0) exp(-k(n -1)t) > 0, ∀t. This contradicts the existence of t 0 and proves that M (t) is always invertible.

Let us now prove that ∆(t) is lower bounded by a positive number. Rewrite equation ( 15) as follows

∆ = -tr(M -1 ) -(n -1)k ∆, (17)
Using the fact that tr(M -1 ) > n ∆ 1/n , this equation shows that ∆ ≥ 0 if ∆ < ( n k(n-1) ) 1/n . Therefore ∆ is ultimately lower bounded by ( n k(n-1) ) 1/n . Finally, since ∆ is lower bounded by a positive number and M is upper bounded, it follows (by direct application of [START_REF] Piazza | An upper bound for the condition number of a matrix in spectral norm[END_REF]) that the condition number This establishes that ẏ is persistently exciting.

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig.1: Evolution of the system/observer pair in 3D space

Fig. 4 :Fig. 5 :

 45 Fig.4: Evolution of the system/observer pair when the output is affected by a noise

Fig. 6 :

 6 Fig. 6: Evolution of the estimate â in the presence of measurement noise

  y(τ )| 2 dτ ≥ δ -µ Let us proceed by contradiction and assume that (5) does not hold, i.e. | ẏ(τ )| < , ∀t ∈ [t, t + δ], then d(δ) < and t+τ t |b T y(τ )| 2 dτ ≥ δ -µ. This contradicts (4) according to which ∀b ∈ S n-1 : t+δ t |π y(τ ) b| 2 dτ = δ-t+δ t |b T y(τ )| 2 dτ ≥ µ Therefore (5) holds true.

4 From 2 -

 42 κ(M ) = ||M ||.||M -1 || ≤ the equation (9) of M one gets d dt (M y ) = y + M ẏ Therefore (M y )(t) = (M y )(0) + t 0 (y + M y )(τ )dτand|(M y )(t)| ≥ t 0 (y + M y )(τ ) -|(M y )(0)| Define: • k M := sup t∈[0,+∞[ ||M (t)||, which implies that |(M y )(t)| ≤ k M , ∀t, • c(t) := sup τ ∈[0,t] | ẏ (τ )|, • w(t) := (y + M y )(t). One has w(τ ) = y (0(M y )(τ ) for some s i ∈ [0, τ ] (i = 1, . . . , n). Therefore | t 0 w(τ )dτ | ≥ | t 0 y (0)dτ |y )(τ )dτ | ≥ f (t, t 1 ) ; t 1 ≥ t with f (t, t 1 ) := t -nc(t 1 ) t 2 2 -k M c(t 1 )t. The function f (., t 1 ) monotically increases for 0 ≤ t ≤ t 2 = 1-k M c(t1) nc(t1) if c(t 1 ) < 1 k M . Define g(c) := t (c) -nc t (c) 2 k M ct (c) with t (c) = 1-k M cnc . The function g decreases monotically on [0, 1k M ] with g(0) = +∞ and g(1/k M ) = 0. Let c > 0 denote the value of c such that g(c) = 2k M , and sett 1 = t 2 = 1-k M c nc . If c(t 1 ) = sup τ ∈[0,t1] | ẏ (τ )| < c then |(M y )(t 1 )| > k M (contradiction). Therefore, sup τ ∈[0,t1] | ẏ (τ )| ≥ c, which proves the existence of a timeinstant t 3 ∈ [0, t 1 ] such that | ẏ (t 3 )| ≥ c.The same proof repeated on every interval [kt 1 , (k + 1)t 1 ] (k ∈ N) shows that | ẏ | is periodically larger than c > 0.

|.| stands for the Euclidean norm of vectors and ||.|| is the induced matrix norm.


Minh-Duc Hua
email: hua@isir.upmc.fr

Daniele Pucci
email: daniele.pucci@iit.it

Tarek Hamel
email: thamel@i3s.unice.fr

Pascal Morin
email: morin@isir.upmc.fr

Claude Samson
email: claude.samson@inria.fr

A novel approach to the automatic control of scale model airplanes

Keywords: aerial vehicles, fixed-wing aircraft, dynamics of flight, nonlinear feedback control

This paper explores a new approach to the control of scale model airplanes as an extension of previous studies addressing the case of vehicles presenting a symmetry of revolution about the thrust axis. The approach is intrinsically nonlinear and, with respect to other contributions on aircraft nonlinear control, no small attack angle assumption is made in order to enlarge the controller's operating domain. Simulation results conducted on a simplified, but not overly simplistic, model of a small airliner illustrate the soundness of the approach.

I. INTRODUCTION

The present paper explores a new approach to the nonlinear control of small underactuated aerial vehicles equipped with fixed planar wings, alike common scale model airplanes. Besides purely scientific issues, this research is motivated by the emerging and rapidly growing market of small automatized aerial drone-type vehicles ranging from VTOL devices (helicopters, quadrotor vehicles,...) to small cruising airplanes used for a multitude of applications (site surveillance, rescue missions, landscape mapping, aerial photography,...). Controlling these devices calls for a unifying theory which is both generic and adaptable to the specificities of each vehicle. As for the usefulness of a nonlinear control approach, it partly resides in the fact that small airplanes are particularly sensitive to wind perturbations and thus confronted to flight conditions involving large and rapidly varying attack angles that existing autopilots, based on linear control techniques, cannot handle properly. The proposed approach departs from traditional approaches based on the linearisation of the vehicle's dynamics about a number of operating conditions and developed primarily for full-size commercial aerial vehicles (airplanes and helicopters) meant to fly along specific trajectories (trim trajectories) that involve a narrow range of attack angles. Many aerial vehicles are required to fly in more diverse conditions that involve large and rapid variations of the attack angle. Examples are given by fighter aircraft, convertible aircraft, or small Unmanned Aerial vehicles (UAVs) operating in windy environments. Vertical Take-Off and Landing (VTOL) vehicles, like e.g. ducted fans and 3D scale model airplanes, are also subjected to large variations of the attack angle when transitioning from hover to horizontal cruising flight. For the automatic control of these devices it matters to ensure a large stability domain that is provably achievable only with a nonlinear control design. Nonlinear feedback control of aircraft can be traced back to the early eighties. Following [START_REF] Singh | Output feedback nonlinear decoupled control synthesis and observer design for maneuvering aircraft[END_REF], control laws based on the dynamic inversion technique have been proposed to extend the flight envelope of military aircraft (see, e.g., [START_REF] Wang | Robust nonlinear flight control of high-performance aircraft[END_REF] and the references therein). The control design strongly relies on tabulated models of aerodynamic forces and moments, like the High-Incidence Research Model (HIRM) of the Group for Aeronautical Research and Technology in Europe (GARTEUR) [START_REF] Muir | Robust flight control design challenge problem formulation and manual: The high incidence research model (hirm)[END_REF]. However, the attack angle is assumed to remain small and away from the stall zone. Should this assumption be violated the system's behavior is unpredictable. Comparatively, nonlinear feedback control of VTOL vehicles is more recent and has been addressed with a larger variety of techniques. Besides dynamic inversion [START_REF] Hauser | Nonlinear control design for slightly non-minimum phase systems: Application to V/STOL[END_REF], other techniques include Lyapunovbased design [START_REF] Marconi | Autonomous vertical landing on an oscillating platform: an internal-model based approach[END_REF], [START_REF] Isidori | Robust autonomous guidance: an internal-model based approach[END_REF], Backstepping [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrotor[END_REF], Sliding modes [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro quadrotor[END_REF], [START_REF] Xu | Sliding mode control of a class of underactuated systems[END_REF], and Predictive control [START_REF] Kim | Nonlinear model predictive tracking control for rotorcraft-based unmanned aerial vehicles[END_REF], [START_REF] Bertrand | Contractive model predictive control of an unmanned aerial vehicle model[END_REF]. A more complete bibliography on this topic is provided in [START_REF] Hua | Introduction to feedback control of underactuated vtol vehicles[END_REF]. However, most of these studies address the stabilization of hover flight or of low-velocity trajectories, and little attention is paid to aerodynamic effects. These are typically either ignored or modelled as a simple additive perturbation the effect of which has to be compensated for by the feedback action. For highly dynamic flight, as in the case of cruising airplanes taking advantage of aerodynamic lift, or in harsh wind conditions, aerodynamic effects cannot be neglected. A proper modelling of aerodynamic effects for the purpose of designing feedback controllers that perform well in a large flight domain, and the design of such controllers whose performance is supported by mathematical analysis are issues that have been seldom be addressed so far by the control and robotics communities. An important difficulty comes from the fact that classical methods used in aerodynamic modelling to precisely describe aerodynamic forces, e.g. computational fluid dynamics (CFD) or wind tunnel measurements, do not provide analytical expressions of aerodynamic characteristics. From a control design perspective these methods are useful to finely tune a controller around a given flight velocity. But they are not best suited to handle the case of large flight envelopes that involve strong variations of either the flight velocity or the attack angle. We instead advocate the use of simple analytical models of aerodynamic characteristics. Although relatively imprecise, these models account for important structural properties of the system in a large flight envelope. The idea behind our approach is to exploit these properties at the control design level and rely on the robustness of feedback controllers to cope with discrepancies between the model and the true aerodynamic characteristics. The control design here proposed follows the one reported in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] and extends it to the case of fixed-wing aircraft whose shape satisfies a symmetry property weaker than the axisymmetry property considered in this previous work. In particular, we show how the stabilization of the airplane's sideslip angle via roll control -a control which is not needed in the axisymmetric case to track reference (velocity or position) trajectories, but with which every airplane's human pilot is familiar-brings the control problem back to the one of the axisymmetric case for the determination of the other control variables. The paper is organized as follows. Section II provides notation and background material. The main results are reported in Section III. First, structural properties of the aerodynamic force induced by the symmetry of the airplane are pointed out and an expression of this force, depending essentially on the attack angle when the sideslip angle is equal to zero, is worked out. We then show how the sideslip angle can be stabilized at zero via roll feedback control. Zeroing this angle achieves what is commonly called the coordinated flight between roll and yaw control actions. Once the sideslip angle vanishes, the aerodynamic force applied to the aircraft is similar to the one applied to an axisymmetric vehicle, and the other control variables can be determined according to the approach presented in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. The control design is complemented with a convergence and stability result for a class of reference trajectories. In Section IV simulations performed on a simplified, but nonetheless physically credible, model of a small AUV illustrate the robust performance of the proposed controller. Concluding remarks and perspectives are proposed in Section V.

II. BACKGROUND AND NOTATION

A. Notation

Throughout the paper, E 3 denotes the 3D Euclidean vector space and vectors in E 3 are denoted with bold letters. Inner and cross products in E 3 are denoted by • and × respectively. Let us consider a rigid body immersed in air, as represented by Fig. 1. The following notation is used.

• G denotes the body's center of mass (CoM); • I = {O; i 0 , j 0 , k 0 } is an inertial frame; • B = {G; i, j, k} is a body-fixed frame;

• ω is the angular velocity of the body-fixed frame with respect to (w.r.t.) the inertial frame; • m is the body's mass;

• v is the CoM's velocity w.r.t. the inertial frame;

• a is the CoM's acceleration w.r.t. the inertial frame;

• g = g k 0 is the gravitational acceleration;

• v a = v-v w is the body's air-velocity (i.e., translational velocity w.r.t. the ambient air), with v w denoting the wind's velocity. The direction of v a in the body-fixed frame is often represented by the two angles α and β denoting respectively the attack angle and sideslip angle, i.e.

v a = |v a |(sin α i + cos α(sin β j -cos β k)) (1)
This definition of attack and sideslip angles slightly differs from the one classically used in aeronautics, but it is locally equivalent when these angles are small. Note also that, alike any parametrization of the unit sphere by two angles, this representation is not uniquely defined everywhere. More precisely, β cannot be defined by continuity at α = ±π/2, nor when |v a |=0. The coordinate vector of any ξ ∈ E 3 w.r.t. the body-fixed frame B is denoted by the ordinary letter ξ, i.e. ξ = ξ 1 i + ξ 2 j + ξ 3 k. Throughout the paper, either ξ or d dt ξ denotes the time-derivative of ξ ∈ E 3 relatively to the inertial frame I.

B. Aerodynamic forces

The resultant aerodynamic force F a applied to a rigid body moving with air-velocity v a is traditionally decomposed into the sum of a drag force F D along the direction of v a and a lift force F L perpendicular to this direction, i.e.

F a = F D + F L (2)
The

F D = -η a |v a |C D v a , F L = η a |v a |C L v ⊥ a ( 3
)
with

• v ⊥ a a vector perpendicular to v a and with same norm, • η a := ρΣ 2 with ρ the free stream air density, and Σ an area germane to the body shape. Another expression of drag and lift forces can be obtained by introducing an air-frame A = {G; i a , j a , k a } defined as follows (See Fig. 1):

k a = - v a |v a | , i a = - v ⊥ a |v a | , j a = k a × i a (4)
Then (3) can also be written as

F D = η a |v a | 2 C D k a , F L = -η a |v a | 2 C L i a (5)
Throughout the paper we assume that the aerodynamic characteristics and the lift-force direction (i.e., the vector i a ) depend on the direction of v a but not on its intensity. We also neglect the influence of Reynolds and Mach numbers variations on these functions. In other words,

C D = C D (k a ), C L = C L (k a )
, and i a = i a (k a ). As commonly done in the dynamics of flight literature, we may also write, with a slight abuse of notation,

C D = C D (α, β), C L = C L (α, β)
when the attack and sideslip angles are well defined.

C. Control problem statement

The control objective here chosen is the asymptotic stabilization of either a desired (or reference) velocity v r or a desired position p r for the vehicle's CoM. It is assumed that the control inputs consist of a thrust force T = -T k, with T denoting the thrust intensity, and a torque vector Γ, independent of the thrust actuation, used to modify the body's angular velocity ω -and thus the body's orientation. Nominally the resultant thrust force passes near the center of mass and thus exerts a small parisitic torque. Aerodynamic forces that apply to the body may also create a resultant parasitic torque Γ a . The control torque Γ then has to precompensate for these parasitic torques and also ensure that (almost) any desired angular velocity is obtained rapidly (instantaneously, in the ideal case). This leads to consider the angular velocity ω as an intermediate control input, with the conceptually important advantage of not having to take the specificities of the physical torque actuation at this stage of the control design. This corresponds to a standard "backstepping" assumption. Once it is made, the vehicle's actuation consists of four input variables, namely, the thrust intensity and the three coordinates of ω.

The control model then reduces to the classical Newton's equation complemented with the orientation kinematics equation:

m v = mg + F a -T k (6a) d dt (i, j, k) = ω × (i, j, k) (6b)

III. MAIN RESULTS

A. Modelling of aerodynamic forces for control design

The following way of writing the aerodynamic force is instrumental for the proposed control approach.

Proposition 1 Assume that k a,1 = 0. Then, the aerodynamic force can be expressed as

F a = -η a |v a | C D -C L i a,1 k a,1 v a - C L |v a | k a,1 i × j a (7)
Proof: Recall the classical vectorial relation:

∀x, y, z ∈ E 3 , x × (y × z) = (x • z)y -(x • y)z (8) Since j a = k a × i a , it follows from (8) that i × j a = (i • i a )k a -(i • k a )i a = i a,1 k a -k a,1 i a Therefore, i a = 1 k a,1 (i a,1 k a -i × j a )
and ( 7) readily follows from ( 2) and ( 5).

Proposition 1 holds independently of the vehicle's shape. In order to further specify aerodynamic forces, we introduce an assumption that is satisfied for almost all aerial vehicles.

Assumption 1

The vehicle is symmetric with respect to the plane {G, i, k}, i.e.

(xi + yj + zk) ∈ S =⇒ (xi -yj + zk) ∈ S
with S the vehicle's surface.

For control design purposes, it is useful to have a more precise expression of the lateral component F a,2 of the aerodynamic force, because this component determines the variations of the lateral velocity component v a,2 , and thus the variations of the sideslip angle β. The following proposition provides such an expression under mild assumptions.

Proposition 2 Assume that 1) C L is a smooth function of (α, β) around α = 0 and C L (0, β) = 0 ∀β; 2)
The vector i a varies smoothly with k a ; Then, the following relationship is satisfied:

F a,2 = -a(k a )|v a |v a,2 + b(k a )v a,1 v a,2 (9)
with a(•) and b(•) denoting continuous functions, and a(•) strictly positive.

The proof of this proposition is given in the appendix.

Remark:

The first condition in Proposition 2 is satisfied when the plane {G, j, k}, which is orthogonal to the plane of symmetry evoked in Assumption 1 and to which the thrust direction is parallel, corresponds to the zero lift-plane, i.e., α = 0 =⇒ F L = 0. As for the second (regularity) condition, one can reasonably expect that it is satisfied in practice, even though a rigorous mathematical justification remains to be worked out.

Another important consequence of the vehicle's symmetry is that if v a belongs to the plane {G, i, k} , i.e. if the sideslip angle is zero, then the lift force F L also belongs to this plane. Also j a = j, k a = -sin α i + cos α k, and one can set i a = cos α i + sin α k since i a is orthogonal to k a and belongs to the plane {G, i, k} in this case.

Let F 0 a denote the aerodynamic force when β = 0. In view of [START_REF] Kim | Nonlinear model predictive tracking control for rotorcraft-based unmanned aerial vehicles[END_REF] and the above relations one deduces that

F 0 a = -η a |v a | C 0 D + C 0 L cot α v a + C 0 L |v a | sin α k ( 10
)
where the 0 superscript denotes evaluation at β = 0, i.e.,

C 0 D (α) = C D (α, 0), C 0 L (α) = C L (α, 0).
It then matters, for control purposes, to know more about the lift and drag coefficients and, in particular, on the way they depend on the attack angle. For small angles, a common model consists of choosing C 0 D constant and C 0 L a linear function of the attack angle, i.e. (C 0 D = c 0 , C 0 L = 2c 1 α) for some positive numbers c 0 and c 1 . With these coefficients one readily verifies that the aerodynamic force F 0 a of relation [START_REF] Pucci | Flight dynamics and control in relation to stall[END_REF] is, in the first approximation, equal to

F 0 a = -η a |v a | c0 v a -2c 1 v a,3 k (11)
with c0 = c 0 +2c 1 , when α is small. The above expression of the aerodynamic force is remarkable because it indicates that this force is the sum of two terms: one along the air velocity v a and independent of the vehicle's orientation, and the other one along the thrust direction. For control purposes, these are good terms. Indeed, since the first term is independent of the vehicle's orientation, the problem of aligning the vehicle's thrust direction with it via the angular velocity control ω is a well posed problem (see, for instance, [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF] or [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] for an (almost) globally exponentially stabilizing solution). As for the second term that is aligned with the thrust direction, it can be compensated for directly via a suitable choice of the thrust intensity T . This is an important remark because it explains how orientation control and thrust control can be decoupled and, to some extent, why human pilots are capable of manually controlling airplanes without the help of control and stability augmentation systems (CSAS). Now, in order to address situations when the attack angle is not small one needs expressions of the drag and lift coefficients over the whole range of possible angles. These expressions have to be compatible with the choice pointed out before for small angles, and they have to be reasonably representative of the physical reality. It is also useful to observe that, in order to preserve the property of independence w.r.t. the vehicle's orientation of the first term of F 0 a , a sufficient condition is the invariance of (C 0 D (α) + C 0 L (α) cot α) w.r.t. the attack angle. Examples of functions that satisfy this condition in a large domain of attack angles and are compatible with measurements performed on NACA wing profiles and axisymmetric shapes are reported in [START_REF] Pucci | Flight dynamics and control in relation to stall[END_REF] and [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. The most important limitation comes from the stall phenomenon corresponding to the brutal reduction of lift occurring beyond a certain value of the attack angle and in a (usually) narrow range of this angle (the stall zone). In the stall zone, the aforementioned invariance condition cannot be satisfied.

For the sake of simplifying the control design exposition, we will consider from now on the particular set of functions specified by the following assumption: Assumption 2 The vehicle's aerodynamic characteristics C D , C L are such that:

C 0 D (α) = c 0 + 2c 1 sin 2 α , C 0 L (α) = c 1 sin 2α (12)
with c 0 and c 1 two positive numbers.

These functions are π-periodic and respectively even and odd w.r.t. α, they generalize the linear approximations evoked previously, and they satisfy the invariance condition since one readily verifies that C 0 D (α)+C 0 L (α) cot α = c 0 +2c 1 (= c0 ). These properties are consistent with the aerodynamic characteristics of a body shape symmetric with respect to both the planes {G, j, k} and {G, i, j}, as in the case of NACA00XX wing profiles. For more general body shapes, C 0 D and C 0 L as defined by ( 12) can still be used as an approximation of many body aerodynamic characteristics in a large flight domain [START_REF] Pucci | Flight dynamics and control in relation to stall[END_REF] [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. The aerodynamic force expression, when the sideslip angle is kept equal to zero, is then given by [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF]. On the basis of this expression, and of the expression (9) of the lateral component of the aerodynamic force, the control design that we propose consists in i) controlling the sideslip angle to zero via the control input ω 3 , and ii) monitoring the control thrust force via T and ω 1,2 so as to compensate for F 0 a -complemented with the force of gravity-and tracking errors.

B. Sideslip angle control

We assume from now on that the wind's velocity is constant w.r.t. the inertial frame, i.e. vw ≡ 0. Since v a is the coordinate vector of v a in body-frame, it follows from (6) that

m va = mv a × ω + mg + F a -T e 3
The control problem of the sideslip angle to zero is equivalent to controlling v a,2 to zero. Then,

va,2 = v a,3 ω 1 -v a,1 ω 3 + g 2 + 1 m F a,2
Using the expression (9) of F a,2 given in Proposition 2, this equation can be written as

va2 = v a,3 ω 1 -v a,1 ω 3 +g 2 - 1 m (a|v a |v a,2 -bv a,1 v a,2 ) (13)
If we calculate the third component of the angular velocity as follows

ω 3 = v a,3 ω 1 + g 2 v a,1 + μsign(v a,1 )v a,2 (14)
with μ ≥ 0 (not necessarily constant), then the resulting closed-loop equation is

va,2 = - a m |v a | - b m v a,1 + μ|v a,1 | v a,2
which implies the asymptotic stability of v a,2 = 0 with an exponential rate of convergence to zero, provided that μ ≥ b m and |v a,1 | is strictly larger than a positive number. The control expression ( 14) calls for a few remarks and adaptations. This control is defined provided that v a,1 , and thus the attack angle α, is different from zero. While this condition is generically satisfied during a "standard" cruising flight, the attack angle may pass through zero during acrobatic manoeuvres or in perturbed airflows. Assuming that the vehicle can hover at zero velocity without wind, then the air velocity is also zero and division by v a,1 is again ill-defined. The corresponding physical problem is that roll modification has no instantaneous influence of the sideslip angle when the attack angle vanishes. It is then better to impose a zero roll control. Note also that the problematic term in ( 14) is the pre-compensation term (v a,3 ω 1 +g 2 )/v a,1 .

In practice it may be sufficient to keep only the feedback part of the control expression, with a large enough gain μ. This is basically what a human pilot does when using a balance indicator (inclinometer), or a "turn-and-slip indicator", to achieve "balanced flight". Another possibility consists in using a term that imperfectly pre-compensates for the term (v a,3 ω 1 +g 2 ) in ( 13), but is always well defined. For example, instead of ( 14), one may consider the following continuous control expression

ω 3 = v a,1 v 2 a,1 + (v a,3 ω 1 + g 2 ) + μ v a,1 |v a,1 | + v a,2 (15)
with > 0 a small constant.

C. Thrust control

Let us assume for the time being, and for the sake of simplifying the presentation of the control design, that roll control is active and fulfills its role of keeping the sideslip angle β equal to zero. Then F a = F 0 a , as specified by [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF]. Consider a reference velocity v r , which may be the time derivative of a reference position p r , and denote the reference acceleration w.r.t. the inertial frame as a r := d dt v r . Velocity and position errors are respectively denoted as ṽ := vv r and p := pp r . From Newton's equation ( 6) one deduces

m v = -T k + m(g -a r ) + F 0 a (16)
Using [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF], this equation can in turn be written as

m v = -T p k + F p ( 17
)
with [START_REF] Hauser | Nonlinear control design for slightly non-minimum phase systems: Application to V/STOL[END_REF] and F p := -η a c0 |v a |v a + m(ga r ). The main interest of this writing is that F p does not depend on the vehicle's orientation. This equation is similar to the one derived in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] for axisymmetric vehicles.

T p := T -2c 1 η a |v a |v a,
It is also formally the same as the equation of a spherical body subjected to an aerodynamic force reduced to the drag component -η a c0 |v a |v a . This is what we have called spherical equivalence. Once it is established, the control design essentially reduces to transposing the control laws derived for the simpler spherical case. Let ( ρ, ṽ) denote a "generalized" error vector. The precise definition of ρ depends on the control objective. For instance, it can be equal to p, if the objective is reference position tracking, or it can be an integral of the velocity error ṽ, if the control objective is reference velocity tracking and an integral correction term is needed to improve the control performance. In the case of reference position tracking it may also be useful to introduce an integral correction. Then ρ = (I p , p), with I p denoting a saturated integral of the position tracking error. The simplest case corresponds to pure velocity control without integral correction, for which ρ = ∅.

The idea is to transform the initial control problem into the much simpler one associated with the linear control system v = ξ, with ξ denoting an intermediary feedback control input in charge of ensuring the exponential stability of the generalized error vector ( ρ, ṽ) to zero for this system. In view of (17), the relation v = ξ can be satisfied if and only if

T p k = Fp with Fp := F p -mξ (18)
This in turn implies that T p = Fp • k, so that

T = (F 0 a + m(g -a r -ξ)) • k = -η a c 0 |v a |v a3 + m(g -a r -ξ) • k , ( 19
)
Also k should be equal to k r := ± Fp /| Fp |. The preference given to a positive thrust yields the choice of positive sign for this latter equality, i.e.

k r = Fp | Fp | (20)
From there, the stabilization of ( ρ, ṽ) to zero essentially relies on the asymptotic stabilization of the body's frame vector k to the unit "reference" vector k r via the choice of ω 1,2 . The following theorem, which extends a result in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] for axisymmetric vehicles, recalls a solution to this problem, and states a convergence and stability property for a class of reference trajectories.

Theorem 1 Assume that 1)

The attack angle along the reference trajectory, denoted as α r , is always strictly positive, i.e., ∃ > 0 : ≤ α r , ∀t; 2) Fp does not vanish along the velocity reference trajectory, i.e., ∃δ > 0 : δ ≤ Fp (v r,a (t), a r (t), 0), ∀t, with v r,a := v rv w . 3) ξ(ṽ) makes ( ρ, ṽ) = 0 a locally exponentially stable equilibrium point of the system:

ρ = f ( ρ, ṽ) (21a) v = ξ (21b)
Let T be computed according to (19), and the angular velocity ω be computed according to Then, (T, ω) ensures local exponential stability of the equilibrium point ( ρ, ṽ, k) = (0, 0, k r ) for the system (21a)-( 16)-(6b).

ω 1 = ω r,1 -(ν + γ γ )k r,2 ω 2 = ω r,2 + (ν + γ γ )k r,1
The proof of this theorem is a straightforward extension of the proof of Proposition 4 in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. The main difference between the two results is the smaller set of reference trajectories for which local asymptotic stabilization is proved, due to the first assumption in the theorem that is needed to ensure local uniform exponential convergence of the sideslip angle to zero. Remarks: R.1) Although only local exponential stability is proved, the practical stability domain can be quite large. Indeed, the difficulty to establish a more global stability result comes from the necessity of avoiding the "thin set" of singularities associated with vehicle's orientations such that k • k r = -1 and the passage through zero of | Fp | and α. A complement to the analysis could consist in modifying the control inputs near these singularities so as to ensure unconditional wellposedness of the control expressions and prove a "practical" stability result that is more global, with ultimate boundedness of the tracking errors instead of their convergence to zero. R.2) The assumptions imposed on the reference trajectory are satisfied by a large class of reference trajectories, including standard cruising flight. However, they can be violated when the reference trajectory involves specific aggressive manoeuvres. Another complement to this study could thus consists in working out more explicit characterizations of these "good" reference trajectories.

IV. SIMULATION RESULTS Fig. 2. Simulated fixed-wing UAV composed of an elliptical fuselage and several disk-shaped wings and tails.

We have simulated a fixed-wing UAV sketched in Fig. 2. Size, mass, and inertia specifications of the aircraft are given in Tab. I and Fig. 2. The vehicle's CoM is located at the ellipsoidal fuselage center. Wings and tails (horizontal and vertical) are decomposed into a set of disk-shaped flat surfaces assembled together as shown in Fig. 2. A reason for this decomposition is that a disk, alike an ellipsoid, is an axisymmetric body for which explicit models of the aerodynamic force exerted on it, in the form of functions depending on the disk's center velocity vector (relatively to the ambient air) and the disk's orientation w.r.t. this velocity vector, are available [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. Using a superposition principle, which does not hold stricto sensu but is commonly used in simulation packages, the resultant aerodynamic force is computed as the sum of the forces applied to all parts (disks+fuselage) taken separately. The resultant parasitic torque is computed in the same manner. Another interesting feature associated with this procedure is that the resultant aerodynamic force is no longer a function of only the aircraft's CoM velocity v a and orientation, contrary to the simplified model used for the control design. It also depends on the aircraft's angular velocity ω, since each disk center velocity is a function of v a and ω. The resulting model is for this reason more complete and precise than the one used for the control design, with the complementary interest that the mismatch between the two models is useful to test the control robustness against modelling errors. In order to simulate the effect of the stall phenomenon, the functions used for the lift and drag coefficients involved in the calculation of the aerodynamic forces applied to the disks composing the wings and tails are chosen so as to closely approximate the measured coefficients of the NACA0018 airfoil at R e = 160000 (see [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF] for complementary details). These functions, and the "ideal" ones defined by [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] and used for the control calculations, are drawn in Fig 3. Discrepancies between these two sets of functions represent another opportunity to test the controller's robustness.

For these simulations, the control variables are composed of the thrust intensity T complemented with a control torque vector Γ, instead of the angular velocity ω. The angular velocity defined in Theorem 1 is used as a desired angular velocity ω , and the role of the torque Γ is to make the aircraft's angular velocity ω closely tracks ω . This is achieved by computing the control torque vector as follows Γ = -K ω J(ωω ) + ω × Jω with K ω denoting a high-gain positive matrix here chosen equal to diag (30,30,30). In this calculation, precompensation of the aerodynamic torque and of the reference angular acceleration is deliberately omitted to further test the robustness of the control design against modelling errors. As for the thrust intensity T , it is computed according to (19) when the result is positive, and otherwise kept equal to zero in order to be consistent with the positive thrust constraint of standard airplanes and also test the recovery possibilities offered by the control once a positive thrust is again required to track the reference trajectory. The control objective is the asymptotic stabilization at zero of the position error p with the reference position p r obtained by integration of the reference velocity defined by

v r (t)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (0, 0, 0) 0 ≤ t < 5 (t -5, 0, 0) 5 ≤ t < 15 (10, 0, 0) 15 ≤ t < 20 (10 cos( 2π(t-20) 5
), 0, -10 sin( 2π(t-20)

5

)) 20 ≤ t < 25 (10, 0, 0) 25 ≤ t < 30 (10 cos( 2π(t-30) 10 ), 10 sin( 2π(t-30) 10

), 0) 30 ≤ t < 40 (10, 0, 0) 40 ≤ t < 45 (55 -t, 0, 0) 45 ≤ t < 55 (0, 0, 0) t ≥ 55
The vehicle is initially hovering at the origin of the inertial frame. This velocity profile has been chosen for the vehicle to perform several complex manoeuvres. The intermediary control vector ξ is a PID controller determined as follows

ξ( ρ, ṽ) = -k p ( p + k i I p ) -k v (ṽ + k i d dt I p )
with ρ = (I p , p), I p a numerically calculated saturated integral of p such that d dt I p = k pz (-I p + sat Δ (I p + p kpz )). The gain parameters used for the calculations are k p = 3.24, k v = 2.5456, k i = 0.648, k pz = 1, Δ = 4/k i . To further assess the robustness of the controller, several parametric errors have been simulated. For instance, we have used m = 2.85 and J = diag(0.04, 0.1, 0.1) in the control calculations, instead of the values reported in Table I.

Simulation results are reported next. The projection on the horizontal and vertical planes of the path followed by the vehicle's CoM is shown in Fig. 4. Variations w.r.t. time of the vehicle's position and orientation, of the attack and sideslip angles, and of the thrust intensity and torque inputs are shown in Figs. 5678. One observes from Figs. 4 and5 that the vehicle closely follows the reference trajectory.

The transient manoeuvre, from hovering to high speed horizontal cruising flight, between t = 5(s) and t = 20(s), is smoothly negotiated (see Figs. 67). During this transient period the attack angle varies in large proportions, from 90(deg) to 5(deg). At horizontal cruising speed of 10(m/s) the attack angle converges to 5.8 (deg) approximately. The passage through the stall zone (i.e., α ∈ [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF][START_REF] Wang | Robust nonlinear flight control of high-performance aircraft[END_REF](deg)), when t ≈ 14(s), does not jeopardize the system's stability. The transient manoeuvre from high speed cruising to hovering, between t = 40(s) and t = 55(s), is also performed smoothly.

Concerning the aerobatic stunt induced by the reference velocity variations between t = 20(s) and t = 25(s), the position error remains small (< 1.2 (m)) during the first part of the loop, then grows to ≈ 5(m) at the end of the loop essentially because of the positivity imposed to the thrust intensity during the time interval t ∈ [22.7, 24.6](s) (see Fig. 8-top). Despite this, the vehicle manages to complete the stunt (see Fig. 4-bottom), then catches up and closely track the reference trajectory when the thrust intensity calculated with the control law becomes again positive.

Finally, the turn manoeuvre between t = 30 (s) and t = 40 (s) is successfully executed, and a small tracking error (| p| < 1.5 (m)) is maintained along a circular horizontal motion. One observes from Fig. 6-bottom and Fig. 7-bottom that the roll angle and the sideslip angle converge to constant values (≈ 33.8 (deg) and ≈ 3.5(deg), respectively). The fact that the sideslip angle does not converge to zero can be explained by the non-compensation of aerodynamic torques by the control torque. Fig. 7-bottom also shows that, at the time instants t = 30 (s) and t = 40 (s), the sideslip angle has short bursts essentially because of the discontinuity of the reference acceleration at these time instants. Finally, one can also observe from Fig. 7-top that the attack angle varies in large proportions during the hovering manoeuvre after t = 54 (s). This is coherent with the fact that the determination of this angle is ill-conditioned when the body's air-velocity v a vanishes.

V. CONCLUDING REMARKS

A novel approach to the automatic control of scale model airplanes, which extends earlier studies by the authors on the control of VTOL and axisymmetric vehicles, has been presented. With respect to other studies on the subject, and in order to extend the control operating domain, a salient feature of the approach is the use of simple analytical models of the aerodynamic forces applied to the aircraft that cover the whole range of attack angles. Implementation aspects (actuation and sensors related issues, in particular) and experiments conducted on physical scale model airplanes will be the subject of forthcoming studies.

ACKNOWLEDGEMENT

This work was supported by the French Agence Nationale de la Recherche through the ANR ASTRID SCAR project "Sensory Control of Aerial Robots" and by the Chaire d'excellence en robotique RTE-UPMC "Mini-drones autonomes".

APPENDIX

Proof of Proposition 2:

From (2), (3), and ( 5 so that ( 9) is satisfied with

a(k a ) = η a C D (k a ) , b(k a ) = -η a C L (k a ) k a,1 i a,2 k a,2
The function a is continuous and strictly positive since C D also satisfies these properties and η a > 0 is constant. Since C L is smooth and vanishes identically at α = 0, π, which corresponds to k a,1 = 0, the function CL(ka) ka,1 can be defined by continuity at these points. There remains to show that ia,2 ka,2 can be defined by continuity at k a,2 = 0. This follows from assumption of smoothness of i a w.r.t. k a and to the fact that k a,2 = 0, implies i a,2 = 0.

Fig. 1 .

 1 Fig. 1. Sketch representation of a scale model airplane

ω 3 :

 3 given by (14)(22) with k r defined by (20), ω r := k r × kr , γ := c 2 + | Fp | 2 with c 2 any positive constant, and ν(k, t) any continuous positive real-valued function such that inf k,t k 1 (k, t) > 0.

Fig. 3 .

 3 Fig. 3. Approximated analytical model (with stall) and "ideal" model of lift and drag coefficients v.s. attack angle at Re = 160000 for NACA0018 airfoil.

Fig. 4 .

 4 Fig. 4. Reference and vehicle trajectories projected on the horizontal and vertical planes.

Fig. 5 .

 5 Fig. 5. Reference and vehicle positions versus time.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Vehicle's orientation (elevation, heading, roll angles) versus time.

  Moment of inertia J [kg.m 2 ] diag(0.033, 0.13, 0.13) Wing total area Σ w [m 2 ] 0.7267 Horizontal tail total area Σ t [m 2 ] 0.1135 Vertical tail total area Σ vt [m 2 ] 0.045 Fuselage total area Σ f [m 2 ] 0.25

	TABLE I
	SPECIFICATIONS OF THE SIMULATED FIXED-WING UAV
	Specification	Numerical Value
	Mass m [kg]	3

  ), F a,2 = -η a |v a |C D v a,2η a |v a | 2 C L i a,2 = -η a |v a |C D v a,2η a |v a | 2 k a,1 CL ka,1 i a,2 = -η a |v a |C D v a,2 + η a |v a |v a,1 CL ka,1 i a,2where the last equality follows from (1). This equality can in turn be written asF a,2 = -η a |v a |C D v a,2 + η a |v a |k a,2 v a,1 = -η a |v a |C D v a,2η a v a,2 v a,1

		CL	ia,2
		ka,1	ka,2
	CL	ia,2
	ka,1	ka,2