Keywords: Fixed-wing aircraft, Nonlinear control design, Aerial robotics, Trajectory-Tracking, Path-Following, Hardware-in-the-loop simulations, Flight experiments. Résumé Drones à voilure fixe, Contrôle non-linéaire, Robotique aérienne, Suivi de trajectoire, Suivi de chemin, Simulation AOA Angle Of Attack. BLDC BrushLess Direct Current. CAS Control Augmentation System. CFD Computational Fluid Dynamics. CPU Central Processing Unit. CRGIS Cultural Resources Geographic Information. EKF Extended Kalman Filter IC Integrated Circuit. xxi IMU Inertial Measurement Unit. INS Inertial Navigation System. LiPo Lithium Polymer. LPV Linear Parameter Varying. LQG Linear Quadratic Gaussian. LQR Linear Quadratic Regulator. MEMS Micro-Electro-Mechanical Systems. MIMO Multi Input Multi Output PI Proportional Integral. PID Proportional Integral Derivative. PWM Pulse Width Modulation. SAS Stability Augmentation System. SISO Single Input Single Output. SLAM Simultaneous Localization And Mapping. SRAM Static Random Access Memory. STOL Short TakeOff and Landing. UAV Unmanned Aerial Vehicle. UDP User Datagram Protocol. VTOL Vertical TakeOff

The present thesis develops a new control approach for scale-model airplanes. The proposed control solutions exploit a simple but pertinent nonlinear model of aerodynamic forces acting on the aircraft. Nonlinear hierarchical controllers are derived on the basis of theoretical stability and convergence analyses. First, the trajectory tracking problem is addressed by extending the thrust vectoring method used for small rotor vehicles to the case of fixed-wing aircraft by compensating for the orientation-dependent aerodynamic forces and by achieving a balanced flight. Then, the path-following control problem is addressed and kinematical guidance and dynamical control laws are developed within a single coherent framework that applies to almost all regular 3D paths. The proposed control laws incorporate integral terms that robustify the control. They are also complemented by addressing several practical issues, and validated via hardware-in-the-loop simulations. Finally, successful flight test results illustrate the soundness and performance of the proposed control design.

duit à l'ère moderne de la robotique aériennes et au développement de drones non militaires à échelle réduite. Ceci est illustré par exemple par le développement de microprocesseurs plus performants, la disponibilité de capteurs inertiels à base de technologies de MEMS et le développement des technologies de batteries. Ces petits appareils se distinguent des avions conventionnels par leur enveloppe de vol élargie et leur interaction avec un environnement dynamique et complexe. Cela a stimulé la recherche dans ce domaine. Les principales plates-formes existantes sont les véhicules à voilure tournante équipés de pales rotatives et capable de vols stationnaires, les véhicules à voilure fixe équipés d'ailes profilées et les véhicules hybrides qui combinent les capacités des deux catégories précédentes.

• Chapitre 2: Bases de la mécanique du vol Les équations de mouvement d'un avion peuvent être dérivées en utilisant la mécanique classique des corps rigides. Dans ce chapitre, nous définissons les variables d'état et les entrées de contrôle utilisées pour modéliser le système, puis nous présentons les équations cinématiques et dynamiques qui décrivent l'évolution de ces variables d'état. Afin de compléter la description dynamique du système, nous présentons également un modèle de la poussée lorsqu'elle est produite par un moteur à hélice.

• Chapitre 3: Aérodynamique des avions Ce chapitre reprend certains concepts de base de l'aérodynamique et introduit les principes qui ont été exploités lors de la conception des véhicules volants plus lourds que l'air. L'objectif est de décrire et de modéliser les forces et des moments aérodynamiques agissant sur les profils aérodynamiques, les ailes finies et enfin sur un avion dans son ensemble. Les hypothèses qui peuvent être faites pour les régimes de vol à faible vitesse sont également expliquées. En conclusion, il est souligné que dans la littérature classique, les coefficients correspondants aux forces aérodynamiques sont modélisés comme étant linéaires en angle d'attaque (attack angle) et en angle de dérapage (sideslip angle). Ce sont des conventions différentes de celles adoptées dans ce travail.

• Chapitre 4: Objectifs de contrôle Les objectifs de contrôle sont définis, et la distinction entre le problème de stabilisation de trajectoire de référence et celui du suivi de chemin est clarifiée. Il est souligné que le suivi de chemin est plus communément envisagé pour les aéronefs à voilure fixe du fait qu'il n'est pas essentiel d'imposer des contraintes de temps strictes dans le cas d'un vol en croisière. Une partie importante de la littérature existante sur le contrôle des aéronefs aborde le sujet du suivi de chemin principalement en le divisant en deux sous-problèmes: le guidage cinématique et le contrôle dynamique. Le guidage cinématique est conçu au niveau cinématique, tandis que le contrôle dynamique est plus compliqué et prend en compte les forces exercées sur l'aéronef, notamment les forces aérodynamiques. Classiquement, les techniques de contrôle linéaire sont appliquées en linéarisant le système autour de trajectoires d'équilibre ("trim trajectories"). D'autres approches non linéaires ont également été développées, certaines basées sur l'inversion dynamique non linéaire combinées avec une architecture hiérarchique de contrôle.

• Chapitre 5: Modèle de commande Le contrôle proposé dans ce travail est basé sur une architecture hiérarchique qui offre de nombreux avantages théoriques et pratiques. Le problème de commande est d'abord simplifié en considérant que la vitesse angulaire joue le rôle d'un terme de contrôle intermédiare, ce qui permet de développer des lois de contrôle indépendantes de la configuration des actionneurs du véhicules. Ensuite, un modèle générique non linéaire de forces aérodynamiques est proposé dans le but de l'intégrer dans la conception de contrôle. Malgré sa simplicité, ce modèle est physiquement pertinent et sera utilisé par la suite pour la synthèse de lois de commande.

• Chapitre 6: Contrôle d'attitude La stabilisation d'attitude correspond à la conception de la boucle interne rapide de l'architecture de contrôle hiérarchique.

Dans ce chapitre, nous dérivons d'abord une loi de contrôle pour la vitesse angulaire comme variable de contrôle intermédiaire avec une analyse de convergence.

Ensuite, nous montrons comment atteindre cette vitesse angulaire désirée avec les surfaces de contrôle d'un avion conventionnel. Enfin, la commande précédente est adaptatée aux avions à deux axes qui manquent de contrôle de gouvernail comme par exemple les ailes volantes.

• Chapitre 7: Stabilisation de trajectoire de référence Une solution au problème de stabilisation de trajectoire de référence est proposée dans ce chapitre en étendant la méthodologie du "thrust vectoring" qui est souvent choisie pour le contrôle des véhicules à voilures tournantes tel que le quadrotor.

Dans le cas des véhicules à voilure fixe, le principal défi est de prendre en compte les forces aérodynamiques d'un avion ce qui complique la conception des systèmes de commande. Afin de résoudre ce problème, le contrôle porposé est conçu et analysé sur la base du modèle de forces aérodynamiques décrit précédemment dans le chapitre 5.

Dans un premier temps, une analyse de l'existence d'orientations d'équilibre pour le véhicle le long d'une trajectoire de référence est présentée. Ceci mène à la définition d'un ensemble de trajectoires de référence admissibles pour lesquelles on peut assurer l'existence d'une orientation d'équilibre. Ensuite, nous présentons la conception du contrôleur de position qui stabilise asymptotiquement la position du véhicule. Enfin nous montrons des résultats de simulation impliquant des trajectoires de référence agressives.

• Chapitre 8: Suivi de chemin Le but de ce chapitre est d'adapter la solution de stabilisation de trajectoire de référence du chapitre précédent au problème du suivi de chemin. L'objectif est de s'approcher d'un chemin géométrique et de le suivre à vitesse donnée, cela sans introduire de contraintes temporelles strictes sur la position de l'avion. L'erreur de position est définie comme étant la projection orthogonale du centre de masse du véhicule sur le chemin. La norme de la vitesse est régulée par la poussée, alors que les surfaces de contrôle permettent de "guider" le véhicule pour assurer la convergence de l'erreur de position à zéro. Des solutions pour les boucles de guidage cinématique et du contrôle dynamique sont présentées. Les lois de contrôle sont théoriquement justifiées par une analyse de stabilité et de convergence.

Des problèmes pratiques complémentaires sont traités, notamment la possibilité de stabiliser la vitesse air plutôt que la vitesse inertielle. Un estimateur du vecteur vitesse air basé sur un modèle dynamique est également proposé.

Les résultats de simulation Hardware-in-the-Loop, impliquant des trajectoires de référence complexes, illustrent la performance et la robustesse de la stratégie de contrôle.

• Chapitre 9: Essais en vol La stratégie de suivi de chemin est mise en pratique dans le cadre d'une série d'essais en vol impliquant un drone à échelle réduite. Nous présentons d'abord les propriétés de l'avion et l'architecture matérielle. Ensuite, nous décrivons le déroulement de l'expérience. Enfin, nous montrons des résultats qui confirment la robustesse de la commande par rapport aux erreurs de modélisation et des mesures, et en présence du vent.

To my parents, JEAN and MARIE

It has been said that it's not so much the destination as the journey. It's also who you travel with! It's about a guide who puts you on the right path. I am very grateful to Claude and Tarek for having been that guide and for all that they offered me.

It is the wish of every student to be tutored by a big scientist. I had the chance of being the student of Claude, an "artist engineer". Thank you Claude for showing me how beautiful engineering can be. We have seen equations coming alive and making an airplane fly! I have enjoyed every single discussion we had, and I have learnt from you a lot! I will keep following your guidelines.

Thank you Tarek for having pushed me forward every time. It is to my best interest to have been guided by a determined "chief", so determined that you race with the sunlight to get to "your office". Thank you for having taught me relevant mathematics.

I would like to thank the members of the Jury for evaluating the present work. I am grateful to Prof. Pascal Morin and Prof. Abdelhamid Tayebi, for reviewing this manuscript and for their constructive criticism. My special thanks go to Prof. Brigitte D'Andrea-Novel whom I first met during my studies at Mines ParisTech and who is examining this thesis, along with Prof. Robert Mahony and Prof. Jean-Marc Moschetta.

"-AA: Bonne nouvelle, l'avion a volé, il s'est mis à tourner en rond". I can never forget those words spoken by the pilot and engineer André after the first successful flight. Thank you André for the very pleasant moments we spent in the lab as well as on the test field, and specially for having saved the airplane from crashing many times.

Playing football/cricket in the corridor, punching some boxing bag, watching football in the robotics room, practicing some magic tricks..That describes our daily work in the team. Thank you Ninad, Simone 1 and Lam-Hung, for your contributions to making our team the coolest and most noisy in the lab. I am also thankful to the "permanent" members, Minh-Duc, Andrew2 , and Guillaume.

Good friends don't let you do stupid things (alone). I am lucky to have a lot of those friends, Fernando, Pedro, Mohit, Davide, Miguel, Melissa, Howard, Henrique, Dino, Myriana, Jonathan, Samuel, Pablo, Danilo, you know what I mean, and since the majority are brazilians I have to say "Obrigado amigos"! The list of the people I knew at I3S is long, but I would like to mention Kallyna, Rodrigo, Pierre, Amina, Ophélie (my favourite magic fan), Vitali, Min, Arnaud, Anca, Cyprien, Gaël, and Kevin.

Conclusion

Bibliography

Notation

• R denotes the set of real numbers.

• The function sign(x) with x ∈ R is defined as:

sign(x) =          +1 if x > 0 0 if x = 0 -1 if x < 0
Euclidean vectors:

• Throughout this thesis, E 3 denotes the 3D Euclidean vector space and vectors in E 3 are denoted by bold letters. Inner and cross products are denoted by the symbols • and × respectively.

• For any vector u ∈ E 3 , Π u denotes the operator of projection on the plane orthogonal to u, i.e. ∀v ∈ E 3 ,

Π u v = v -(v • u)u.
• The norm of a vector ξ is denoted by |ξ|.

Vectors of coordinates:

• The coordinate vector of any ξ ∈ E 3 with respect to (w.r.t.) a specified frame is denoted by the ordinary letter ξ, and ξ i with i = 1, 2, 3 represents a coordinate, i.e. ξ = [ξ 1 , ξ 2 , ξ 3 ] ∈ R 3 .

• The norm of a vector ξ is denoted by |ξ|.

xix

• sat ∆ (y)(∆ > 0, y ∈ R n ) is the classical vector-valued saturation function sat ∆ (y) = min(1, ∆ |y| )y

• s at ∆ (y)(∆ > 0, y ∈ R n ) denotes a twice differentiable adaptation, with bounded derivatives, of the classical saturation function sat ∆ (y). More precisely s at ∆ (y) = α ∆ (|y|)y with α ∆ : [0, +∞) → (0, 1] a decreasing twice differentiable function such

that α ∆ (0) = 1, d dx α ∆ (0) = d 2 dx 2 α ∆ (0) = 0, α ∆ (x) ≤ ∆ x , lim x→∞ (α ∆ (x)x) = ∆. A typical example is α ∆ (x) = ∆
x tanh( x ∆ ). From these definitions s at ∆ (y) y when |y| is small and | s at ∆ (y)| ≤ ∆, ∀y.

• The scalar product of two coordinate vectors u and v is denoted by u v • S(.) is the skew symmetric matrix associated with the cross product of vectors of coordinates (i.e. S(u)v = u × v, ∀u, v ∈ R 3 ).

Automation in Aviation

T he idea of an airplane capable of flying autonomously without a pilot goes back to the early era of the invention of aerial vehicles. In this chapter, some historical events linked to the development of automatic flight controls are presented, followed by a brief description of today's research activities in the field of aerial robotics.

Historical examples of autonomous flights

In December 1903 (USA), the Wright brothers achieved the first powered manned flight.

They invented three-axis flight control which allowed to steer the aircraft. However, they quickly realized the difficulties of manually controlling their early aircraft which lacked inherent stability. In Europe, similar pioneering achievements in aeronautics can be mentioned, such as the works of Santos-Dumont in France, who succeeded in October 1906 in flying his "14-bis" design; this event was considered as the first heavier-than-air flight in Europe.

Prior to this, aviation pioneers attempted to fly models of unmanned aerial vehicles The american inventor and entrepreneur Elmer Ambrose Sperry, founder of the Sperry Gyroscope Company, worked on gyrocompass 2 technologies for maritime applications which led him to develop gyrostabilizers for airplanes. His first idea was to mount three gyros on the airplane's pitch, roll and yaw axes and couple them to the aircraft's control by servomotors to achieve automatic control of the orientation of the aircraft. He then decided to mount all three gyros on a single platform, thus allowing to measure the bank angle of the vehicle (the inclination of the aircraft with respect to a horizontal reference plane). On the 18th of June 1914, Sperry's son flew a gyrostabilized seaplane down the Seine River outside of Paris, where he stood up and let go of the controls while his mechanic climbed out on a wing. preplanned path defined by altitude and course, before commanding it to dive into its target. In order to be able to take away the pilot from the cockpit, Sperry worked on radio technologies and obtained, in 1917, the first patent for a radio-control system. This allowed to carry out flight tests with unmanned vehicles. On the 6th of March 1918, a

Curtiss-Sperry aerial torpedo catapulted into the air, successfully flew its preplanned path before diving into the water. However, not every flight was successful and many crashes were still happening, due to the lack of a reliable technology at that time.

Following these works, Charles Franklin Kettering formed the Dayton Wright Air-plane Company which took in charge the production of Liberty Eagle aerial torpedoes.

Kettering's design came to be known as the "bug", and represented numerous innovations in flight controls. Barometers where first used to control the climbing phase of the flight via a link to the elevators until a preset altitude was reached and where the system switches to turn the control to Sperry's gyrocompass to maintain level flight. Horizontal heading was also maintained by the gyrocompass which was also linked to small pneumatic valves for the actuation of the rudder. However this first concept led to consecutive deep stall situations and crashes. This was first interpreted as an over-actuation of the elevators, and adjustments were accordingly made, but this did not solve the problem.

Engineers went into searching for other solutions, until it was decided to reverse the roles of the barometer and the gyrocompass. The gyrocompass now controlled the pitch angle to achieve a climb out until a preset altitude was reached, then the barometers assumed complete control of level flight. With this new concept, and on the 4th of October 1918, the Liberty Eagle climbed normally and leveled off at an altitude of 200f t, then flew a preset distance of 1500f t, stopped its engine and dived to its planned target [START_REF] Newcome | Unmanned Aviation: A Brief History of Unmanned Aerial Vehicles[END_REF].

Many similar scenarios followed as the successful UAV flights of 1924 (England)

taken by the Royal Aircraft Establishment, and Germany's V-1 "buzzbombs" that were used in world war II, as well as the first French reconnaissance UAV, the R.20 that was developed in the late 1950s.

By 1941, UAVs were being equipped with television cameras and transmitters in order to enable their operation beyond the controller's sight. Then radar-guidance systems were also developed which enabled the control of UAVs in bad weather and at night. Positioning System (GPS). All these technologies made fully autonomous flight possible with higher precision. In the current day, pilots put airliners on autopilot within 20 minutes after takeoff, and it is even safe to perform automatic landings in bad weather. These small UAVs can perform various applications such as power line inspection, environmental monitoring, crop monitoring, aerial photography, search and rescue, and security surveillance. These tasks used to be carried by conventional helicopters and aircraft, but the use of smaller and relatively cheaper vehicles has dramatically decreased the cost of such operations. Moreover manned flights can be dangerous in some situations and even forbidden. Thus sending a small unmanned vehicle offers a practical and affordable solution. For instance, following the eruption of mount Usu in 2000, the Japanese government enlisted UAVs to fly reconnaissance missions to observe the advance of volcanic mud threatening surrounding villages, whilst it was too dangerous for manned aircraft to approach the active volcano.

Modern aerial robotics

Clearly, such small vehicles rely on traditional flight technology but as explained in [START_REF] Kumar | Springer Handbook of Robotics, chapter Modeling and Control of Aerial Robots[END_REF], their applications distinguish them from conventional aircraft by the fact that they need to interact autonomously with a complex dynamic three-dimensional environment.

This has complexified the sensing and control tasks of these vehicles and has boosted the research in robotics departments in universities as well as in industries and start-up companies. Many active research topics are related to aerial robotics, such as control theory for under-actuated vehicles, observer theory, path planning, computer vision, Simultaneous Localization And Mapping (SLAM), and multi-agent systems.

Categories of small UAVs

A variety of UAV designs exist, however the major categories of existing platforms are the rotor vehicles, fixed-wing vehicles and hybrid vehicles. Other types of vehicles are the airships (aerostats), flapping wing vehicles, ducted fans, and rockets.

Rotor vehicles

These vehicles rely on the presence of rotating aerodynamic blades to generate actuating forces and torques. They are particularly capable of achieving hovering flights, which requires generation of enough thrust to sustain their own weight and payload.

Perhaps the most famous platform of this type is the quadrotor which is equipped with four planar rotor-propellers, two of them rotating clockwise and the two others rotating counterclockwise. By commanding the speed of rotation of each motor individually it is possible to generate a desired thrust orthogonal to the plane of the rotors and desired torques along each of the three body-axis [START_REF] Hamel | Dynamic modelling and configuration stabilization for an x4-flyer[END_REF]. Other variants of these vehicles exist, with additional or less rotors like the hexacopter equipped with six rotors or the tricopter equipped with three rotors. In other less conventional configurations, the rotor-propellers are augmented with thrust-tilting capabilities [54][20], or are all mounted with different thrust orientations [START_REF] Rajappa | Modeling, control and design optimization for a fully-actuated hexarotor aerial vehicle with tilted propellers[END_REF] allowing the vehicle to behave as a fullyactuated vehicle.

Fixed-wing vehicles

Fixed-wing vehicles take advantage of profiled wings (permanently fixed to the aircraft's body) to generate aerodynamic lift forces, and achieve forward flights with thrust magnitudes less than the weight of the vehicle. The torque control is achieved through control surfaces at the trailing edges of the surfaces. Fixed-wing vehicles are the main subject of this thesis and will be detailed later on.

Categories of small UAVs

E

quations of motion of an aircraft can be derived using classical mechanics of rigid bodies. In this chapter, we define the state variables and actuation terms that are used to model the system, then we present the kinematics and dynamics equations that describe the evolution of these state variables via a set of ordinary differential equations. A simple modeling technique of the thrust produced by a motor-propeller is presented as well.

Airplane configuration and actuators

The basic components of a general aviation airplane are shown in figure 2.1. The fixed parts are the following:

• Fuselage: the center body of the airplane.

• Wing: this is the main wing of the airplane, and the main lift producing component.

• Horizontal tail: also named horizontal stabilizer.

• Vertical tail: also named vertical stabilizer. • Ailerons: These two control surfaces are mounted at the trailing edge of the wing near the wing tips at both the left and right sides. They usually deflect air in opposite directions and generate a rolling motion.

Airplane configuration and actuators

• Elevator: This control surface is located at the trailing edge of the horizontal tail, it is used to control the pitching motion.

• Rudder: This control surface is located at the back of the vertical tail, it is used to control the yaw motion.

The angular deflections of the ailerons, the elevator and the rudder are denoted by δ a , δ e and δ r respectively. The sign conventions of each of these angles are described later in chapter 3. We also define the vector of surface controls δ = [δ a , δ e , δ r ] .

Another actuator is the motor that produces the thrust vector T . Additional details on the expression of the thrust are reported in section 2.4.

Frames and state variables 2.2.1 Inertial frame and translational state variables

In this thesis, we assume a non-rotating earth, and a flat-Earth approximation. This is a valid approximation for vehicles flying near the surface of the Earth for short range navigation. Using these assumptions, the inertial frame denoted by I is considered to be an earth-fixed frame. The unitary vectors ı 0 ,  0 and k 0 are used to designate an orthonormal base for this frame, i.e., I = {O; ı 0 ,  0 , k 0 }. The choice of these vectors is done according to the North-East-Down (NED) convention: the origin O is a fixed point on the surface of the earth, the vectors ı 0 ,  0 are horizontal and point respectively in the north and east directions, and the vector k 0 is vertical and points downwards.

• The position of the center of mass G of the airplane is denoted by the position vector p = -→ OG.

• The inertial speed v of the center of mass of the aircraft is defined as the rate of change of p as seen in the inertial frame, i.e., v = d dt p. If the norm of the inertial speed |v| is always greater than a positive number, the direction (heading) of v can be defined as:

h = v |v| (2.1)
So that one can write v = |v|h.

The heading h is a unitary vector, hence it can be parameterized by two angles, the flight path angle γ, and the course angle χ (see figure 2.3), and can be written as:

h = cos γ cos χı 0 + cos γ sin χ 0 -sin γk 0 (2.2)

Body-fixed frame and rotational state variables

The aircraft is assumed to be a rigid body, and the body-fixed frame B is an aircraft-fixed frame. The three unitary vectors ı,  and k are used to designate an orthonormal base associated with this frame, i.e., B = {G; ı, , k}. We choose these vectors according to the following convention: the vector ı is in the plane of symmetry of the airplane along the fuselage, the vector  is parallel to the main wing, and the third vector k = ı ×  points in the direction below the aircraft belly (see figure 2.2).

• The orientation of the vehicle is represented by the body-fixed unitary vectors (ı, , k). It can also be represented by a rotation matrix, element of the Special Orthogonal group, R ∈ SO(3) : B → I, where the columns of R are the vectors of coordinates of (ı, , k) expressed in the inertial frame. An additional possibility is to use the three Euler angles parametrization known as the yaw ψ, pitch θ and roll φ, these angles can be related to the elements of the rotation matrix R, for more details see for example [START_REF] Stevens | Aircraft control and simulation[END_REF].

• The vector ω is the angular velocity of the body-fixed frame B w.r.t. to the inertial frame I, i.e., d dt (ı, , k) = ω × (ı, , k).

Wind frame and air-velocity

Let v w denote the wind velocity. The air-velocity v a is defined as the body's relative velocity to the wind i.e. v a = vv w . Assuming that the air-velocity is not zero, define the unitary vector h a = va |va| as the direction of air-velocity, and h a its coordinate vector in the body-fixed frame, i.e. h a represents the relative direction of the air-velocity w.r.t. the body. The vector of coordinates h a can be parametrized by two angles, and the commonly chosen ones are the angle-of-attack α and the sideslip angle β. These angles are shown in figure 2.4 and are defined as follows: A wind frame W is commonly defined when studying airplanes2 , with three associ-

h a = cos β(cos α ı + sin α k) + sin β  (2.3) h a (α, β) = [cos β cos α, sin β, cos β sin α] (2.4) α = atan( h a • k h a • ı ) (2.5) β = asin(h a • ) (2.6)
ated unitary vectors (h a , h ⊥,1 a , h ⊥,2 a ), where h ⊥,1
a is in the plane of symmetry (ı, k) and orthogonal to h a , and h ⊥,2 a completes the orthonormal base.

h a = v a |v a | (2.7) h ⊥,1 a =  × h a | × h a | (2.8) h ⊥,2 a = h a × h ⊥,1 a (2.9)

Kinematics and Dynamics

The aircraft is considered to be a rigid body 3 , so that the classical laws of rigid mechanics can be applied, and the study of the translational and rotational motions can be separated as follows:

Translation

Translation deals with the motion of the center of mass G of the airplane in the inertial frame I = {O; ı 0 ,  0 , k 0 }.

Kinematics of translational motion

The kinematics translation equation is the expression of the inertial speed, which is by definition the derivative of the position vector:

v = dp dt (2.10)
The second order kinematics equation is the expression of the linear acceleration a, which is the derivative of v in the same frame:

a = dv dt (2.11)
According to the definition of the heading h in (2.1), one can write v = |v|h and differentiate it to get an alternative expression for the linear acceleration:

a = d|v| dt h + |v|(ω h × h) (2.

12

)

ω h = h × dh dt (2.13)
where the vector ω h is the angular velocity associated with the unitary vector h.

Dynamics of translational motion

According to the flat-Earth approximation explained in section 2.2.1, the laws of conservation of momentum can be applied in the earth-fixed inertial frame I. Working in this (NED) frame, also leads to a constant gravity vector g = g 0 k 0 , with g 0 = 9.81m/s 2 as the standard gravitational acceleration.

Applying Newton's second law to the motion of the center of mass of the airplane, and considering its mass m to be constant, we get:

ma = mg + F a + T (2.14)
where F a is the total aerodynamic force 4 , and T is the thrust vector or the propulsive force.

Rotation

The study of the kinematics and dynamics of the rotational motion of a rigid body can be separated from its translational motion by taking as a reference point the center of mass G of the body i.e. the origin of the coordinate system associated with the body fixed frame B.

Kinematics of rotational motion

The derivative of the body fixed unitary vectors ı,  and k taken in the inertial frame, are related to the angular velocity vector ω according to the following kinematics equation:

d dt (ı, , k) = ω × (ı, , k) (2.15)
One can equivalently derive this equation by differentiating the rotation matrix R : B → I. Recalling that ω is the vector of coordinates of the angular velocity vector in the body frame, equation (2.15) is equivalent to the following:

Ṙ = RS(ω) (2.16)

Dynamics of rotational motion

Let H be the angular momentum vector taken about the center of mass. Applying the law of conservation of rotational momentum we can write:

dH dt = M a + M T (2.17) H = - ˆP∈body -→ GP × ( -→ GP × ω)dm = J .ω (2.18)
where M a is the total aerodynamic torque about G, and M T is the torque created by the propulsion system due to the displacement of the point of application of the thrust force from the center of mass. In the case of a propeller motor, additional aerodynamic and gyroscopic torques have to be included. J denotes the inertia operator at G. 4 The expression of the aerodynamic force Fa will be detailed in the next chapter. [START_REF] Farrell | Backstepping-based flight control with adaptive function approximation[END_REF]

Propulsion

Expressing the previous relations using vectors of coordinates in the body fixed frame offers advantages when working with rigid bodies. Indeed the inertia operator J becomes a constant matrix J ∈ R 3×3 called the inertia matrix, and equation (2.17) can be rewritten as:

J dω dt = -S(ω)Jω + M a + M T (2.19)

Propulsion

The propulsive mechanisms that produce the thrust force are the source of energy5 that powers the airplane through the air. Many types of propulsive systems exist. The reciprocating engine-propeller combination is the oldest propulsion device and still equips many general aviation aircraft. The invention of the more powerful jet engines created a milestone in the history of aerospace engineering, and made high-speed flight possible.

However, small scale fixed-wing airplanes are generally equipped with a motorpropeller combination, which due to the light weight of these vehicles generates a high thrust-to-weight ratio. The motor can be an internal combustion engine or more frequently an electrical brushless motor. For the applications of this thesis, we are interested in small airplanes and we will present the model of thrust generated by a propeller. Consider that the motor's shaft and therefore the propeller rotates at an angular velocity Ω = Ωı T , with ı T denoting the unitary vector parallel to the shaft of the motor.

The standard accepted thrust model, is a function of the angular velocity of the shaft, and is expressed as follows:

T = k Ω (Ω 2 )ı T (2.20)
with k Ω a positive constant.

More advanced models take into account the dependency of thrust on both the angular velocity Ω and the air-velocity v a by combining momentum conservation and blade element theories. These models also predict the presence of H-forces acting orthogonally to the motor's axis ı T , which makes the total thrust force expression three dimensional.

One can write for instance:

T = (k Ω (Ω 2 ) + ∆ T )ı T + T H (2.21)
where T H is orthogonal to ı T , and ∆ T is generally a negative term reducing the ideal thrust magnitude in equation (2.20). Such theories are relevant when studying nearhovering maneuvers for multicoptors as in [START_REF] Bangura | Aerodynamics of rotor blades for quadrotors[END_REF], [START_REF] Kai | Nonlinear feedback control of quadrotors exploiting first-order drag effects[END_REF] and helicopters as in [START_REF] Prouty | Helicopter Performance, Stability, and Control[END_REF]. However for the study of fixed wing aircraft traveling at relatively higher speeds, aerodynamic forces vary with the square of the airspeed |v a | and dominate first-order drag forces such as those that appear in T H . The standard expression in (2.20) is therefore an acceptable model in our case.

Let Θ be the point of application of the thrust force, and j prop ∈ R + the inertia of the motor/propeller about the shaft axis, then the total torque M T created by the motorpropeller about the center of mass of the aircraft can be written as follows:

M T = -→ GΘ × T -sign(Ω)k Ω Ω 2 ı T -j prop ω × Ω (2.22)
with k Ω a positive constant. The second term represents a resisting aerodynamic torque and the third term is the gyroscopic effect of the spinning propeller. If the motor is mounted so that the thrust direction passes through the center of mass of the aircraft, the first term in (2.22) is equal to zero. Furthermore, the effect of the last two terms is usually comparatively small, and the torque M T can be neglected in general because it can be easily compensated for by torque control actions.

Aerodynamics of Airplanes

A erodynamics is a sub-field of fluid mechanics that studies the interactions between solid objects and their surrounding air. This chapter revisits some basic concepts of this field, and introduces the principles that were exploited to invent heavier-than-air flying vehicles. The objective is to understand the behavior and modeling of aerodynamic forces and moments acting on airfoils, finite-wings and finally on an airplane as a whole. Assumptions that can be made for low airspeed flying regimes are explained as well. The aim of this chapter is to provide a control engineer with the necessary information needed to design control systems for airplanes. For additional material on aerodynamics see [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF][START_REF] Anderson | Introduction to flight[END_REF].

Sources and variations of aerodynamic forces and moments

Consider an experiment in which a fixed aerial vehicle is immersed in a steady and uniform three-dimensional flow of air. This rush of air generates a force on the body called the aerodynamic force. The force exerted by the airflow at each point on the surface of the vehicle is the effect of two forces per unit area, a shear stress and a pressure distribution:

• Friction shear stress f f : This force per unit area is tangential to the surface of the body. It represents the viscous friction that is the effect of opposition to motion between the fluid and the solid surface.

• Pressure f P : This force per unit area is normal to the surface. It is the result of repetitive "striking" of air molecules1 on the surface of the solid.

The total aerodynamic force exerted on the body is obtained by integrating these two vector quantities on the entire external surface (S) of the solid:

F a = ‹ (S) f f ds + ‹ (S) f P ds (3.1)
The point of application of this force is called the aerodynamic center which we denote by C. The associated aerodynamic torque, evaluated at the center of gravity G is:

M a = -→ GC × F a + M a,0 (3.2)
where M a,0 is a pure torque called the zero-lift moment. Then comes the question on how to predict the distribution of pressure and shear stresses. The solution involves a number of physical principles that lead to a set of coupled nonlinear partial differential equations called the Navier-Stokes equations. They can be solved using advanced numerical computation techniques, leading to the description of the evolution of the fluid in space and time. However, even to these days, some fluid properties cannot be totally represented in mathematical forms and their quantitative determination relies mostly on experiments. In the case of aerodynamic forces and torques, they can be measured in wind tunnel tests, which will be sufficient at least from a control engineer's point of view to determine their characteristics. Simplifying the analytical expressions of these forces is possible. Indeed, applying a dimensional analysis based on the Buckingham PI theorem shows that aerodynamic forces can be expressed as a function of a minimum number of independent dimensionless variables.

Steady aerodynamic coefficients

Each component F a,i of the aerodynamic force F a can be written in the following form (see [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF] p. 39):

F a,i = η a |v a | 2 C(R e , M, h a ) (3.3) η a = ρS 2 (3.4) R e = ρl|v a | µ (3.5) M = |v a | V s (3.6)
Where ρ is the freestream 2 density 3 , S is a characteristic surface commonly chosen as the wing planform area, l is a characteristic length commonly chosen as the wingspan, µ is the viscosity of air4 , and V s is the speed of propagation of pressure waves (sound) in air. The two dimensionless numbers R e and M are called respectively the Reynolds and the Mach numbers and are of prime importance in aerodynamics. C(•) is a dimensionless static coefficient that depends on the Reynolds and Mach numbers, as well as on the direction of the air-velocity. The complexity of the friction and pressure distribution in (3.1) is now hidden in dimensionless aerodynamic coefficients whose variations will be discussed further in next sections.

Dimensional analysis can be applied in the same manner to the aerodynamic torque M a and yields the following expression:

M a,i = η a l|v a | 2 C M (R e , M, h a ) (3.7)
With C M another dimensionless static coefficient.

It is well known that, for subsonic flows with a Mach number M smaller than 0.3, and for aerial vehicles flying at low Reynold numbers, as in the case of small UAVs, the aerodynamic coefficients are almost independent of the Mach and Reynolds numbers 5 .

This property is assumed in many works related to modeling aerodynamics for small UAVs such as [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF] and [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF], and is summarized in the following assumption:

Unsteady aerodynamic coefficients

The previous discussion treated the case of a fixed airplane in a steady airflow. However, spatial and temporal variations of the flow pattern around the object appear when maneuvering the airplane in rotational or accelerated motions, and the aerodynamic effects become different than those experienced in a steady flight.

It is therefore considered by aerodynamicists, that an acceptable mathematical approximation of aerodynamic coefficients would take the form C(α, β, ω, α, β, ω), where ω is the vector of coordinates of the angular velocity expressed in the body frame.

Additionally, the motion of control surfaces δ creates variations of forces on parts of the wings or tails. Unsteady motion is also associated with high control deflection rates, therefore transient control forces and torques are also functions of δ and δ.

Combining the preceding discussions, we can write as a summary the following approximations:

C(.) C(α, β, ω, δ, α, β, ω, δ) (3.8)
C M (.) C M (α, β, ω, δ, α, β, ω, δ) (3.9)

Aerodynamic forces of airplanes

The total aerodynamic force F a can be resolved into components along the body frame axis:

F a = F X + F Y + F Z (3.10) F X = -η a |v a | 2 C X (.)ı (3.11) F Y = -η a |v a | 2 C Y (.) (3.12)
F Z = -η a |v a | 2 C Z (.)k (3.13) (3.14)
The coefficient C X , C Y and C Z are called respectively the axial, side and normal forces coefficients 6 .

Another possibility is to express the components of F a in the wind frame:

F a = F D + F L + F C (3.15) F D = -η a |v a | 2 C D (.)h a (3.16) F L = η a |v a | 2 C L (.)h ⊥,1 a (3.17) F C = -η a |v a | 2 C C (.)h ⊥,2 a (3.18)
where F D is called the drag force, F L the lift force and F C the crosswind force. And C D , C L and C C are their corresponding non-dimensional coefficients. The axial and normal components were first used by Otto Lilienthal in 1889, however the Wright brothers preferred to use the drag and lift decomposition [START_REF] Anderson | Introduction to flight[END_REF]. In these days, the drag and lift aerodynamic coefficients are almost always used to describe aerodynamic forces on airplanes while the axial and normal components are still used in the case of bodies of revolution like rockets.

It is possible to work out an algebraic transformation between the aerodynamic coefficients of both systems as follows:

   C D C C C L    =    cos β cos α sin β cos β sin α -sin β cos α cos β -sin α sin β -sin α 0 cos α       C X C Y C Z    (3.19)
An airplane is designed to achieve a sufficient lift force during forward flight as opposite to non lifting bodies. This is primarily achieved by the design of the main wings, which strongly characterizes to the variations of the airplane's aerodynamic coefficients.

It is common to study first the aerodynamics of a wing's sections, and then to extend it to the case of an airplane. This is the purpose of the next sections.

Aerodynamics of airfoils

An airfoil is the cross-sectional shape of a wing. Figure 3.6 shows some nomenclature associated with airfoils. The measurements of aerodynamic forces and torques acting on an airfoil in a wind tunnel test are usually carried out on constant-chord wings with a span larger than the size of the wind tunnel's airflow. These measurements correspond to infinite wings, in the sense that any additional effects at the wing tips are not taken into account. The air-velocity is also considered to be parallel to the section, i.e. v a lies in the plane of the airfoil. With these assumptions we are brought to a planar study. Figure 3.7 shows the axis and force representations in the 2D case, note also that the direction of the airvelocity is now parametrized by the angle of attack α alone. In order to entirely describe the aerodynamic force F a , we will show measurements of the aerodynamic coefficients C L and C D for such infinite airfoils. The measured coefficients correspond to a steady case, therefore they depend only on α.

In order to produce a lifting force in the most efficient manner, airfoil designers try to find optimal shapes ensuring a compromise between high lift and low drag. The objective is to maximize the lift-to-drag ratio. One of such popular airfoils is the positively cambered NACA2412 shape 7 which is the wing section of most Cessna aircraft 8 . The following aerodynamic characteristics were obtained using the XFOIL tool 9 , which runs a data fusion process combining available experimental data with numerical predictions. 7 This airfoil has a camber of 2% of the chord 8 Some modern airfoils are superior to the NACA2412 airfoil in terms of efficiency, we mention the LS(1)-0417 airfoil designed by NASA during the 1970s [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF] 9 http://web.mit.edu/drela/Public/web/xfoil/ In figure 3.9, it can be seen that C L varies almost linearly with α and attains a peak value at an angle that we will denote by α s , then it drops drastically as α is increased further. This situation where the lift component drops is called the stall region and α s is called the stall angle. The stall phenomenon is caused by the flow separation on the upper surface of the airfoil that is more accentuated for α > α s .

The drag coefficient is almost constant below the stall angle α s , then increases dramatically beyond stall. The drag component of the aerodynamic force is mainly caused by viscosity of the fluid or what is known as skin friction drag which is obtained by integrating the shear stress over the surface of the body. However in the stall region, the flow separates early on the top surface causing a drop in the pressure of air in the wake, the difference of pressure between the air "in front" and at the "back" of the airfoil leads to an additional form drag which contributes in increasing the total net drag force. Note that the zero lift line (ZLL) 10 is not aligned with the chord line but it is inclined upwards and corresponds to an angle of attack of α 0 -2 • . In other words, for α = 0 a lifting force is still generated and we get C L (0) > 0. This is the case for most conventional cambered airfoils. However, throughout this thesis we choose the convention to align the body axis ı with the zero lift line as in figure 3.12 so as to preserve the property C L (0) = 0. This convention will be extended to the 3D case when working with finite wings and airplanes.

Classical modeling functions of airfoil's aerodynamic coefficients

The variations of the aerodynamic coefficients C L (α) and C D (α) are classically approximated by their first order Taylor expansions as follows (see [START_REF] Stengel | Flight Dynamics[END_REF]):

C L (α) C L,α α (3.

Aerodynamics of finite wings

In this section we extend the study to the aerodynamics of finite wings. The air velocity is still considered to be parallel to the wing's section. An important phenomenon happens at the wing tips where due to the pressure gradient, an air flow tends to appear from the region of high pressure under the wing to the relatively lower pressure on the top For a finite wing, the slope of the lift curve is smaller than that of the corresponding infinite airfoil with the same cross section. The aspect ratio AR of a wing is defined as follows:

AR = b 2 S (3.22)
Where b is the wing span and S is the planform area of the wing. Denote by C L,α the slope of the lift curve of the corresponding infinite airfoil. Then an approximation of the slope of the finite wing's lift curve is given by the following equation, as in [START_REF] Anderson | Introduction to flight[END_REF]:

C L,α = C L,α 1 + C L,α πeAR (3.23)
Where e is called the span efficiency factor and is typically equal to 0.95 for subsonic aircraft.

The drag force coefficient of the finite wing becomes the sum of the infinite airfoil's drag coefficient previously described and an induced drag term that can be considered to be equal to

C 2 L
πeAR . Writing the second order Taylor expansions of the lift and drag coefficients for a finite wing gives:

C L (α) C L,α α (3.24) C D (α) = C D0 + C 2 L πeAR C D0 + C 2 L,α πeAR α 2 (3.25)
Note again that those approximations are valid for values of α below the stall region (small value of α).

A numerical application on the wing whose section is supposed to be similar to the NACA2412 airfoil, and whose wingspan is equal to b = 1.5m and wing planform area is

equal to S = 0.36m 2 gives: C L (α) = 4.835α and C D (α) = 0.0095 + 1.218α 2 .

Extension to the case of airplanes

Now, the unit vector ı is chosen in the plane of symmetry of the airplane and parallel to the zero lift plane.

Steady forces with fixed control surfaces

We first consider the steady aerodynamic forces on the airplane. We also consider the control surfaces (ailerons, elevator and rudder) to be fixed in their neutral positions i.e. δ a = δ e = δ r = 0. Hence, the steady coefficients depend only on α and β.

Large values of the sideslip angle β makes the analysis very complex, however airplanes are usually designed to operate in balanced flights which consist of flying with almost zero sideslip angle. This has the advantage of exploiting the optimized aerodynamic characteristics that are obtained when the air-velocity is in the plane of symmetry of the airplane (highest lift and maximal lift-to-drag ratio). Another advantage of the balanced flight is the comfort it brings to the people on board, by zeroing lateral accelerations. The analysis of aerodynamic coefficients is limited to small values of β, for which

we can write C C C Y .
We first discuss the variations of the lift coefficient C L . For subsonic speeds, aerodynamic wind tunnel tests show that the lift coefficient of the airplane can be approximated to that of the complete wing by itself, including the portion of the wing that is masked by the fuselage (see [2, chapter6]). So the discussion about C L in section 3.2.2 concerning finite wings is still valid for the case of a fuselage-wing combination. Other components of the plane such as the horizontal tail can contribute to the lift but with smaller values. In general we can expect the variations of C L for an airplane to resemble to the case of a finite wing, and for small values of α and β, C L can still be written as:

C L (α, β) C L,α α (3.26)
with C L,α calculated according to 3.23.

Concerning the drag coefficient, it can be shown again that for small values of α and β, C D still takes the following form (see [2, chapter 6]):

C D (α, β) = C D0 + 1 πe os AR C 2 L C D0 + C D,α 2 α 2 (3.27) C D,α 2 = C 2 L,α πe os AR (3.28)
where the coefficient C D0 takes into account the skin friction drag of all the parts of the airplane. The term e os is called the Oswald efficiency factor, it is different from the efficiency factor e previously defined for finite wings, it is typically smaller than e (e os < e),

and is introduced to take into account parasitic drag forces due to aerodynamic interactions between the different parts.

The side force coefficient for a symmetric airplane can be approximated by its first order Taylor expansion as follows:

C C C Y C C,β β (3.29)
with C C,β > 0. The principal contributors to the side force are the vertical tail and the fuselage.

Total aerodynamic forces

It is commonly assumed that the unsteady effects can be modeled as linear perturbations to the steady aerodynamic coefficients previously described. In general for an aerodynamic coefficient C(.), one can write:

C(.) = C st (α, β) + ∆ C ( α, β, ω, ω, δ, δ) (3.30)
where C st (α, β) is a steady aerodynamic force coefficient.

Aerodynamic torques of airplanes

The components of the aerodynamic torque M a are commonly resolved in the body frame:

M a = M a,1 ı + M a,2  + M a,3 k, where M a,1 is called the rolling torque, M a,2 is
the pitching torque and M a,3 is the yawing torque. We define the torque aerodynamic coefficients C l , C m and C n according to the following:

M a,1 = η a b|v a | 2 C l (.) (3.31) M a,2 = η a c|v a | 2 C m (.) (3.32) M a,3 = η a b|v a | 2 C n (.) (3.33)
Note that the reference length used for the rolling and yawing torques is the wingspan b of the wing, but the chord length c is used in the case of the pitching torque 11 .

The purpose of the control surfaces δ = [δ a , δ e , δ r ] is to exploit the surrounding flow of air in order to modify the components of the generated torque. When a control surface deflects, it effectively changes the camber of the wing, which results in a change of the local force produced by that wing, hence affecting the aerodynamic torque about the center of mass of the aircraft.

The left and right ailerons are designed to deflect in opposite directions to generate opposite incremental forces on the two sides of the wing, which results primarily in a rolling moment M a,1 . The right aileron's deflection δ ar is positive when the control surface deflects upwards, and the left aileron's deflection δ al is positive when the control surface deflects downwards. The combined term δ a is commonly chosen as δ a = δar+δ al 2 , and with the previous conventions, a positive value of δ a is expected to generate a positive rolling moment. Similarly, a positive elevator deflection δ e corresponds to a control surface deflecting upwards and creates a positive pitching moment. A positive deflection of the rudder corresponds to a deflection of the control surface towards the right side of the airplane and is expected to create a positive yawing moment. (α,β,ω,δ,α,β,ω,δ). Some of these variables only affect a component of a torque lightly and can be neglected with respect to other more effective terms. Moreover, it is commonly considered that these coefficients can be modeled with acceptable accuracy using linear approximations.

Rolling moment coefficient

For small sideslip angles and small control surface deflections, the aerodynamic coefficient C l can be approximated as following [START_REF] Stevens | Aircraft control and simulation[END_REF]:

C l C l,β (α)β + C l,δa (α, β)δ a + C l,δr (α, β)δ r + b 2|v a | (C l,ω 1 (α)ω 1 + C l,ω 3 (α)ω 3 ) (3.34)
The sign of the term C l,β (α) is important for the lateral stability of the airplane, and is mainly dependent on the dihedral of the main wings of the aircraft. A positive dihedral angle gives one wing a positive angle of attack to the lateral component of the air-velocity (v a,2 ), and the other wing receives a similar negative angle of attack, thus generating a negative rolling moment when β is positive. Therefore C l,β (α) is generally negative for an airplane with dihedral, and this is a stabilizing configuration for the aircraft since it helps reducing the sideslip angle by rolling away from the relative wind.

A negative dihedral angle is called an anhedral and creates the opposite effect.

The ailerons generate torque through the term C l,δa . The rudder action also generates unwanted roll effect through the term C l,δr because the center of pressure of the rudder is normally above the level of the airplane's center of mass.

The term C l,ω 1 (α) represents the damping moment resisting the rolling motion of the aircraft and is a negative term. Note that in order for C l,ω 1 and C l,ω 3 to be dimensionless they were multiplied by the term b 2|va| in equation (3.34).

The values of the introduced terms in equation (3.34), are usually hard to estimate precisely, however some empirical formulas can be found in references like [START_REF] Nelson | Flight stability and automatic control[END_REF] or [START_REF] Stengel | Flight Dynamics[END_REF] and are used for pre-design stages and to evaluate the performance and stability of an aircraft.

Pitching moment coefficient

The aerodynamic coefficient C m can be approximated as follows [START_REF] Stevens | Aircraft control and simulation[END_REF]:

C m C m,0 + x cp C l,α c α + C m,δe (α)δ e + c 2|v a | (C m,ω 2 (α)ω 2 + C m, α(α) α) (3.

Yawing moment coefficient

The aerodynamic coefficient C n can be approximated as follows [START_REF] Stevens | Aircraft control and simulation[END_REF]:

C n C n,β (α)β + C n,δr (α, β)δ r + C n,δa (α, β)δ a + b 2|v a | (C n,ω 1 (α)ω 1 + C n,ω 3 (α)ω 3 ) (3.36)
C n,β (α) is strongly influenced by the design of the vertical tail. It provides directional stability when it is positive and helps the aircraft to align with the relative wind, hence reducing the sideslip angle.

The rudder action appears in the term C n,δr . And the negative term C n,ω 3 corresponds to the damping effect on the yawing motion.

Total aerodynamic torque

Using equations (3.31)-(3.33), we can write the total aerodynamic torque in the following compact form:

M a = η a |v a | 2 b    C l (.) c b C m (.) C n (.)    η a |v a | 2 b(C M,0 + C M,α α + C M, α(α) α + C M,β (α)β + A(α, β)δ + B(α)ω) (3.

Control Objectives

F

light control is the technology that assists a pilot in operating the aircraft, or even allows for a fully autonomous flight to take place without human intervention.

Various operating modes are encountered in practice. In this chapter, some control problems are discussed, then a survey of some existing advancements in guidance, linear and nonlinear control strategies is presented.

General control model

Using Equations (2.10), (2.14), (2.15), and (2.19), we can write the following set of ordinary differential equations describing the time evolution of the aircraft's motion:

               dp dt = v dv dt = g + Fa m + T m ı d dt (ı, , k) = ω × (ı, , k) J dω dt = -S(ω)Jω + Γ(δ) + Γ (4.1)
where in the second equation, the thrust direction is considered to be parallel to the body axis (ı T = ı). The total aerodynamic torque is written as the sum of a torque Γ(δ) generated by the control surfaces δ and a residual torque Γ , these terms can be deduced from equation (3.37). The torque due to propulsion is neglected.

In these equations, T and δ constitute control terms that should be used to control the state of the vehicle (p, v, (ı, , k), ω).

This typical model of a rigid airplane has six main degrees of freedom for the position and orientation1 . On the other hand, the dimension of the set of control terms (T, δ) is equal to four. This makes the aircraft an under-actuated system. However by looking at equations (4.1), it can be seen that the rotational dynamics are fully actuated via the control surface δ, while the translational dynamics are under-actuated via the monodimensional thrust T . Thus at the translational level, controllability results from the nonlinear coupling between the thrust and the orientation of the vehicle.

Objectives

A human pilot can fly an inherently stable airplane. This is done by directly requesting a commanded thrust T using a throttle. Control surface deflections δ are also under the control of the pilot who usually uses a joystick or a yoke to command the ailerons and the elevator, while the rudder is commanded by pedals. Small scale airplanes can also be piloted manually using radio control transmitters. However the stability properties of the airplane can sometimes be poor. Such as the case of highly maneuverable airplanes that are designed to gain maneuverability at the expense of their passive stability characteristics. The manual control of such vehicles can be demanding for the pilot in terms of concentration, a solution is to add embedded computers between the pilot and the actuators that will execute control algorithms to assist the pilot in keeping control of his aircraft. A step further to fully autonomous flight consists in giving the computer the complete task of controlling the vehicle. In the following we present the mathematical formulation of the attitude control, speed control, path following and trajectory tracking tasks.

Attitude and speed control systems

Instead of interpreting the pilot's commands directly as actuators setpoints, one can design a system in which the throttle is replaced by a speed setpoint v * ∈ R + . One may also use the joystick2 to specify a desired attitude equivalent to three desired body-fixed unitary vectors (ī, , k).

The speed setpoint v * can be considered as a desired magnitude of the inertial speed |v| or a desired magnitude of the air-velocity |v a |. This first task can be achieved by stabilizing the error |v| -v * (or alternatively |v a | -v * ) to zero using the thrust intensity T as a control variable.

The second task of stabilizing the attitude is achieved using the control variable δ to get a convergence of the aircraft's frame B = {G, ı, , k} to the desired mobile frame B = {G, ī, , k}.

Path Following

We consider now a fully automated case where the airplane has to approach and follow a pre-defined path without any time constraints. Consider a three-times differentiable curve C in 3D-space. The path following control objective consists in stabilizing the aircraft's speed |v| (or alternatively the airspeed |v a |) at a desired value v * , and in zeroing the "distance" between the vehicle's position and the desired geometric path.

These tasks are commonly decomposed into sub-problems: kinematical guidance and dynamical control.

Kinematical Guidance

Guidance involves body kinematic equations only and usually consists in determining the desired heading direction h * for the aircraft given the path to follow 3 . It is mostly independent of the aircraft characteristics.

Dynamical Control

This stage of control takes into account the specificities of the aircraft and of aerodynamic forces and torques that influence its dynamics. Dynamical control is in charge via the production of torques and thrust of making the actual aircraft heading direction h converge to the direction h * specified at the guidance level, and of stabilizing the aircraft velocity |v| or the airspeed |v a | to the desired value v * .

This separation between kinematical guidance and dynamical control is conceptually

Trajectory Tracking

Trajectory tracking refers to the problem of stabilizing a time-parametrized reference position trajectory. Let p r (t) denote a reference trajectory in E 3 , the control objective is to stabilize the position tracking error p = pp r about zero using the same control inputs T , and δ.

Trajectory tracking is best adapted to highly maneuverable rotor vehicles like quadrotors and helicopters, whose positions can be precisely and timely monitored near hover-

A review of guidance techniques

ing. Its applicability to fixed-wing aircraft, although possible is more limited because the position-timing issue for cruising vehicles is less essential than following a preplanned path with a given velocity needed for the production of strong air-lift forces on profiled wings.

Other control systems modes

Additional functions of automatic control systems require the stabilization of a part of the vehicle's state, like the altitude hold control system which monitors the altitude of the aircraft while allowing the pilot to freely steer the airplane in the horizontal plane.

The control systems may also run a heading hold mode, which means that the pilot specifies the desired heading vector h * without necessarily having a curve to follow, this can be done by parameterizing the direction h * using the course and the flight path angles, the objective would then be to achieve a constant course angle, or a constant climb or descent rate. Automatic control systems can also achieve automatic take-off and landing.

In some other control schemes, the pilot inputs are considered to be the angle of attack, and the stability axis roll rate 5 , while the sideslip angle is automatically regulated at zero.

In addition, if the airplane is marginally stable or unstable, stability augmentation systems (SAS) provide artificial stability to improve the flying qualities by adding for example artificial damping on the rotational motion of the aircraft. Augmentation system may also provide the pilot with a particular type of response to the control inputs, this is known as a control augmentation system (CAS).

A review of guidance techniques

Guidance concerns the kinematical part of the path-following problem, it constitutes the outer loop of an automated flight control system.

Many works study the case of 2D path following in the horizontal plane, by considering that the altitude is maintained or controlled by a seperate altitude controller.

Other studies investigate the possibilities to extend guidance strategies to the 3D case [START_REF] Breivik | Principles of guidance-based path following in 2d and 3d[END_REF] [11] [START_REF] Guilherme V Pelizer | Comparison of 3d path-following algorithms for unmanned aerial vehicles[END_REF]. The reference path may be fixed or, in more advanced applications, it can be attached to a frame that is moving with respect to the inertial frame [39] [38]. Basic studies consider the case of straight-line paths, and circular paths while others formulate the problem of following any regularly parametrized path [9] [11].

The guidance problem is generally formulated as the task of zeroing the distance between the vehicle's position and a virtual target on the reference path. This virtual target is generally considered as the orthogonal projection of the center of mass of the vehicle onto the path [START_REF] Samson | Path following and time-varying feedback stabilization of a wheeled mobile robot[END_REF]. Other variant characterization of the virtual target has been proposed to avoid projection non-uniqueness issues in some particular situations. For instance, the curvilinear abcissa of the virtual target may be considered as a free degree of freedom [START_REF] Breivik | Principles of guidance-based path following in 2d and 3d[END_REF] [11], or the virtual target is considered to be at the intersection of the path with a circle centered on the vehicle with a specified radius as in [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF]. Let us denote by p the distance from the virtual target to the vehicle.

The desired path angle γ * and course angle χ * are then calculated as a function of the components of p expressed in a frame attached to the reference path. The expressions of γ * and χ * are sometimes chosen as simple stabilizing saturated proportional controllers such as Line-of-sight strategies sometimes combined with pure pursuit tactics, or more generally as stabilizing nonlinear functions of p such as lookahead, or vector field strategies [START_REF] Derek R Nelson | Vector field path following for miniature air vehicles[END_REF]. In order to avoid the singularities associated with the definition of the path angle γ, the authors in [START_REF] Venanzio Cichella | Geometric 3d path-following control for a fixed-wing UAV on SO(3)[END_REF] use the Special Orthogonal group SO(3) in the formulation of the desired heading vector.

These desired path and course angles can be fed to the inner control loops as inputs, however in some guidance strategies, the desired values of their variations are also computed so that inner loops have to generate desired accelerations determined at the guidance level. Accordingly the desired path angle rate ω * γ and course angle rate ω * χ are computed, either as linear or nonlinear controllers (see for example [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF]) for the stabilization of the errors γ = γ -γ * and χ = χ -χ * .

A review of linear control techniques

Inner-loop controllers have historically, and to these days, been essentially designed on the basis of linearized modeled dynamic equations about so-called trim trajectories.

They are also reported in all major flight dynamics textbooks (as in [START_REF] Stevens | Aircraft control and simulation[END_REF], [START_REF] Nelson | Flight stability and automatic control[END_REF], [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF]), and for this reason, are often taken for granted in path following studies that focus only on the simpler generic guidance part of the problem such as the techniques presented in the previous section. In this section, the main principles of linear control techniques are briefly described.

Trim trajectories and linearization

Trim trajectories correspond to steady flight conditions, during which the components of the total forces and moments in the body-fixed frame are zero or constant. This corresponds to a constant steady state equilibrium along the trajectory. The most common trim conditions are the steady level flight, steady turning flight and the steady climb and sink flights.

The nonlinear dynamics of the system in (4.1) can be written equivalently in the form

Ẋ = f (X, U ) (4.2)
where X is the state vector and U = (T, δ a , δ e , δ r ) is the control input vector.

Let X * and U * be the solutions of the system along trim trajectories, that satisfy the system's equation Ẋ * = f (X * , U * ) 6 .

The nonlinear equations Ẋ = f (X, U ) can be linearized around the previously calculated equilibrium using the small-disturbance theory, where it is considered that the motion of the aircraft consists in small deviations about a steady flight. The state and control input vectors can be written as X = X * + ∆X, and U = U * + ∆U . And one can verify that

∆ Ẋ = A ∆X + B ∆U + H.O.T (4.3)
where the matrices A and B are defined as A = ∂f ∂X and B = ∂f ∂U evaluated at the equilibrium. These partial derivatives involve many aerodynamic derivatives as a result of differentiating aerodynamic coefficients. They are generally estimated from geometrical properties, from the variations of aerodynamic coefficients, or from perturbed motion of an aircraft during a flight test or a wind tunnel test. Higher order terms (H.O.T) can be neglected for small variations of ∆X and ∆U .

In order to simplify the problem, it is common to decouple the longitudinal and lateral dynamics. This means that the equations in ∆ Ẋ = A ∆X + B ∆U can be separated into two decoupled sets:

• Longitudinal dynamics described by:

∆ Ẋlon = A lon ∆X lon + B lon ∆U lon (4.4)
With for example

X lon = (v 1 , v 3 , θ, ω 2 )
where v 1 and v 3 are the first and third components of the speed in the body axis, and U lon = (T, δ e ) .

• Lateral dynamics described by:

∆ Ẋlat = A lat ∆X lat + B lat ∆U lat (4.5)
With for example, X lat = (v 2 , φ, ω 1 , ω 3 ) where v 2 is the second component of the speed in the body axis, and U lat = (δ a , δ r ) .

Open loop dynamics and stability of the uncontrolled motion

In order to study the stability of the uncontrolled motion, we can consider the longitudinal and lateral systems in equations (4.4) and (4.5), with the terms ∆U lon and ∆U lat put to zero. Then by computing the eigenvalues and eigenvectors of the system matrices A lon and A lat , we can get a clear idea about the inherent stability of the system and the different dynamical modes of motion close to the equilibrium. These modes are described by most of the flight stability studies as follows:

• The longitudinal modes associated with the system ∆ Ẋlon = A lon ∆X lon :

-Short-period mode: This is generally a short-period oscillatory mode associated with variations of the pitch angle θ (or equivalently the angle of attack α). This mode also characterizes the response of the pitch motion to the elevator input. Therefore it is of interest to design the airplane in such a way that this mode has a high convergence rate.

-Phugoid mode: This is a lightly-damped long-period oscillatory mode associated with the interplay of velocity and altitude (or path angle). This mode occurs slowly that even a pilot can easily correct for it manually.

• The lateral modes associated with the system ∆ Ẋlat = A lat ∆X lat :

-Roll mode: This is a highly convergent motion associated with the response of the rolling motion of the vehicle to the ailerons input.

-Spiral mode: This is a slowly convergent or divergent mode and is associated with the tendency of the vehicle to sideslip into a turn.

-Dutch-Roll mode: This a lightly damped oscillatory motion with a low frequency, and is associated with a rolling and yawing motion with some sideslipping. This mode can be excited by a rudder pulse.

The spiral and dutch-roll modes can be made stable by using appropriately sized vertical tail and wing dihedral.

MIMO linear control

The multi-input multi-output (MIMO) control design techniques, sometimes referred to as modern control techniques, are generally time-domain techniques that apply to state-space models as in equations (4.4) and (4.5).

Consider in general a linear system of the form ∆ Ẋ = A ∆X + B ∆U , if the system is controllable, a feedback control ∆U = -K∆X can be designed by assigning stable eigenvalues and eigenvectors to the resulting system ∆ Ẋ = (A -BK) ∆X. One then has to apply the control U = U * +∆U . In general the state ∆X is not fully available from measurements. It is possible however to design output feedback control for some measured output vector Y = C∆X, or design a state estimator that estimates the vector X using the output measurements Y and then apply the full state feedback. Additionally one can design dynamic regulators, such as integrators, in which case the state vector ∆X can be augmented with additional intermediate variables.

Instead of assigning desired eigenvalues and eigenvectors for the system, other popular approaches can be used to design a feedback control ∆U = -K∆X, such as linear quadratic regulators (LQR) which consist in selecting a performance criterion in the time domain in the form J = 0.5 ´∞ 0 (∆X Q∆X + ∆U R∆U )dt where Q is symmetric positive semi-definite, and R is symmetric positive definite, then finding the optimal gain K that minimizes J. Linear quadratic methods can be adapted to solve output feedback control problems, or they can be designed with full state feedback combined with a Kalman observer. This combination of a Kalman filter and a full state linear quadratic regulator is known as Linear-Quadratic-Gaussian (LQG) control and has certain guaranteed robustness properties as explained in [START_REF] Stevens | Aircraft control and simulation[END_REF]. Another popular MIMO control design is H ∞ control, this is a frequency domain technique that minimizes the maximum singular value of a transfer function matrix corresponding to a chosen error signal, this technique falls under the category of "robust control" design.

SISO linear control

Rather than studying high order systems, it might be convenient to close several simple feedback loops in succession and study lower order systems, this is known as hierarchical or nested loops control.

Additionally, it is simpler to deal with first order single-input single-output (SISO) systems, where each state variable is associated (controlled) with the most effective actuator, while the effect of the variations of other control inputs and states is neglected or treated as a disturbance.

The MIMO techniques that were described in the previous section were developed more recently and may theoretically guarantee the stability of the whole system which is not the case of the classical designs based on successive loops closure and SISO systems. However SISO techniques are still being employed in some cases, because of their simplicity.

For SISO systems, time-domain techniques can be easily applied to the resulting first order systems, or one can compute Laplace transfer functions and apply frequencydomain techniques such as root locus analysis to calculate control gains. -Sideslip control: The rudder (δ r ) is used to regulate the sideslip angle (β) and to achieve a balanced flight.

• Longitudinal control systems:

-Path angle control: In this outer loop, the pitch angle (θ) or the attack-angle (α) of the aircraft, are used as intermediate control variables to track the desired path angle (γ) specified at the guidance level.

-Pitch or angle of attack (AOA) control: In this inner loop, the elevators (δ e ) are used to control either the pitch angle (θ) or the attack-angle (α) of the aircraft.

-Velocity control: This objective can be achieved independently using the thrust (T ) of the aircraft as a control input and is known as auto-throttle control.

The linear techniques described so far apply on linearized models near specific flight conditions. In order to control the airplane on a larger flight envelope, the feedback gains should be adjusted according to the flight condition. This can be done using gain scheduling techniques which consists in interpolating gains between many values calculated at different equilibrium flight conditions.

A review of nonlinear control techniques

Nonlinear control techniques do not rely on linearization of the system's equations and may thus cover a wider flight envelope. They can yield enlarged domain of stability including high and rapidly varying angles of attack. Some techniques are based on feedback linearization that transforms the nonlinear system into an equivalent linear one using a state transformation. Or, the system can be divided into subsystems, and hierarchical controllers can be applied by closing nested control loops.

Nonlinear Dynamic Inversion (exact input-output linearization)

This method consists in differentiating the output vector successively until independent input control variables, allowing for the linearization of the output, are obtained (see for example [START_REF] Lane | Flight control design using non-linear inverse dynamics[END_REF]).

A balanced flight condition is generally assumed and where the rudder regulates the sideslip angle close to zero independently of the rest of the system. The output is then commonly taken as the inertial speed v = |v|h(χ, γ) or the air velocity v a , whose desired values and variations are specified at the guidance level.

Accordingly, the output vector that we will denote by z is the tracking error on the velocity and has a dimension of three. Differentiating it gives the error on the acceleration of the vehicle who is a function of the thrust T and the aerodynamics forces.

Assuming that the aerodynamic forces are independent of the angular rates and the control surface deflections, ż becomes a function of the thrust T and the attitude of the vehicle. Differentiating two more times allows to obtain three independent control variables u = ( T , δ a , δ e ) . These computations are generally cumbersome since the nonlinear equations of flight are complex, but it can be shown that the final equations can be written in the following form:

d dt    z 1 z 2 z 3    =    z 2 (x, T ) z 3 (x, T, Ṫ ) U (x, T, Ṫ , u)    (4.6)
where z 1 = z, x represents the original internal state variables, u = ( T , δ a , δ e ) is the input vector, and U (x, T, Ṫ , u) has the following form:

U (x, T, Ṫ , u) = Φ 1 (x, T, Ṫ ) + Φ 2 (x, T, Ṫ )u (4.7)
If Φ 2 (x, T, Ṫ ) and det(Φ 2 (x, T, Ṫ )) are nonzero, the application u → U (x, T, Ṫ , u) is bijective, this allows U to be viewed as a new control vector. System 4.6 is linear and controllable, and any MIMO linear control design can be applied and the stabilization of a desired reference velocity is possible. Finally, the values of the control terms u are computed as, u = Φ -1 2 (x, T, Ṫ ) (U -Φ 1 (x, T, Ṫ )). In order to complete such a nonlinear control design, one has to examine the stability of the internal variables (apart from the controlled output variables). One way of examining the boundedness of internal variables is to set the output error identically to zero, and examine the stability of the zero dynamics.

However this technique raises the following issues:

• The second derivative of the thrust T is generally not available as a control variable, indeed one has to incorporate the motor's dynamics into the control design, which complicates the task. Additionally, the previous control design requires the knowledge of T and Ṫ which are not measured and are hard to estimate.

• The condition det(Φ 2 (x, T, Ṫ )) = 0 is not guaranteed along the trajectory.

• This nonlinear dynamic inversion (NDI) controller requires the full knowledge of all the nonlinear dynamics of the aircraft. In practice, many of the aerodynamic coefficients may not be available, and the robustness of such a controller is not guaranteed.

Backstepping Control

This design procedure is in some ways similar to the nonlinear dynamic inversion technique presented in the previous section, however, instead of linearizing the entire system, a stabilizing feedback control is constructed for each subsystem by considering at each step a state variable as a virtual control input. Lyapunov functions are used to prove the stability of the global system. This technique falls in the category of hierarchical controllers, such as the control designs that will be developed in this thesis. It Each of these tasks, forms a nonlinear system, whose control terms can be designed using backstepping techniques.

In this section we will present an example for stabilizing the flight path angle γ as it is commonly done using the classical approximation of the lift aerodynamic coefficient C L .

We will also use this example to point out the difficulties encountered when using such classical models, in order to motivate the proposed modeling and design approaches that will be proposed in later chapters.

Assuming first that the task of zeroing the sideslip angle is achieved, the simplest path angle control design considers that the roll angle is close to zero. These assumptions reduce the system equations to a two dimensional longitudinal situation, therefore neglecting lateral effects. Accordingly, we have

α θ -γ (4.8)
In the absence of wind, it can be shown using the translational dynamics equation and the definition of γ that:

γ = T sin α m|v| + η a |v| 2 C L m|v| - g 0 cos γ |v| (4.9) = T sin α m|v| + η a |v|C L,α α m - g 0 cos γ |v| (4.10)
Where we have used the classical modeling expression of the lift aerodynamic coef-

ficient, C L = C L,α α.
In equation 4.10, α or equivalently θ = α + γ should play the role of a virtual control input. As can be seen it is complicated to work out an explicit expression for the control term due to the appearance of both terms α and sin α in equation 4.10. In some works, the term T sin α is considered to be negligible comparing to the lift force and therefore is omitted from the equation (as in [START_REF] Peng | Backstepping control for longitudinal dynamics of hypersonic flight vehicle[END_REF]). In other works, it is assumed a low angleof-attack flying regime and sin α is replaced by its first order Taylor approximation, i.e.

sin α α while T is replaced by its desired value T * whose computation is supposed to be known from the airspeed regulation subsystem. Let's consider here the former case, and replace α by θ -γ. The system can be written in a strict feedback form as follows:

γ = a 1 |v|(θ -γ) - g 0 cos γ |v| (4.11) θ = ω 2 (4.12) ω2 = Φ 1 (|v|, α, α, ω 2 ) + Φ 2 (|v|, α)δ e (4. 13
)
With

a 1 = ηaC L,α m
a positive constant, and Φ 2 a positive function. Backstepping control techniques can be used to solve this problem (see for example [START_REF] Peng | Backstepping control for longitudinal dynamics of hypersonic flight vehicle[END_REF]), they consist of a sequence of design problems for low order systems, and can be briefly described as follows:

• Extending this idea to the 3D case is not straightforward. Indeed the assumption that the roll angle is zero doesn't hold, and equation 4.8 is not valid. The angle α becomes a highly nonlinear trigonometric function of the attitude and the speed direction, i.e. α = α(φ, θ, ψ, γ, χ). Therefore in general, it is considered that α is an intermediate control input for the stabilization of γ (as in [START_REF] Farrell | Backstepping-based flight control with adaptive function approximation[END_REF]). This however raises another issue.

Indeed keeping the assumptions on a zero sideslip angle and in the absence of wind, the expression of γ in a three dimensional case takes the following form:

γ = T sin α m|v| cos 2 γ -sin 2 φ cos γ + η a |v|C L,α α m cos φ - g 0 cos γ |v| (4.14)
Therefore computing a value of a desired angle of attack α * will depend on the roll angle φ. The desired attitude of the vehicle is hard to compute and becomes itself a function of the current orientation of the vehicle, which makes the problem of the attitude control loop ill-posed.

Some works consider the case where the pilot input commands consist of a desired angle of attack (or pitch angular rate), a stability-axis roll rate, and a desired sideslip angle.

For example in [START_REF] Escande | Nonlinear dynamic inversion and LQ techniques[END_REF], this problem was solved by combining time-scale separation, NDI techniques for the computation of feedforward terms for each loop, and LQ techniques for the computation of feedback terms. A similar work is found in [START_REF] Snell | Nonlinear inversion flight control for a supermaneuverable aircraft[END_REF] and in [START_REF] Wang | Robust nonlinear flight control of a high-performance aircraft[END_REF] where dynamic inversion is combined with stochastic robust design. Nonlinear adaptive control techniques were also investigated in order to take into account uncertainties and parametric changes in the nonlinear model of the aircraft as in [START_REF] Sonneveldt | Nonlinear adaptive trajectory control applied to an f-16 model[END_REF] and [START_REF] Bouadi | Contribution to Flight Control Law Design and Aircraft Trajectory Tracking[END_REF]. In other works such as [START_REF] Marcos | Development of linear parameter varying models for aircraft[END_REF] and [START_REF] Duan | A velocity-based lpv modeling and control framework for airbreathing hypersonic vehicle[END_REF], nonlinearity is dealt with by applying Linear Parameter Varying (LPV) modeling and control.

As was explained in section 4.2, trajectory tracking applications were mainly developed for vehicles operating near hovering modes. There is indeed a large amount of research and results concerning trajectory tracking solutions using rotor vehicles (mainly quadrotors), see for example [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], [START_REF] Lee | Geometric tracking control of a quadrotor uav on se (3)[END_REF] and [START_REF] Hamel | Dynamic modelling and configuration stabilization for an x4-flyer[END_REF]. However trajectory tracking applications for fixed-wing vehicles are less common and few works can be found in this area [START_REF] Sonneveldt | Nonlinear adaptive trajectory control applied to an f-16 model[END_REF] [51]

[21] [START_REF] Núñez | Towards 4d trajectory tracking for transport aircraft[END_REF]. The robotics community has been mainly involved in projects related to small rotor vehicles (such as the quadrotor) and the well accepted hierarchical control design based on the thrust direction control paradigm has gained popularity due to its simplicity and practical performance (see [START_REF] Hua | Introduction to feedback control of underactuated vtol vehicles: A review of basic control design ideas and principles[END_REF]). On the other hand, control design for fixedwing aircraft is more than a century years old and has been extensively developed in the aeronautics community, or more specifically in what is known as the field of Guidance, Navigation and Control (GNC). The traditional and adopted techniques for the control of aircraft are mainly based on the linearization of the aircraft dynamics around trim trajectories. Linear control methods are then applied to the resulting system. Other nonlinear techniques were also discussed in section 4.5.

Recently, the control community had to deal with projects related to hybrid vehicles 1(see chapter 1) that can perform stationary flights like multicopters and helicopters, and also fly at high speeds while taking advantage of the lifting properties of profiled wings like airplanes. The control of these vehicles was initially addressed by combining the two existing approaches of control and by switching between these control laws.

Hierarchical control design

This transition between the two modes is delicate and does not handle the high nonlinearity of the aerodynamic models at high angles of attack. The question that arises is whether it is possible to design similar control architectures that can deal with all classes of vehicles using a unique generic control law that does not involve a switching policy. Previous research works in [START_REF] Pucci | Nonlinear control of PVTOL vehicles subjected to drag and lift[END_REF][START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] investigated the possibility of adaptation of the hierarchical control design approach based on the thrust direction control paradigm to bodies subjected to aerodynamic forces such as airplanes. It was shown that a class of aerodynamic models allows indeed the development of such feedback control theories with stability and convergence analysis.

More precisely, a common control architecture based on a hierarchical control decomposition can be applied to a large class of vehicles such as rotor vehicles (e.g. quadrotors), axisymmetric vehicles (e.g. rockets) [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF], fixed wing vehicles [START_REF] Hua | A novel approach to the automatic control of scale model airplanes[END_REF] [23], and hybrid vehicles. A generic control model that can serve as a basis to design autopilots for trajectory tracking and path following is also common to these vehicles. A synthesis of these control approaches in the case of fixed-wing airplanes is presented next.

In this chapter, the concept of hierarchical control design is presented. Then it is shown how the fully actuated attitude dynamics allows to consider a lower-order model for control design, independently of the actuation configuration of a specific vehicle. Finally, a simplified "model for control" of aerodynamic forces is presented and compared with classical aerodynamic modeling techniques.

Hierarchical control design

Consider again the control model of section 4.1:

dp dt = v (5.1) dv dt = g + F a m + T m ı (5.2) d dt (ı, , k) = ω × (ı, , k) (5.3) J dω dt = -S(ω)Jω + Γ(δ) + Γ (5.4)
Hierarchical control consists in decomposing the system into a cascade of subsystems of lower orders. For instance, in this case, the system (5.1)-(5.4) can be decomposed into two main subsystems as follows:

• Position control: An outer-loop for the subsystem composed of equations (5.1) and

(5.2), where the position p and/or the velocity v are controlled using the thrust T and the orientation of the vehicle as virtual intermediate inputs. This consists in determining a desired thrust, and a desired orientation of the vehicle equivalent to the determination of three body-fixed unitary vectors that we denote by (ī, , k).

• Attitude control: Inner control loops concerning equations (5.3) and (5.4), where first in equation ( 5.3) the angular velocity is used as an intermediate control input that we denote by ω * to make the orientation of the vehicle track the desired orientation specified at the position control level. Then in equation (5.4), the control surfaces δ should generate a torque control term that makes the actual angular velocity track the desired angular velocity previously defined.

For convergence analysis, either backstepping (based on classical Lyapunov-type approaches) or high-gain techniques (based on singular perturbation theory) can be applied. In this latter case, it is important to keep in mind that the rate of the convergence of inner-loops must be higher than that of outer-loops, so that an assumption of timescale separation would be valid. In other words, each loop should converge relatively faster than its outer-loop so that the intermediate control variables can be considered to be applied almost instantaneously. This is what gives the position control the name of "slow" outer-loop, and the attitude control the name of "fast" inner-loop.

This control scheme has many advantages, let us mention some of them:

• Working with systems of low orders (mainly first and second order systems) facilitates practical implementation, tuning and failure diagnosis procedures.

• Position and velocity estimation relies on low frequency measurement rates, like global positioning systems (GPS), pitot tubes and cameras, while for attitude estimation higher frequency measurement rates are available from onboard Inertial Measurement Units (IMU) for example. This means that the rate of convergence for the estimators of translation variables is more limited than that of the attitude estimators. This has direct implication on control design, and justifies the separation of the system's dynamics to slow outer-loops and fast inner-loops.

Control model

Designing control laws along with convergence analysis is done on the basis of a mathematical model. The chosen model should be close enough to the actual physical system and yet simple enough to lend itself to analysis. In this section we explain the choices that are made in order to simplify the control model. Following these assumptions, we choose the following control model:

         dp dt = v dv dt = g + Fa m + T m ı d dt (ı, , k) = ω × (ı, , k) (5.5)
With T and ω taken as control inputs.

It is interesting to note that these equations are representative of a large class of aerial vehicles, by being independent of the actuation specificities of a vehicle. The resulting system can be used to design general hierarchical control principles that apply to a large number of vehicles.

Modeling aerodynamic forces for control

System (5.5) has to be complemented with a model for the resultant aerodynamic force F a . Chapter 3 presented a physical description of aerodynamic forces, and their classical approximation as linear functions of the attack and sideslip angles. In this section, an alternative nonlinear model of aerodynamic forces will be presented which is inspired from previous works in [START_REF] Hua | A novel approach to the automatic control of scale model airplanes[END_REF], [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] and [START_REF] Pucci | Flight dynamics and control in relation to stall[END_REF].

The model that we propose to use is,

F a = -η a |v a |(c 0 v a,1 ı + c0 v a,2  + c0 v a,3 k) (5.6)
With c 0 , c0 and c0 denoting positive coefficients.

This model is compatible with relations 3.10-3.18, and corresponds to, C X = c 0 h a,1 = c 0 cos α cos β, C Y = c0 h a,2 = c0 sin β and C Z = c0 h a,3 = c0 sin α cos β.

Let c 1 be a positive constant coefficient such that, c0 = c 0 + 2c 1 . In the case of zero sideslip (β = 0), one easily verifies using equation (3.19) that: is that the aerodynamic coefficients do not grow unbounded when the angles α and β get large, but they vary periodically with the relative motion of air. For instance if the drag coefficient c 0 were equal to zero then, in the case of zero sideslip angle, the resultant aerodynamic force would be orthogonal to the zero-lift plane with an amplitude proportional to sin α|v a | 2 .

C D (α, 0) = c 0 + 2c 1 sin 2 (α) (5.
Note that this model represents forces acting on bodies with symmetries. More precisely, (ı, ), (, k) and (ı, k) should all be planes of symmetry for the expression 5.6 to hold for all values of v a,1 , v a,2 and v a,3 . This can be seen by noticing that F a (-v a ) = -F a (v a ). One can also notice that the equivalent drag and lift coefficients satisfy the properties C D (-α, 0) = C D (α, 0), and C L (-α, 0) = -C L (α, 0). This of course is not exact in the case of an actual aircraft whose geometry presents only a vertical plane of symmetry which is (ı, k). However, we expect that the model approximates at best the aerodynamic forces for h a,1 > 0 and h a,3 > 0, i.e. for forward flights and positive angle-of-attack regimes, which cover almost all the intended maneuvers of an aircraft. This model is also coherent with experimental data performed on a variety of wing profiles and axisymmetric bodies as explained in [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF] and [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. However, for lift-optimized wing profiles, this model fails to account for stall phenomena occurring at large attack angles. Modeling a large flight envelope is still possible by combining different classes of modeling functions, the interested reader is referred to [START_REF] Pucci | Towards a unified approach for the control of aerial vehicles[END_REF], and [START_REF] Pucci | Flight dynamics and control in relation to stall[END_REF]. However, for the applications in this thesis, we are mostly interested in controlling the aircraft in the linear regime without dealing with the stall region.

It is clear that in accordance with assumption 3.1, the aerodynamic coefficients in expression 5.6 do not vary with the Reynolds and Mach numbers. Additionally this model implicitly assumes that the effects of rotational and unsteady motions and of the control surfaces are neglected. This latter assumption has the advantage in a mathematical point of view of making system 5.5 triangular, so that a hierarchical controller can be designed without dealing with zero-dynamics. On the other hand, whether the unsteady effects can really be neglected without destabilizing the system is a matter of investigation that has to be verified in simulations and practical experiments.

Case of spherical bodies

For a spherical body, and considering the symmetry of the problem along all the body axis, the aerodynamic force do not depend on the body's orientation and is reduced to its drag component. The aerodynamic force is thus parallel to v a . The general model (5.6) can be adapted to this case by zeroing the lift coefficient, i.e. taking c 1 = 0, and since the coefficients c 0 , c0 and c0 should all be equal one can write,

F a = -η a |v a |(c 0 v a,1 ı + c 0 v a,2  + c 0 v a,3 k) = -η a |v a |c 0 v a (5.9)
with c 0 a positive and constant drag coefficient.

Case of airplanes

The shape of the airplane does not present any axial symmetry, and the coefficients c 0 , c0 and c0 are in general all different. The general aerodynamic model in (5.6), can be written equivalently as follows,

F a = -η a |v a |( c0 v a -2c 1 v a,1 ı -( c0 -c0 )v a,2 ) (5.10)
Where we used again the definition c0 = c 0 + 2c 1 .

As a conclusion, despite the simplicity of these models, they are physically pertinent and locally similar to the classical known approximations of aerodynamic forces2 . They also allow for elegant stability analysis and control design as will be seen in the next chapters. Similar models were recently proposed in [START_REF] Leandro R Lustosa | Longitudinal study of a tilt-body vehicle: modeling, control and stability analysis[END_REF] and [START_REF] Leandro | The φ-theory approach to flight control design of hybrid vehicles[END_REF], where aerodynamic forces and torques are expressed in the body frame of the vehicle and are globally nonsingular (i.e. they are not written as functions of attack and sideslip angles), the authors have also noticed the advantages of employing these models for simulation and stability analysis.

6

Attitude Control

A

ttitude stabilization corresponds to designing the fast inner-loop of the hierarchical control architecture. Once this loop is active and achieving its objective, outer-loops can be designed separately.

In this chapter, we first derive a control law for the angular velocity as an intermediate control variable along with a convergence analysis. Then we show how to track this desired angular velocity with the control surfaces of a conventional airplane. Finally, an adaptation of the previous control to the case of a flying wing is shown.

Problem statement

The attitude control problem is associated with the following subsystem taken from section 5.1:

d dt (ı, , k) = ω × (ı, , k) (6.1) J dω dt = -S(ω)Jω + Γ(δ) + Γ (6.
2)

The objective is to achieve a convergence of the aircraft's frame B = {G; ı, , k} to a desired time-varying frame B = {G; ī, , k}.

Attitude control design

According to section 5.2, we divide this subsystem further to two subsystems:

• An orientation control associated with equation 6.1, and where a desired angular velocity ω * plays the role of the intermediate control variable that ensures the convergence B → B

• A torque control problem where the control surfaces δ generate the required torque that ensures the convergence ω → ω * , where ω * is the vector of coordinates of ω * expressed in the body frame.

Attitude control design

Consider a time-varying desired body frame B = {G; ī(t), (t), k(t)}. It is also assumed that ī(t), (t) and k(t) vary smoothly with time so that their derivatives i(t), (t), and k(t) are well-defined.

The angular velocities of ī and  are defined respectively as:

ω ī = ī × i (6.3) ω  =  ×  (6.4)
One can then deduce the angular velocity of the frame B as:

ω = ω ī + (ī • ω )ī (6.5) = ω  + ( • ω ī) Proposition 6.1.
Assuming that the frame B and its angular velocity ω are well-defined, an angular velocity control that almost globally asymptotically (locally exponentially) stabilizes

B = B is ω = ω + k ω (t)((ı × ī) + ( × ) + (k × k)) (6.6)
with k ω (t) > > 0.

Proof.

Let (u θ , θ) denote the axis-angle representation of the rotation between the frames B and B. Consider the following candidate Lyapunov function:

V θ = 0.5 tan 2 ( θ 2 ) ≥ 0
Taking the derivative of V θ gives

Vθ = 0.5 tan ( θ 2 ) θ cos 2 ( θ 2 )
By definition we have θ = ( ωω) • u θ . Replacing ω by the control law of equa-

tion (6.6) yields θ = -k ω (t)((ı × ī) + ( × ) + (k × k)) • u θ .
Next, we show that (ı ×ī) + ( × ) + (k × k) = 2 sin θ u θ . According to Rodrigues' formula for the rotated vectors we have,

ī = cos θ ı + sin θ (u θ × ı) + (1 -cos θ)(u θ • ı) u θ (6.7)  = cos θ  + sin θ (u θ × ) + (1 -cos θ)(u θ • ) u θ (6.8) k = cos θ k + sin θ (u θ × k) + (1 -cos θ)(u θ • k) u θ (6.9)
From 6.7, one deduces that:

ı × ī = sin θ (ı × (u θ × ı)) + (1 -cos θ)(u θ • ı) (ı × u θ ) = sin θ u θ -sin θ(ı • u θ ) ı + (1 -cos θ)(u θ • ı) (ı × u θ ) (6.10)
Similarly,

 ×  = sin θ u θ -sin θ( • u θ )  + (1 -cos θ)(u θ • ) ( × u θ ) (6.11)
and,

k × k = sin θ u θ -sin θ(k • u θ ) k + (1 -cos θ)(u θ • k) (k × u θ ) (6.12)
Summing (6.10), (6.11) and (6.12) gives:

(ı × ī) + ( × ) + (k × k) = 3 sin θ u θ -sin θ u θ + (1 -cos θ) (u θ × u θ ) = 2 sin θ u θ (6.13)
Therefore θ = -2k ω (t) sin θ. Replacing θ in the expression of Vθ gives:

Vθ = -4k ω (t)V θ ≤ 0
Therefore V θ converges exponentially to zero, and the almost global asymptotic stability of B = B follows from the definition of V θ, with the domain of attraction { θ(0) : θ(0) = π}. Also note that orientations such that θ = π are unstable equilibria.

Torque control

Denote by ω * the desired angular velocity determined previously. The objective is to calculate the expression of control surfaces δ that achieves the convergence of ω to ω * .

We write equation (6.2) again for convenience, J dω dt = -S(ω)Jω + Γ(δ) + Γ (6.14)

Recall that Γ(δ) is the torque due to the control surfaces, and Γ is the residual torque, and both of these terms can be deduced from equation (3.37).

Proposition 6.2. The following expression of the desired torque Γ(δ) * produced by the control surfaces makes

ω track ω * Γ(δ) * = J ω * + S(ω)Jω * -Γ -k γ J(ω -ω * ) (6.15)
with k γ a positive gain.

Proof. Let ω = ω -ω * denote the error on the angular velocity. One verifies using the control expression in equation (6.15) that

J ω = -S(ω)J ω -k γ J ω (6.16)
Consider the following cost function

V ω = |J ω| 2 ≥ 0 (6.17)
which is equal to zero only when ω is zero. Taking the derivative of V ω and using equation (6.16) gives

Vω = 2(J ω) (J ω) = -2(J ω) (S(ω)J ω + k γ J ω) = -2k γ |J ω| 2 = -2k γ V ω ≤ 0
From which results the global exponential convergence of V ω and hence of ω to zero.

In practice however, it is sometimes hard to get a good estimate of the inertia matrix J. The cancellation of the residual torque Γ produced by the non-moving surfaces is also challenging and involves identification of many aerodynamic coefficients and parameters, as can be seen in equations (3.34)-(3.36). However using a simple proportional feedback with a high gain is sufficient in practice to make ω track ω * . Indeed the fully actuated nature of the attitude dynamics and the availability of the estimate of the angular velocity at a high rate makes the use of high proportional gains possible1 , allowing us to take advantage of the robustness associated with the use of high proportional gains.

The modified control law becomes

Γ(δ) * = -K γ (ω -ω * ) (6.18)
with K γ a high gain matrix, which can be chosen as a diagonal matrix for example if one assumes that J is almost diagonal.

It should be noted that for a regular airplane, the neglected torque Γ has a tendency to achieve passive stability by its restoring and damping terms as explained in section 3.3. Additionally the maximum active torque Γ(δ) that can be achieved is limited by the mechanical stops of the actuators (the possible range of the components of δ) and the flying regime (the value of |v a |). This may limit the maximum achievable angular velocity ω * .

The last step of the attitude control design is to compute the angles of deflection δ of the control surfaces that achieve the desired control torque Γ(δ) * . In view of (3.37), the relation used to model the production of the control torque takes the following form,

Γ(δ) |v a | 2 Ā(α, β)δ (6.19)
The matrix-valued function Ā(α, β) is aircraft specific, it depends in particular on the placement of the control surfaces with respect to the aircraft center of mass and on their dimensions. Typically around nominal air velocities and for small attack and side-slip angles, Ā can be approximated as a constant matrix which in view of equations (3.34)-(3.37) can be written as

Ā = η a b    C l,δa 0 C l,δr 0 c b C m,δe 0
C n,δa 0 C n,δr    (6.20)
Inverting relation (6.19) suggests to make δ track the desired angles of deflection δ *

given by

δ * = 1 |v a | 2 Ā-1 Γ * = - 1 |v a | 2 Ā-1 K γ (ω -ω * ) (6.21)
The matrix Ā may not be known precisely in practice. One may then calculate the desired deflection angles according to the following relation,

δ * = - 1 |v a | 2 K δ (ω -ω * ) (6.22)
with K δ a high gain matrix. If one assumes further that Ā is almost diagonal, K δ can be chosen as a high gain diagonal matrix.

In general the actuators have their own dynamics that should be included in the analysis in order to make δ track δ * . However small scale airplanes are generally equipped with electrical servomotors that accept setpoints in angular deflections and whose proper dynamics can be supposed to be sufficiently fast so that δ * can be considered to be instantaneously tracked. The maximum deflection rate δ is surely limited by the technology of the actuators and this is usually determined in the specifications of the motors. This in turn has implications on the maximum rate of change of the angular velocity ω * that can be achieved. However, these issues are not addressed in this thesis.

Adaptation to two-axis pitch-roll autopilots

In the preceding section, we considered the case of three-axis autopilots monitoring angles of the ailerons, elevator and rudder. Some small scale vehicles such as flying-wings involve ailerons and tail elevator actions only. The autopilot associated with their configuration is called a two-axis autopilot. The possibility of designing such less sophisticated autopilots is discussed next.

Pitch-roll attitude control

The desired attitude of the the vehicle is commonly chosen to achieve a balanced flight by zeroing the lateral air-velocity component v a,2 , which is equivalent to zeroing the sideslip angle β. Provided that an adequate bank angle is created, the tail vertical surface is generally very efficient at maintaining the side-slip angle small without an active yaw control. This corresponds to a positive value of the coefficient C n,β in equation (3.36) which as explained in section 3.3 corresponds to a restoring yawing torque created passively due to the presence of the vertical tail.

This explains why active yaw control via the use of a rotating rudder surface is of secondary importance for most common airplanes. Following the computation of the desired angular velocity ω * given by 6.6, it thus essentially suffices to create, via elevator and ailerons actions, pitch and roll torques that asymptotically stabilize ω 1 -ω * 1 and ω 2 -ω * 2 at zero. This can be achieved without creating a yaw torque with the rudder, and makes k track k. The vertical tail will then make ı and  converge to ī and  respectively, via the generated passive torques.

Case of a flying-wing aircraft

A flying-wing is an unconventional aircraft configuration that lacks a tail and has no fuselage, with most of the payload and equipment being housed inside the main wing (see figure 6.2). It is sometimes regarded as a practical concept due its excellent aerodynamic efficiency, and its light weight. However, because this airframe lacks conventional fixed and control surfaces, it suffers from the absence of inherent passive stability and it is in general difficult to control. Some precautions should be taken at the design level. For instance, in order to ensure longitudinal static stability (as in the case of an ordinary airplane), the equipment should be adequately distributed within the limited space, so that the center of gravity lies ahead of the wing aerodynamic center. It is clear that directional stability would be difficult to ensure due to the absence of a vertical tail, and among many of the proposed solutions, it is considered that a good flying-wing design includes swept-back wings with twisted tips [START_REF] Tonti | Development of a flight dynamics model of a flying wing configuration[END_REF]. Finally, a carefully designed flying-wing would still be marginally stable, and its manual control requires a lot of concentration, therefore it is common to equip such vehicles with onboard controllers that provide stability augmentation.

The actuation configuration of this vehicle consists of two "elevons", which are control surfaces at the trailing edges of both sides of the wing. These elevons replace the ailerons and the elevator. If they deflect together in the same direction, they create a pitching moment, and if they deflect differentially, they create a rolling moment.

The right elevon's deflection that we denote here by δ ev,r is considered to be positive when it deflects upwards, while the left elevon's deflection that we denote by δ ev,l is considered to be positive when it deflects downwards.

Similarly to equation (6. [START_REF] Hua | Time sub-optimal nonlinear PI and PID controllers applied to Longitudinal Headway Car Control[END_REF]) and (6.20), we can relate the active rolling torque Γ 1 and pitching torque Γ 2 to the elevons' deflections according to relations of the following form,

Γ 1 Γ 2 |v a | 2 Ā δ ev,r δ ev,l (6.23)
with,

Ā = η a b C l,δ C l,δ c b C m,δ -c b C m,δ (6.24)
The attitude control design that we presented in this chapter can thus be applied to the flying-wing aircraft, by inverting equation (6.23) for the computation of the desired elevon's deflections (similarly to equation (6.21)). Clearly, no active yawing control torque is intentionally produced with this strategy, and the situation corresponds to the two-axis (pitch-roll) control that was previously discussed.

Trajectory Tracking

T rajectory tracking control is designed in this chapter by extending the methodology of thrust vectoring that is often chosen for the control of low-velocity small scale rotor vehicles with reduced lift surfaces such as the quadrotor. The main challenge is to take into account the dynamics of aerodynamic forces that significantly complicate the design of flight control systems for fixed-wing vehicles. To this aim, the control solution is designed and analyzed on the basis of the model of aerodynamic forces previously described in section 5.3.

We begin by presenting the position control problem. An equilibria analysis follows, where we define a set of admissible reference trajectories. Then, we present the design of the position controller that asymptotically stabilizes the reference position. Finally we show simulation results involving challenging reference trajectories.

Problem statement

Let p r (t) denote a (three times differentiable) reference trajectory in E 3 , with bounded time-derivatives at all orders. In particular v r (t) and a r (t) denote the reference velocity and reference acceleration respectively. The control objective is to stabilize the position tracking error p = p-p r . The control model presented in section 5.2 can now be written as follows,

         d p dt = ṽ dṽ dt = g + Fa m + T m ı -a r d dt (ı, , k) = ω × (ı, , k) (7.1)

Equilibrium analysis

A necessary condition for the asymptotic stabilization of the reference position p r (t), is the existence of an equilibrium of the state on the trajectory. In the case of the proposed hierarchical control architecture, the existence of the equilibrium can be proved by finding locally unique time-functions of the thrust1 T (t) and of the body-fixed frame of the vehicle B r = {G; ı r (t),  r (t), k r (t)} that satisfy the second equation in (7.1) with dṽ dt ≡ 0. This corresponds to the following equilibrium equations,

T ı r = F r (7.2) F r = -mg -F a + ma r (7.3)
In order to give a clear view on the methodology for extending the thrust vectoring control, we start by considering the classical case of vehicles flying in a near hovering mode such as quadrotors, then we show how to extend the analysis for vehicles whose shape can be approximated by a sphere, and finally we solve the equilibrium equations for an airplane whose aerodynamic forces are assumed to take the form proposed in 5.3.

Case of a quadrotor

The control of a quadrotor is classically addressed by neglecting external aerodynamic forces 2 and taking F a = 0. With this assumption, equation (7.3) becomes,

F r = -mg + ma r (7.4)
Now one can apply equation (7.2). The underlying geometrical interpretation is that the thrust vector should be chosen equal to the "apparent external force" F r , which is equivalent to saying that at equilibrium the thrust direction should be parallel to F r , and the magnitude of the thrust counterbalances the external force intensity.

As long as F r is different from zero, there exist only two solutions for the equilibrium equation:

T = |F r | ı r = F r |F r | (7.5) T = -|F r | ı r = - F r |F r | (7.6)
However in practice, the thrust is constrained to be positive and the solution in (7.5) is chosen.

Remark 7.1. With the previous choice of ı r and the condition F r = 0, the equilibrium equations are satisfied for any choice3 of the unitary vector  r which can be used as a free degree of freedom for secondary objectives. And finally k r can be computed according to k r = ı r ×  r to complete the orthonormal base associated with the body frame B r .

Remark 7.2. If the reference trajectory corresponds to a uniform motion with a constant velocity, i.e. a r (t) ≡ 0, one gets for the thrust direction ı r = -k 0 . Which means that the plane of the quadrotor remains horizontal. This is surely not realistic, and is the consequence of neglecting the aerodynamic forces. A more rigorous analysis should include a model for F a , a simple case is to assimilate the vehicle to a spherical body as shown in the next section.

Case of a spherical body

In the case of a spherical body, the aerodynamic force is reduced to its drag component.

According to equation (5.9), we have F a = -η a |v a |c 0 v a , which gives,

F r = -mg + ma r + η a |v a |c 0 v a (7.7)
The air-velocity in this case is equal to

v a = v r -v w .
Notice that F r is independent of the orientation of the vehicle, and as long as F r = 0, we have the same solution for the equilibrium with a constraint positive thrust:

T = |F r |,
and ı r = Fr |Fr| . The choice of  r is still a free degree of freedom that can be controlled independently to achieve a complementary objective.

Case of airplanes

In the case of airplanes, the aerodynamic force F a is no longer reduced to its drag component, and the existence of lifting forces makes F a in general dependent on the orientation of the vehicle. Hence, the solution for ı r is no longer systematic. However, using the model proposed in section 5.3, the aerodynamic force can be written as

F a = -η a |v a |( c0 v a -2c 1 v a,1 ı-( c0 -c0 )v a,2 ).
Note that the first term -η a |v a | c0 v a is independent of the orientation of the vehicle, while the second term +η a |v a |2c 1 v a,1 ı is along the thrust direction and can be lumped with the thrust. More precisely, the equilibrium equation dṽ dt ≡ 0 can now be written as following,

T ı r = Fr -η a |v a |( c0 -c0 )v a,2  r (7.8) T = T + η a |v a |2c 1 v a,1 (7.9
)

Fr = -mg + η a |v a | c0 v a + ma r (7.10)
An additional term (-η a |v a |( c0 -c0 )v a,2  r ) appears now in the equilibrium equation, and the solution for ı r is not independent of the choice of  r anymore, but the entire equilibrium attitude should be chosen simultaneously. Many solutions to this equation might exist, however we choose the one that corresponds to a balanced flight i.e. such that v a,2 = 0.

Using equations (7.8) and (7.9), and setting v a,2 = 0, the corresponding solutions for ı r and for a positive thrust T are the following:

T = | Fr | -η a |v a |2c 1 v a,1 (7.11) ı r = Fr | Fr | (7.12)
As for the unit vector  r , it has to be orthogonal to both ı r and v a (so that v a,2 = 0 is satisfied). Therefore we have  r = ± va×ır |va×ır| . These two possible orientations correspond to flying either cockpit/up or cockpit/down with the aircraft nose facing the incoming air. The common situation is the cockpit/up situation and corresponds to choosing:

 r = v a × ı r |v a × ı r | (7.13)
Finally, the third unit vector k r is just the cross-product of the other two unit vectors.

k r = ı r ×  r (7.14)
These solutions for the equilibrium exist provided that Fr = 0 and v a × ı r = 0 along the trajectory, this leads us to define a set of admissible trajectories as follows, With the above conditions satisfied, | Fr | and |v a × ı r | never cross zero, ensuring the existence of an equilibrium orientation (ı r (t),  r (t), k r (t)) along the reference trajectory as given by equations (7.12),(7.13) and (7.14). It is not difficult to verify that this set of trajectories is much larger than the set of classically defined trim trajectories for which the aircraft translational and angular velocities expressed in the body frame are constant.

Control design and convergence analysis

Let us now focus on control design. We denote the acceleration tracking error as ã = aa r = dṽ dt . Using the expression (5.10) of F a in (7.1) then yields,

mã = (mg -η a |v a |c 0 v a -ma r ) + (T + η a |v a |2c 1 v a,1 )ı + η a |v a |( c0 -c0 )v a,2  (7.16)
Let I p denote a saturated integral of the position tracking error p, and ξ( p, ṽ, I p) denote a bounded PID-like control law that asymptotically (and locally exponentially) stabilizes ( p, ṽ, I p) = (0, 0, 0) for the linear control system ã = ξ, and such that ξ(0, 0, 0) = 0. In view of (7.16),

ã = ξ( p, ṽ, I p) - F m + T m ı + η a m |v a |( c0 -c0 )v a,2  (7.17) F = -mg + η a |v a |c 0 v a + ma r + mξ (7.18) T = T + η a |v a |2c 1 v a,1 (7.19)
We would like the terms Set T = F • ı so that in view of the relation in (7.19), the desired thrust is calculated according to

-F m + T m ı + ηa m |v a |( c0 -c0 )v a,2  in equation (7.
T = F • ı -η a |v a |2c 1 v a,1 (7.21)
Let us assume that | F | is always larger than some positive number, and set the following desired body axis orientation,

ī = F | F | (7.22)  = v a × ī |v a × ī| (7.23) k = ī ×  (7.24)
Assuming that the desired thrust is instantaneously applied, equation (7.17) may also be written as:

ã = ξ( p, ṽ, I p) + | F | m (ı × (ı × ī)) + η a m |v a |( c0 -c0 )(v a • ( -)) (7.25)
In order to avoid useless theoretical complications, we assume from now on that the thrust T applied to the aircraft is bounded. This in turn implies, by virtue of energy dissipation in the air, that |v| is itself bounded, and if the wind speed v w is bounded, then v a is also bounded. This also implies that

| F | = | -mg + η a |v a |c 0 v a + ma r | is bounded.
Therefore it suffices to work out an angular velocity control ω that makes |ı -ī| and | -| converge to zero to ensure the convergence of ( p, ṽ, I p) to (0, 0, 0). This is the core of the control strategy, it implies in particular that the body frame B = {G; ı, , k} converges to the frame B = {G; ī, , k}. This latter problem is well posed because F and v a and, subsequently, the frame B, do not depend on (are not functions of) the airplane orientation. This point is important to properly justify the proposed control design.

As an application of proposition 6.

ω = ω + k ω (t)((ı × ī) + ( × ) + (k × k)) (7.26)
with k w (t) > > 0 renders the equilibrium B = B exponentially stable.

For the sake of simplification it is thereafter assumed that there is no wind so that v w = 0 and v a = v.

Theorem 7.1. Assuming zero wind velocity, if p r (t) is an admissible trajectory, then the control law (7.21)-(7.26) locally exponentially stabilizes ( p, ṽ, I p, B) = (0, 0, 0, Br ).

Sketch of proof.

The angular velocity control law in 7.26 makes o(t) in relation (7.20) converge asymptotically (locally exponentially) to zero a , provided that the thrust control T is chosen according to (7.21). Using the assumption that ξ( p, ṽ, I p) asymptotically (and locally exponentially) stabilizes ( p, ṽ, I p) = (0, 0, 0) for the linear control system ã = ξ, one deduces that under the same conditions stated in 7.1, the control laws (7.21)-(7.26) locally exponentially stabilizes ( p, ṽ, I p, B) = (0, 0, 0, B). In the case where p r (t) is an admissible trajectory, these conditions are, in view of the definition of such a trajectory, satisfied when ξ( p, ṽ, I p) = 0. Moreover B then coincides with the frame Br = {G; ı r ,  r , k r }.

a A more formal proof requires to use an angular velocity control ω that stabilizes | F |(|ı-ī|+|-|) at zero. For more details see [START_REF] Kai | A nonlinear approach to the control of a disc-shaped aircraft[END_REF] and the proof of proposition 4 in [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF]. is always larger than some positive number.

Examples of bounded feedback terms

Saturated PD controller

The feedback part ξ of the controller can be designed as a nonlinear saturated PD controller, as follows:

ξ( p, ṽ) = -k p s at ∆p ( p) -k d s at ∆v (ṽ) (7.27) = -k p α ∆p (| p|)( p) -k d α ∆v (|ṽ|)(ṽ) (7.28)
Where ∆ p and ∆ v are positive numbers, and k p , and k d are positive gains. Using the Lyapunov function V ( p, ṽ) = 1 2 |ṽ| 2 + k p ´| p| 0 α ∆p (s)sds ≥ 0, it can be shown that its time derivative is equal to V = -k d sat ∆v (|ṽ|)|ṽ| 2 ≤ 0. Applying LaSalle's invariance principle shows that the origin of this system is globally asymptotically stable and locally exponentially stable.

Saturated PID controller

To compensate for modeling errors and slowly time-varying perturbation terms, an integral term should be added. This term should be however bounded in order to avoid singularities in specifying the desired attitude. Among different possibilities, the bounded conditional integrator I p defined in [START_REF] Hua | Time sub-optimal nonlinear PI and PID controllers applied to Longitudinal Headway Car Control[END_REF] is chosen. It is the solution of the following differential equation (with I p = İ p = 0),

Ï p = -k dI İ p + s at ∆ Ï 2 (k pI (-I p + s at ∆ I (I p + p k pI ))) (7.29)
Where k dI , ∆ Ï , k pI and ∆ I are positive numbers. It is shown in [START_REF] Hua | Time sub-optimal nonlinear PI and PID controllers applied to Longitudinal Headway Car Control[END_REF] that |I p|, | İ p| and | Ï p| are bounded respectively by

∆ I + ∆ Ï 2k 2 dI , ∆ Ï 2k dI and ∆ Ï .
An expression of a stabilizing nonlinear PID controller using the previous integrator is the following: For which in analogy to the PD controller defined previously, we can deduce the convergence of ( p, ṗ) to zero which also ensures the convergence of ( p, ṽ) to zero even in the presence of a constant perturbation term (for details see [START_REF] Hua | Time sub-optimal nonlinear PI and PID controllers applied to Longitudinal Headway Car Control[END_REF]).

ξ( p, ṽ, I p) = -k p s at ∆p ( p + k I I p) -k d s at ∆v (ṽ + k I İ p) -k I Ï p (7.

Simulations

The object of this section is to illustrate the tracking performance of the control (T, ω)

given by (7.21) and (7.26) applied to an aerial vehicle weighting 3Kg and modeled similarly to the control model presented in 5.3, with η a = 0.55, c 0 = 0.01, c 1 = 1.5, c0 = 3.

In order to test the robustness of the control against modeling errors, we choose model parameters for the control computation that are slightly different than those of the simulated model: m = 2.7, η a = 0.48, c 0 = 0.01 and c 1 = 1.5. The thrust calculated with these values is applied to the aircraft with the multiplication factor of 0.8.

Furthermore the applied thrust is constrained to be nonnegative.

The reference trajectory used for this simulation consists first of a straight accelerated trajectory at 2.4m/s 2 followed by a half-circle left-turning maneuver with a speed of 20m/s and a radius of 80m, then a right-turning maneuver with the same speed and radius, after which the trajectory is decelerated until the speed reaches 15m/s. Finally,

The trajectory consists of a vertical circular loop with a speed of 15m/s and a radius of 15m then a deceleration along a straight line.

The bounded feedback term ξ is calculated according to the following:

ξ( p, ṽ, I p) = -5 s at 20 ( p + I p) -5 s at 20 (ṽ + İ p) -Ï p Ï p = -0.5 İ p + s at 1.5 2 (20(-I p + s at 4 (I p + p 20
)

))

The gain k ω (t) involved in the expression (7.26) of ω is chosen constant and equal to 10.

The initial distance between the aircraft and the reference trajectory is equal to 52m. Hardware-in-the-loop simulation results involving a scale-model aircraft and challenging reference paths, with large initial tracking errors, and air-velocity measurements approximations are reported in section 8.7.

Problem Statement

Consider a three-times differentiable curve C in 3D-space parametrized by its curvilinear abscissa s, and let Q(s) represent the point on the curve closest to the center of mass G of the airplane. Depending on the curve, this point can be always unique (as in the case of a straight line) or only locally unique, depending on the position of the airplane w.r.t the curve (as, for instance, in the case of a circle for which uniqueness is granted provided that G does not belong to the circle axis passing through the origin and perpendicular to the circle's plane). Let q denote the position of the point Q, and let v * ∈ R + denote the desired magnitude of the airplane's velocity. Define the position error vector p = pq.

A way to achieve the path following objectives consists in regulating e v = |v| -v * and p at zero. In order to structure the control design according to the hierarchical architecture explained in section 5.1, the control model of section 5.2 is divided to the following subsystems, where each is associated with a specific control task:

• Kinematical guidance: This is the outer loop that involves body kinematics equations only. It corresponds to the following equation,

ṗ = v - q = |v|h -q (8.1)
In this equation, we consider that the heading vector h = v |v| (provided that v is different than zero) plays the role of an intermediate control variable for the guidance loop.

• Velocity control: This control block monitors the speed |v| via thrust adaptation. It corresponds to an equation of the form,

ėv = f (T ) (8.2)
that will be detailed later. The corresponding control variable is thus the thrust T . We will also discuss the possibility to regulate the airspeed |v a | instead of the inertial speed.

• Heading stabilization: This control block is related to the monitoring of the heading h and to achieving a balanced flight (zeroing v a,2 ) via the aircraft attitude as an intermediate control input. The corresponding equation is the translational dynamics equation, written here using the proposed aerodynamic model (5.10) as follows,

a = ḡ + T m ı + v a,2 Ō(v a ) (8.3) with
, Ō(v a ) = η a |v a | c0 -c0 m , ḡ = g - η a c0 m |v a |v a (8.4)
and,

T = T + 2η a c 1 v a,1 |v a | (8.5)
• Attitude control: This is the inner loop subsystem related to the stabilization of the attitude using the angular velocity ω. This control block involves the following equation,

d dt (ı, , k) = ω × (ı, , k) (8.6)
The angular velocity is then stabilized using the control surface angles δ according to the control design in chapter 6.

In accordance with this control architecture, the decomposition of v into the product of |v| and h is all the more justified that convergence to the desired path can be per- formed at various speeds with the same heading policy and, vice versa, that changing the heading policy does not imply modifying the speed.

Velocity Control

The proposed control design involves the interconnection of these different stages as shown in figure 8.2. Each of these blocks will be detailed next.

Velocity Control

Using equations (

) and the definition of e v = |v| -v * , one finds that,

d dt e v = ḡ • h + (ı • h) T m -v * + v a,2 Ō(v a ) • h (8.7) with, Ō(v a ) = η a |v a | c0 -c0 m .
This relation in turn suggests to set,

T = m -ḡ • h + v * -k T,1 e v -k T,2 α e I ev /ı • h (8.8)
with k T,1 and k T,2 denoting positive gains and I ev a bounded integral of the velocity error e v defined as,

İev = k T,3 -I ev + s at ∆e v (I ev + e v /k T,3 ) (8.9)
with ∆ ev denoting a positive number, and k T,3 a typically large positive number. The scalar function α e is defined by,

α e (e v , I ev ) = α ∆e v (|I ev + e v /k T,3 |) (∈ (0, 1]) (8.10)
The integrator I ev in 8.9 is ultimately bounded by ∆ ev , and | İev | is ultimately bounded by 2k T,3 ∆ ev .

Proof. From equation (8.9), we can deduce that,

0.5 d dt I 2 ev = -k T,3 (I 2 ev -I ev s at ∆e v (I ev + e v /k T,3 )) ≤ -k T,3 |I ev |(|I ev | -∆ ev ) (8.11)
Therefore |I ev | decreases when |I ev | ≥ ∆ ev , and ∃t 0 :

∀t ≥ t 0 , |I ev | ≤ ∆ ev
From equation (8.9), and the result on the ultimate boundedness of |I ev |, we have, ensures the exponential convergence of e v to zero. When v a,2 ≡ 0, the origin of the feedback controlled system (8.12) is globally asymptotically and locally exponentially stable.

∀t ≥ t 0 | İev | ≤ k T,3 |I ev | + k T,3 ∆ ev ≤ k T,3 ∆ ev + k T,3 ∆ ev = 2k T,3 ∆ ev
Proof. Consider the candidate Lyapunov function L = 0.5e 2 v + 0.5k T,2 I 2 ev . Along any solution to the system one verifies that,

L = -k T,1 e 2 v -k T,2 K T,3 (1 -α e )I 2 ev + v a,2 e v Ō(v a ) • h (8.13) • If ∀t v a,2 = 0 then, L = -k T,1 e 2 v -k T,2 K T,3 (1 -α e )I 2 ev ≤ 0 (8.14)
Note that we are stating the result assuming that the perturbation vanishes, it is possible however to prove that with the chosen control term, and with the presence of a constant bounded and non-vanishing perturbation, the system converges to an equilibrium (ev, Ie v ) = (0, c) with c = 0, a constant that depends on the magnitude of the perturbation. For details see [START_REF] Hua | Time sub-optimal nonlinear PI and PID controllers applied to Longitudinal Headway Car Control[END_REF].

From which we can deduce the global asymptotic stability of (e v , I ev ) = (0, 0).

When e v and I ev are close to zero, one can verify that the variations of e v and I ev satisfy the equations of a second order linear system. Indeed one gets ëv + k T,1 ėv + k T,2 e v = 0 from which we can deduce that the equilibrium is locally exponentially stable.

• If ∃t : v a,2 = 0, and assuming that |v a,2 | converges to zero exponentially, the ultimate exponential convergence follows from the fact that an additive, exponentially vanishing, perturbation applied to a system whose origin is exponentially stable does not prevent the solutions to this system from converging to zero exponentially.

How to control the aircraft attitude so as to make v a,2 converge to zero, and thus achieve a so-called "balanced" flight with zero sideslip angle, is addressed later on via the control of the aircraft attitude. In order to avoid theoretical complications of little practical importance, boundedness of the aircraft velocity is assumed. This assumption can also be justified because of thrust physical limitation and energy dissipating drag forces. Note also that the integral action is not only useful in practice to compensate for imprecisely known pre-compensation terms, but also to compensate for the imperfect knowledge of the physically applied thrust.

Kinematical Guidance

Kinematical relations

Consider again the definition of the point Q(s) on the reference path. We define at the point Q(s) on this path, an associated parallel transport frame F = {Q; u, ū, ū} with u the vector tangent to the path C at the point Q (see figure 8.1). An advantage of a parallel transport frame over the more conventional Frénet frame is that it is continuously welldefined at points where the path curvature vanishes. It does not suffer from ambiguity and sudden orientation changes when the curves straightens out [START_REF] Bishop | There is more than one way to frame a curve[END_REF] [START_REF] Andrew | Parallel transport approach to curve framing[END_REF]. Its relative drawback is that it is not defined from the sole curve characteristics (curvature and torsion). More precisely, it is uniquely defined only once the vectors ū and ū are arbitrarily chosen at some point on the curve. The corresponding variational frame equations are,

d ds ū = -γ 1 u ; d ds ū = -γ 2 u ; d ds u = γ 1 ū + γ 2 ū (8.15)
A formal proof on the boundedness of |v| can be found in [START_REF] Kai | A unified approach to fixedwing aircraft path following guidance and control[END_REF].

to u, and therefore, p • u = 0. Differentiating this equality gives:

ṗ • u + p • u = 0 v • u -ṡ + ṡ p • (γ 1 ū + γ 2 ū) = 0
From which we can deduce that,

ṡ = (u • v) 1 -γ 1 y 1 -γ 2 y 2 (8.21)
For the sake of avoiding problems of little practical relevance, we will assume from now on that the chosen path is such that γ 1 and γ 2 are uniformly bounded. The uniqueness of the projection of the aircraft center of mass on the path is then granted provided that 1 -γ 1 y 1 -γ 2 y 2 is larger than some positive constant. In the case of a straight line, for which γ 1 = γ 2 = 0, this condition is thus satisfied independently of the aircraft position. If the condition is satisfied and |v| is bounded, then the time-derivative of u is also bounded.

Guidance

By viewing the aircraft as a point moving in 3D-space with a given speed |v|, the problem is to determine a desired heading direction h * that yields the convergence of this point to the desired path and ensures that the point moves thereafter along the path with the desired direction given by ±u, i.e. such that h * converges to sign vu u with sign vu chosen in advance and equal to either 1 or -1. There are obviously a multitude of solutions to this problem as was discussed in section 4.3. The solution that we propose here consists in setting,

h * = sin(θ h )l + (cos θ h )sign vu u (8.22)
with l denoting some unit vector orthogonal to u and θ h an angle depending on the position error p and converging to zero when | p| tends to zero. For instance, a simple possible choice for l and θ h is,

ȳ = k 1 D s at ∆ h (y)/|v| (8.23a) l = - ȳ1 ū + ȳ2 ū |ȳ| (8.23b) θ h = arctan(|ȳ|/ 1 -|ȳ| 2 ) (8.23c) with D = diag{d 1 , d 2 } a diagonal matrix with d i ∈ (0, 1] (i = 1, 2), ∆ h = µ|v| k 1 max(d 1 ,d 2 ) ,
µ ∈ (0, 1), and k 1 a positive gain. Note that |ȳ| < 1 so that θ h is well defined and belongs to [0, π/2). Note also that sin(θ h ) = |ȳ| ≤ µ so that the term sin(θ h )l entering the expression of h * is always well defined. Moreover, if d 1 = d 2 then sin(θ h ) tends to µ when |y| tends to infinity i.e. when the aircraft is far from the path. Therefore, in this case arcsin(µ) characterizes the angle of incidence of the desired heading direction w.r. Also note that the previous choices of ȳ, l and θ h are equivalent to writing the desired heading vector as follows:

h * = -ȳ 1 ū -ȳ2 ū + 1 -|ȳ| 2 sign vu u (8.24)
The following proposition summarizes the stability and convergence properties associated with this desired heading direction.

Proposition 8.2. Assume that 1 -γ 1 y 1 -γ 2 y 2 and |v| are always larger than some positive number so that the position error p and the aircraft heading vector h are always well defined.

Assume that |v| is bounded and that h = h * +o, with o denoting a "residual" vector such that the integral ´t 0 |o(s)|ds is bounded. Then |y| converges to zero, and h tends to sign vu u. The convergence is ultimately exponential if the convergence of |o| is itself ultimately exponential.

Moreover | ẏ| is ultimately upper bounded by µ|v| and | ẏi | is ultimately upper bounded by

d i max(d 1 ,d 2 ) µ|v|(i = 1, 2).
In the case where o ≡ 0 the equilibrium p = 0 is locally exponentially stable.

Proof. Define ō = |v|Π u o. The time-variation of pF is given by, That is to take into account the interconnection between different subsystems and the fact that although h is controlled to converge to h * , it may not be equal to h * at all times. ω h * of h * . This suggests to use a bounded integral term calculated according to,

ṗF = Π u v = |v|Π u (h * + o) = -|v|(ȳ 1 ū + ȳ2 ū) + ō (8.
ż = ω h * × z + k z -z + s at ∆z (z + h/k z ) ; z(0) = 0 (8.28)
with ∆ z > 0 the chosen upper bound for |z(t)|, and k z denoting a positive number.

The expression (8.27) of ωh is then modified to

ωh = ω h * + k h,1 h + k h,2 α h z (8.29)
with k h,1 denoting a positive gain and the scalar function α h defined by

α h ( h, z) = α ∆z (|z + h/k z |) (∈ (0, 1])
From now on, the arguments of this function are omitted for the sake of legibility.

Proposition 8.3. Assume that 1 -γ 1 y 1 -γ 2 y 2 and |v| are always larger than some positive number so that the position error p, the desired heading vector h * , and the aircraft heading vector h are always well defined. Assume also that ω h = Π h ωh + o, with ωh given by (8.28)-(8.29) and o a "residual" vector such that the integral ´t 0 |o(s)|ds is bounded. We distinguish two cases:

case 1: ∀t : o(t) = 0.
• In this case, the system (8.26), (8.28), (8.29) has two equilibria, namely (h, z) = (h * , 0) and (h, z) = (-h * , 0). The first of these equilibria is locally exponentially stable, whereas the second one is unstable;

• (h, z) converges to the first (desired) equilibrium provided that h(0) = -h * (0).

case 2: ∃t : o(t) = 0.

• If (h, z) does not converge to the unstable asymptotic equilibrium (-h * , 0) then it converges to the desired asymptotic equilibrium (h * , 0).

• If |o| converges ultimately exponentially to zero, then the convergence of (h, z) to (h * , 0)

is also ultimately exponential.

Proof. Forming the time-derivative of the positive function

V 0 = (1 -h • h * ) + 0.5k h,2 |z| 2 yields, V0 = -k h,1 | h| 2 -k h,2 k z (1 -α h )|z| 2 -o • h (8.30) case 1: ∀t : o(t) = 0.
In the absence of wind, this frame coincides with the wind frame defined in section 2.2.3.

In this case V0 ≤ 0, and V0 = 0 if and only if (h, z) = (h * , 0) or (h, z) = (-h * , 0). It is simple to verify that these two points, for which V 0 is stationary, are indeed equilibria of the system (8.26), (8.28). Because V 0 is non-increasing, and α h (0, z) = 1

only if z = 0, the convergence of | h| and |z| to zero follows. This in turn implies that (h, z) converges to one of the system's equilibria. Because

1-h(t)•h * ≤ V 0 (t) ≤ V 0 (0)
is always smaller than two when h(0) = -h * , h cannot converge to -h * , and thus necessarily converges to the desired equilibrium.

Another way of establishing the stability, or instability, properties of the system's equilibria consists in studying linear approximations of the system about these equilibria. One verifies that, in the first order approximation about ( h, z) = (0, 0), the variations of h and z satisfy the following equations:

ḣ = ω h * × h ∓ (k h,1 h + k h,2 z) ż = ω h * × z + h (8.31)
With the sign in the right-hand side of the first equality depending on the chosen equilibrium, i.e. (h, z) = (h * , 0) or (h, z) = (-h * , 0). The minus sign goes with the first equilibrium, and the plus sign with the second one. Consider a frame rotating with the angular velocity ω * h , and let x 1 (resp. x 2 ) denote the two-dimensional vector of Cartesian coordinates (in this frame) of the projection of h (resp. z) onto the plane orthogonal to h * . System (8.31) becomes equivalent to the following linear system,

ẋ1 = ∓(k h,1 x 1 + k h,2 x 2 ) ẋ2 = x 1 (8.32)
The characteristic polynomial associated with the first (resp. second) one is

(λ 2 + k h,1 λ + k h,2 ) 2 = 0 (resp. (λ 2 -k h,1 λ -k h,2 ) 2 = 0).
All poles of the first system have negative real parts, whereas two poles of the second system are real positive. This in turn proves that the equilibrium (h, z) = (h * , 0) is exponentially stable, and that the other equilibrium is (exponentially) unstable.

case 2: ∃t : o(t) = 0.
Let us first establish that ( h, z) converges to (0, 0). The assumed boundedness of ´t 0 |o(s)|ds implies that ´t 0 o(s) • h(s)ds is also bounded. Since V 0 is positive and uniformly bounded from above, |V 0 (t) -V 0 (0 The convergence of ( h, z) to (0, 0) in turn implies that (h, z) converges either to (h * , 0) or to (-h * , 0). Therefore, non-convergence to the unstable point (-h * , 0)

) + ´t 0 o(s) • h(s)ds| = k h,1 ´t 0 | h(s)| 2 ds + k h,2 k z ´t 0 (1 -α 1 (s))|z(s)|
implies convergence to the desired stable point (h * , 0). The ultimate exponential rate of convergence when |o| converges to zero exponentially follows from the fact that an additive, exponentially vanishing, perturbation applied to a system whose origin is exponentially stable does not prevent the solutions to this system from converging to zero exponentially.

Let us now define the desired mobile frame B. Recall the expression of acceleration in equation (2.12), a = |v|h + |v|(ω h × h). This relation suggests to define a "desired" acceleration as follows,

a * = v * h + |v|( ωh × h) (8.33)
with ωh given by (8.28), (8.29). From relation (8.3) we note that

ı = a -ḡ -v a,2 Ō(v a ) |a -ḡ -v a,2 Ō(v a )| (8.34)
The desired acceleration is in turn used to define ī as follows (compare with (8.34))

ī = a * - ḡ |a * -ḡ| (8.35)
The side component v a,2 of the airspeed should be regulated to zero, therefore, for the vector , we set

 = v a × ī |v a × ī| = v a × (a * -g) |v a × (a * -g)| (8.36)
so that v a and  are orthogonal, and the third vector k is just calculated as the cross product of ī and , i.e. k = ī ×  (8.37)

An important property is that, like p, v, v a , g, ω * h and ω h , the unit vectors (ī, , k) so defined do not depend on the aircraft attitude. Therefore, their time-derivatives do not depend on the aircraft angular velocity ω. Let ω ī := ī × i and ω  :=  ×  denote the angular velocities of ī and  respectively. The angular velocity of the frame B is then given by ω = ω ī + (ī.ω )ī = ω  + (.ω ī), and this vector does not depend on ω either.

The problem of stabilizing B = B is thus well-posed. Proposition 8.4. Assume that 1 -γ 1 y 1 -γ 2 y 2 and |v| are always larger than some positive number so that the position error p, the aircraft heading vector h, and the desired heading vector h * are always well defined. Assume also that |a * -ḡ|, |ī × v a |, are always larger than some positive number so that the frame B and its angular velocity ω are well-defined. Then an angular velocity control that almost globally asymptotically (locally exponentially) stabilizes

B = B is ω = ω + k ω (t) (ı × ī) + ( × ) + (k × k) (8.38)
with k ω (t) > > 0.

A proof of this proposition is similar to the proof of proposition 6.1.

Prior to stating overall stability result, the following theorem summarizes convergence properties that can be established from the partial results obtained so far. • h| is, by assumption, always larger than some positive number, (|a -ḡ| -|a * -ḡ|) converges to zero. Now combining this with the fact that Π h ξ converges to zero, one deduces that ( ḣωh × h) also converges to zero. Because ω h := h × ḣ and ω h • h = 0 one then infers that ω h converges to Π h ωh . All convergence rates evoked so far -the rate of convergence of ω h to Π h ωh in particular-are ultimately exponential.

Therefore ω h = Π h ωh +o, with |o| vanishing ultimately exponentially, and thus such that the integral ´t 0 |o(s)|ds is bounded. By application of Proposition 8.3,if (h,z) does not converge to the asymptotically unstable point (-h * , 0), then it converges to (h * , 0) with a rate of convergence also ultimately exponential. In this latter case, by application of Proposition 8.2, the ultimate exponential convergence of the path following error y to zero, and of h to sign vu u, follows.

The conditions pointed out in Theorem 8.1, under which convergence to the desired path is granted, may seem restrictive at first glance; but they are in fact inherent to the control problem at hand. They are also related to the existence of particular trajectories along which the linearized equations of the system are not controllable. Although they are not satisfied in only very specific situations, they nonetheless rule out the possibility of global convergence results. However, it remains possible to state local asymptotic stability results when these conditions are satisfied on the desired path. For instance, in the case of zero wind velocity, and when |v| = v * is constant, one verifies that these conditions are satisfied on the desired path if |u×ı| and u•ı are positive (and larger than a small number) on the path. Since, for a balanced flight, ı = a-ḡ |a-ḡ| with a = v * 2 (γ 1 ū+γ 2 ū) and ḡ = gηac 0 m v * 2 u, these conditions are themselves satisfied if

A1: |g × u -v * 2 (γ 1 ū -γ 2 ū)| > 1 > 0 (8.39)
and

A2: | η a c0 m v * 2 -g • u| > 2 > 0 (8.40)
We can then state the following local exponential stability result. With â denoting an estimate of the acceleration a. For this latter estimate, one may use onboard accelerometers that measure ga or, alternatively, assume that the variations of v in the body frame are slow so that â = ω × v. An even cruder estimation is â = 0. Finally assuming that the sideslip velocity v a,2 is kept small by virtue of the lateral energy dissipation and the passive lateral stability of the plane (via its vertical stabilizer and wing dihedral) and re-enforced by the attitude controller, one can assume that va,2 = 0. The resulting estimate of v a , used for the implementation in this work is va = v a,1 ı + va,3 k with â = 0.

Airspeed control

Instead of stabilizing the inertial velocity |v| one may wish to monitor the air velocity in the direction of ı, i.e. the component v a,1 = v a .ı, which can be measured, for instance, with a pitot tube. This is a more common situation and a more secure choice for flight tests, since maintaining the airspeed at nominal values ensures the availability of enough lifting forces and decreases the chances of stall. This choice is even more crucial for small UAVs for which the difference between the inertial velocity and airspeed can easily become significant in the presence of wind.

Define now the speed error as e v = v a,1 -v * . Using (8.3) the time-derivative of e v Via the balanced flight policy.

For simulations and for flight tests.

satisfies the relation,

d dt e v = d dt ((v -v w ) • ı) -v * = (a -vw ) • ı + ω • (ı × v a ) -v * = (g -vw ) • ı + ω • (ı × v a ) -ηac 0 m |v a |v a,1 -v * + T m + v a,2 Ō(v a ) • ı
Assuming that v a,2 converges exponentially to zero (balanced flight), exponential convergence of v a,1 -v * to zero is then obtained by setting

T = T * -mk T,1 e v -mα e k T,2 I ev (8.42)
with In the case of airplanes equipped with electrical motors.

T * = -m((g -vw ) • ı + ω • (ı × v a ) -ηac 0 m |v a |v a,1 -v * ), k T,

Thrust bounds and attack angle monitoring

We have so far assumed that the aircraft could produce the desired thrust T calculated according to (8.5)- (8.8). In practice this desired value of T may leave the physical thrust interval [T min , T max ]. When this happens at least one of the control objectives -i.e. convergence of the aircraft's heading direction h to the desired one h * , or stabilization of |v| at the desired speed v -cannot momentarily be achieved with the available thrust.

For instance, descending from a high altitude to a horizontal path with the convergence dynamics specified by h * may require a negative thrust (to slow down the aircraft) that a common aircraft cannot produce. Similarly, climbing with these convergence dynamics and velocity may require a thrust that exceeds T max .

To avoid this situation, one of the possibilities consists in reducing the rate of convergence to the desired path in order to have the control action focused on the stabilization of the aircraft velocity. This can be done, for instance, by choosing the parameter µ in the expression of the saturation ∆ h = µ|v| k 1 max(d 1 ,d 2 ) small enough. In the case of a horizontal path one may alternatively set the parameter d 2 small enough so as to impose a small rate of climb or descent, during the transient phase of convergence to the path, that is compatible with the thrust limitations.

One may also temporarily accept an increase of the aircraft's velocity beyond the desired speed v * during a descending transient phase when a negative thrust (needed to slow down the aircraft) is calculated and cannot be produced (i.e. when T min = 0). In this case only the objective of convergence of the heading direction h to the desired one h * is maintained. In this case the desired speed v * (t) is supposed to be equal to the actual resulting speed of the aircraft |v(t)|, this in turn yields, just as in the case of air-velocity control, to setting v * = d dt |v| in relation (8.33).

A different issue concerns positive thrust limitations. Indeed, while unlimited thrust power theoretically allows one to control an aircraft at any speed and attack angle, a finite value T max automatically limits the aircraft's speed. Moreover, when the maximal thrust is significantly smaller than the aircraft's weight it critically matters to keep the attack angle small and under the stall value. Without this safety feature the direction ī that is calculated (or imposed by a pilot) without taking this limitation into consideration may yield a large angle of attack and a loss of lift leading to a continuing descent even at full thrust. To automatically integrate this safety feature in the proposed control design let us rename the unit vectors {ī, , k} calculated without taking thrust limitations into account as {ī * ,  * , k * }. Let also α max denote the desired upper-bound for the attack angle.

Whenever the predicted attack angle α * = arcsin( va |va| . k * ) is larger than α max we propose to determine new desired attitude directions for the aircraft such that the associated attack angle is equal to α max . More precisely, we propose to set

ī = rot(α max  * ) v a |v a | ,  =  * , k = ī ×  (8.43)
with rot(α max  * ) va |va| denoting the vector va |va| rotated by the angle α max about the unit vector  * . In this case the pre-compensation velocity ω is calculated with ω ī = ī × i and, in view of (8.43)

i = α max ω  * × ī + rot(α max  * )Π va |va| va |v a |

Transition between reference paths

Suppose that the vehicle is following a path (C 1 ), and that once it reaches (or gets close to) a point W 1 ∈ (C 1 ), it has to follow another curve (C 2 ).

Many switching policies can be employed. In this work we decide that the switching occurs as soon as the vehicle enters a sphere centered on W 1 with a specified "acceptance radius" that we denote by r ac . With this policy, the guidance algorithm switches to following (C 2 ) once the distance of the center of mass of the vehicle to W 1 is less or equal to r ac : | --→ GW 1 | ≤ r ac . However, this strategy has a drawback. If the vehicle fails to enter the sphere, it will continue following (C 1 ) beyond the point W 1 . Therefore another condition should be added to the switching policy, such as a half-plane switching criteria (see [5, chapter 11]). These particular situations are not addressed in this thesis, and may be the subject of future developments.

Controller hardware and software

A "Pixhawk" controller, equipped with a 168 MHz ARM CPU and 256 KB of RAM, is used as the autopilot hardware. Our code implementation is based on the open-source PX4 flight stack that runs on top of a NuttX real-time operating system and uses the PX4 middleware [START_REF] Lorenz | Px4: A node-based multithreaded open source robotics framework for deeply embedded platforms[END_REF]. This software architecture runs different threads with a modular approach and handles inter-process communications, allowing for the development and use of off-the-shelf control code. We took advantage of this possibility to replace the pre-existing position and attitude control modules by our own libraries in order to implement the proposed control algorithms. However we kept the other pre-existing functionalities and, in particular, the extended Kalman Filter (EKF) state observer. The software logs all the flight data on a SD card, allowing the flight information to be analyzed after the flight.

Ground control station

We use the "Qgroundcontrol" software (qgroundcontrol.com) as a control station installed on a base computer, to design missions and change parameters online during flights. In our simulation setup, Qgroundcontrol establishes a communication with the "PIXHAWK" controller using the "MAVLINK" protocol. In parallel, Qgroundcontrol establishes a UDP link with the X-plane simulator to send and receive data. This dual connection allows Qgroundcontrol to perform hardware-in-the-loop simulation by allowing an indirect communication between the controller and the simulator. Simulated position and attitude of the aircraft are used to create virtual sensory measurements (GPS, IMU, barometer, pitot) that are purposefully corrupted with artificial noise and produced (sent) at a realistic rate. These simulated sensor's measurements are handled effectively by the state estimator whose outputs are very satisfactory when compared with the real state of the vehicle obtained in the simulator. Since sate estimation is not the objective of this work, we will not detail the algorithms regarding sensor data generation or state estimation, interested readers can refer to the source codes of Qgroundcontrol and PX4. The controller calculates PWM control values and sends them to the ground station, which in turn sends them back to the simulator that uses them as setpoints for the generation of thrust and for control surfaces deflections.

Which is in this case a simulation.

Simulation results

For these simulations we consider, as in a realistic case, that a priori unknown wind may be present, and we choose to stabilize the airspeed v a,1 at 12m/s. Besides the three-axis controller, we have also tested the two-axis (pitch-roll) adaptation proposed in section 6.4.

The aerodynamic coefficients used for the control calculations are: η a c 0 = 0.006, and η a c 1 = 0.5. The choice of these coefficients need not be precise, indeed the integral action included in the control design is supposed to compensate for these modeling errors. However in order to use coefficients with the correct order of magnitude, one can refer to measurements data for NACA profiles as in section 3.2.1, then apply the finite-wing correction according to section 3.2.2, where the coefficients are adapted to the characteristics of the wing in question (surface and span).

The control parameters that are used are:

• Guidance: k 1 = 1, µ = 0.5, d 1 = 1, d 2 = 0.5 • Airspeed control: k T,1 = 1.8, k T,2 = k T,1 /2, k T,3 = 10, ∆ ev = 10 • Heading stabilization: k h,1 = 1.4, k h,2 = 0.49, ∆ z = 0.6, k z = 10 • Attitude control: k ω = 7.0
During implementation, the saturation functions s at ∆ (x) = α ∆ (|x|)x entering the control laws were replaced by the classical vector-valued saturation functions, i.e. with α ∆ (|x|) = max(1, ∆ |x| ). As explained in section 6.3, control surfaces angles are calculated

according to δ * = -K δ |va| 2 (w -w * ) = -1 |va| 2 diag([90, 120, 90])(w -w * ).
For the thrust calculation we used k T = 18.

Simulation 1: a custom 3D path

The chosen closed reference curve (see Fig. 8.5) consists, for the first part, of a horizontal segment connected to a horizontal half-circle of radius equal to 40m, followed by another horizontal segment. The second part of the reference path is similar to the first one except that it is inclined with an angle of 15 • w.r.t. the horizontal plane. Just recall here that an inclined circular path does not qualify as a trim trajectory.

The first set of simulations are carried out with no wind so that one can appreciate the controllers performance in this "ideal" case.

Simulations results are reported in Fig. 8.5 and Fig. 8.6. They show that the controller allows the aircraft to approach the desired path and then follow it closely. Despite the approximations involved in the estimation of the air-velocity v a , the imperfections of the model used for the control design, and the noisy state estimates produced by the embedded EKF observers, the controller nominally maintains the magnitude of the position error well under the accepted norm of a wingspan. Larger errors only occur when the aircraft has to switch from one piece of the curve to the next one. The switching is done according to an acceptance radius as explained in section 8.6.5.

The chosen acceptance radius is chosen differently for every point according to the situation, it is most of the time equal or close to 5m. The same simulation is repeated using the two-axis pitch-roll attitude autopilot. This controller shows similar behavior to the three-axis autopilot (see figures 8.8 and 8.9). The aerodynamic coefficients used for the control calculations are: η a c 0 = 0.003, and η a c 1 = 0.08, and the control parameters are: The airspeed v a,1 is chosen to be stabilized at 15m/s. In order to test more aggressive maneuvers, we choose as a reference path an inclined circle with a 25 • inclination w.r.t.

• Guidance: k 1 = 1, µ = 0.5, d 1 = 1, d 2 = 0.4 • Airspeed control: k T,1 = 0.9, k T,2 = k T,1 /2, k T,3 = 10, ∆ ev = 10 • Heading stabilization: k h,1 = 1, k h,2 = 0.25, ∆ z = 0.6, k z = 10
the horizontal plane, and a radius of 35m. After some time the reference path is suddenly switched to a vertical circular path with a radius of 35m as well, and situated far from the first one.

The flying-wing is initially far from the inclined circle, and it is seen in figure 8.17

that it approaches the path and manages to follow it with an acceptable error of magnitude less than 2.5m. The vehicle then quits the first path and approaches the vertical circle, then follows it in a looping maneuver while keeping the position error less than 4m. The airspeed v a,1 is also regulated efficiently at 15m/s, it increases above the setpoint during descending phases where the theoretical thrust reaches its lower saturation at zero. We first present the airplane properties and the hardware architecture. Then we describe the experimental setup and the progress of the flight tests. Finally we show results involving challenging reference paths, with large initial position errors.

Demonstration vehicle

The vehicle used for the tests is an E-flite "TIMBER" airplane shown in figure 9.1. This airplane is designed to fly at low airspeeds. It is equipped with wheels as well as a steerable tail wheel which gives it taxiying capabilities and also short takeoff and landing (STOL) capabilities.

Some of its basic features are the following,

• Wingspan 1.5m

• Wing area 0.36m 2

• Fuselage length 1m

Avionics and actuators

The autopilot hardware that we used is a "PIXRACER" controller board, another controller of the pixhawk family. It is equipped with a 180M hz ARM CPU, and 256KB

of SRAM. This controller is dedicated to mobile robotics and is equipped with embedded Inertial Measurement Units (IMU) which provide measurements redundancy. The models of these sensors modules are the following:

• Invensense ICM-20608 Accel / Gyro

• Invensense MPU9250 Accel / Gyro / Mag

• Honeywell HMC5983 magnetometer

• MS5611 barometer

This board also has a slot for a microSD card that is used for logging data. It has an interface to a safety switch and buzzer. An I2C bus allows to connect to peripheral ICs that uses the I 2 C serial communication protocol, like for example differential pressure sensors and external magnetometer modules. Radio modems can be connected to telemetry ports to exchange live data with a ground base computer. PPM-Sum or S.Bus receivers can be connected to the board to achieve manual radio control. And finally PWM output ports are connected to servomotors and ESCs, the pulse modulated signal they generate has a period of 2000µs and an amplitude above 2 Volts, with the minimum value corresponding to a pulse width of about 900µs. The software code is the same as the one tested in HITL simulations, except that it is now compiled for the architecture of the Pixracer.

Experimental setup

Prior to going to the field, the sensors were calibrated using routines provided by the ground station (qgroundcontrol). The controller was also loaded with verified initial

The Pixracer belongs to the FMUv4 generation, while the early version of pixhawk boards belongs to the FMUv2 generation. As in the case of HITL simulations, reference paths were prepared using two different methods:

• The classical waypoint tools provided by the ground control station. These are a set of points in space for which the position is determined by longitude, latitude and altitude values. Standard waypoints are connected with straight segments and constitute straight line paths to be followed. Other "loiter" type waypoints are interpreted as the centers of a circular path, whose radius should be additionally specified along with the direction of loitering (clockwise or anti-clockwise as seen from the top).

• Custom three-dimensional path were also coded in the autopilot software. These paths are made of a series of continuously connected lines and arcs of circles.

The flight tests were carried out in a large field in the south-east of France, and care was taken during the design of reference paths in order to respect the dimensions of the field, keep a safe distance from inhabited areas and fly in an authorized range of altitude.

The pilot and the ground operator agreed on a test scenario. The pilot had to manually takeoff and land the airplane, and also to switch the flight mode to fully autonomous mode at a suitable time, and switch back to manual control in case of a loss of control or a technical failure. The ground operator was in permanent communication with the pilot, monitoring airspeed and battery level, and detecting anomalies sent to the ground station through telemetry. He also made online parameter changes (tuning) according to visual information of the airplane, telemetry data, and pilot observations. For this second test, the chosen closed reference path (see Fig. 9.8) is similar to the one tested in HITL and consists, for the first part, of a horizontal segment connected to a horizontal half-circle of radius equal to 40m, followed by another horizontal segment.

Results

The second part of the reference path is similar to the first one except that it is inclined with an angle of 15 • w.r.t. the horizontal plane. The commanded airspeed was 10m/s and the path was traveled twice. During this experiments, the weather was relatively calm, with a varying east wind of magnitude less than 5m/s. The vehicle was being flown manually and the controller was switched to autopilot mode at about 80m of the reference path. The vehicle then approached and followed the desired path closely with an error of magnitude ±1m and RMS value of 0.9m (see Fig.

Context of the thesis:

The present work is an extension of prior works on the control of flying vehicles, including rotor vehicles [START_REF] Hua | Introduction to feedback control of underactuated vtol vehicles: A review of basic control design ideas and principles[END_REF], rocket-like vehicles [START_REF] Pucci | Nonlinear feedback control of axisymmetric aerial vehicles[END_REF] and fixed-wing vehicles [START_REF] Hua | A novel approach to the automatic control of scale model airplanes[END_REF]. It contributes to the development of a novel and unified control framework for aerial vehicles.

This "unification" is motivated by the fact that the control community developed specific control designs for the class of fixed-wing aircraft that are different from nonlinear control strategies employed for rotor vehicles such as quadrotors. The principles that these works have in common, and on which this thesis is based, are the following:

• The proposed control design is based on a hierarchical scheme which offers theoretical as well as practical advantages. One of which, is the possibility to exploit the fully actuated attitude dynamics, and consider a lower-order model for control design that is independent of the actuation configuration of a specific vehicle.

• The flight control community has always relied on conventions and approximations of aerodynamic forces developed by aerodynamicists: the aerodynamic forces are classically decomposed along the wind axis, and their corresponding coefficients are considered to be linear in attack and sideslip angles. Those are different conventions than the ones adopted for this work. Indeed, an alternative generic 123 nonlinear model of aerodynamic forces was proposed in [48] [49], for the purpose of integrating it in control design. This model is physically pertinent and equivalent to classical approximations for small attack and sideslip angles. One of its assets is that the proposed aerodynamic coefficients are periodic nonlinear models that do not grow unbounded when the attack and sideslip angles become large, so that they can cope with larger aerodynamic envelopes. This model is nonetheless a simplification of the physical reality since it implicitly neglects the effects of rotational and accelerated motions, as well as the deflections of control surfaces. These assumptions make the system's equations triangular, so that nonlinear hierarchical controllers can be designed without dealing with zero-dynamics.

• For a vehicle that is not subject to aerodynamic forces (apart from the thrust generation mechanism) such as the quadrotor, trajectory tracking control can be designed using "thrust vectoring" control, where the thrust direction and its magnitude play the role of an intermediate control variable. This solution has been extended to other vehicles subject to aerodynamic forces such as spherical bodies, axisymmetric vehicles and fixed-wing aircraft.

Results and contributions:

The main results and contributions of this thesis, can be summarized as follows:

• Attitude control laws are first designed, and serve as the fast inner-loops of the control scheme. The study considers the actuation configuration of a conventional aircraft, and some methods for computing the desired deflection angles of control surfaces are explained. A two-axis adaptation to the case of airplanes that lack rudder control is also proposed, and the case of a flying-wing configuration is considered.

• A trajectory tracking controller is designed, which exponentially stabilizes any reference trajectory belonging to a large set of admissible trajectories. The proposed control design methodology is inspired from previous work [START_REF] Hua | A novel approach to the automatic control of scale model airplanes[END_REF] on the control of scale-model airplanes, but is different in the way the sideslip angle is stabilized to zero. While in [START_REF] Hua | A novel approach to the automatic control of scale model airplanes[END_REF] the control of this angle is decoupled from the control of the thrust direction, a more elegant solution is proposed in this thesis, where a complete desired vehicle's orientation is specified to ensure exponential stability of zero tracking errors. The results of this part are published in [START_REF] Kai | A nonlinear approach to the control of a disc-shaped aircraft[END_REF].

• The solution is then adapted to the path-following problem. It combines kinematic guidance and dynamic control into a single framework. An extended nonlinear 124 analysis is performed that takes into consideration the interconnections between the different subsystems. These results are published in [START_REF] Kai | A unified approach to fixedwing aircraft path following guidance and control[END_REF] and [START_REF] Kai | Design and experimental validation of a new guidance and flight control system for scale-model airplanes[END_REF].

• Many practical aspects are investigated in order to put the path-following solution into practice. In particular, a method to estimate the air-velocity vector is proposed, which is suitable for small UAVs that are not equipped with attack and sideslip angles sensors. The guidance and control strategies were implemented on a flight controller running a real time operating system. It was first extensively tested in realistic HITL simulations, which validated the application of the control algorithm and the embedded system.

• The path-following strategy was finally tested in flight experiments using a small RC airplane. The reported results, involving aggressive maneuvers, confirmed the robustness of the control solution w.r.t. modeling, measurements, and estimation

approximations.

An asset of the resulting nonlinear controllers (trajectory tracking and path-following), is that they are designed to operate in a large spectrum of operating conditions. In particular, they avoid singularities associated with the parametrization of the vehicle's attitude and heading, and they overcome the limitations associated with classical methods based on linearization along trim trajectories. Indeed, they go further in terms of convergence and stability analyses over an extended flight domain and allow performing aggressive maneuvers beyond trim trajectories. The proposed laws are also complemented with bounded integral actions to compensate for inevitable modeling errors.

Challenges, recommendations and possible extensions:

Many challenges were encountered during this project. Some of the system's parameters were difficult to estimate, mainly the BLDC motor's parameters, and the aerodynamic coefficients of the airplane. This has been effectively handled by the robustness of the controller thanks to the integral terms. However, a closer look at the airspeed regulation results shows a degraded performance in some parts of the aggressive maneuvers. A reason for this might be a poor identification of the motor's parameters that led to nonoptimal control gains. This suggests that the thrust control can be enhanced if a better identification of the motor's parameter is done. A thrust test bench can be constructed for this purpose. The flight simulator used for the HITL simulations was sufficiently realistic, which was an enabling factor for the success of the flight experiments. However, it was noticed after comparing with experiments, that the motor-propeller behavior was not similar in the two cases (simulation and experiments). This suggests enhancing the simulator (by creating custom plug-ins for instance).

125

An additional challenge was the necessity to estimate the air-velocity vector v a . The proposed model-based estimator worked well in practice, however its design relies on many assumptions such as the balanced flight condition, which makes its reliability uncertain. Future improvements of this estimator might take advantage of additional sensors that can measure the direction of the air-velocity.

Additionally, it is recommended that other issues such as the actuator's time constants and saturation limits, are dealt with.

The next logical stage of these studies will involve experiments with more demanding flight scenarios, practical implementation of the trajectory tracking solution, and extension to the control design of hybrid vehicles.

Finally, to my

  smart parents, you made sure I can get what life didn't give you, this achievement is yours. 9.1 Demonstration vehicle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Avionics and actuators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

( 1 1. 1 .

 11 UAVs), as the first flight of an unmanned model glider performed by George Cayley in 1804 (England), and the powered flights of unmanned aircraft performed by John Stringfellow in 1848 (England), and Du Temple in 1857 (France). As explained in [36], Historical examples of autonomous flights the development of unmanned aircraft has depended on the availability of three critical technologies: automatic stabilization, remote control and autonomous navigation 1 .

Figure 1 . 1 :

 11 Figure 1.1: Demonstration of the Sperry Gyrostabilizer, Bezons, France, June 1914 (source: HistoricWings.com)

A

  complete autonomous flight was relying, up until then, on limited navigation capabilities. The traveled distance for instance was measured using a wing-mounted anemometer tied to a pneumatically operated subtracting counter composed of clockworks and timing gears, combined with gyrocompasses to measure the attitude and barometers to measure the altitude. Turbulences and aggressive aircraft maneuvers easily perturbed these mechanisms. In the end, long flight times meant hours of accumulating navigational errors, and engineers wanted to develop more accurate navigation systems. This led Charles Stark Draper, a professor at MIT, to develop the first Inertial Navigation System (INS) for aircraft navigation in 1949, which solved many of the previously evoked problems. Satellite constellations emitting radio signals for navigation appeared in the 1960s and would enable later in the 1990s what is known today as the Global

Today's technology

  is enabling advanced and fast development of flight control systems, leading to the modern era of aerial robotics and the development of small non-military UAVs. This is exemplified by the cheap and more capable microprocessors, the availability of low-cost and robust Inertial Measurement Units (IMUs) based on Micro-Electro-Mechanical Systems (MEMS) technology, and the development of battery technologies such as Lithium Polymer (LiPo) batteries.

Figure 1 . 2 :

 12 Figure 1.2: Flame Wheel F450 DJI Quadrotor at I3S Laboratory Figure 1.3: Yamaha R-max Helicopter (source: Wikipedia)

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.4: The Aerosonde UAV, designed to collect weather data, (source: Wikipedia)

Figure 1 . 6 :Figure 1 . 7 :

 1617 Figure 1.6: SkyTote Tailsitter (source: Wikipedia)

Figure 2 . 1 :

 21 Figure 2.1: Airplane basic components, three views of a CESSNA 172 airplane (source: Technika Lotnicza i Astronautyczna nr 9/1983)

Figure 2 . 2 :

 22 Figure 2.2: Inertial and Body-fixed frames

Figure 2 . 3 :

 23 Figure 2.3: Path and Course angles

Figure 2 . 4 :

 24 Figure 2.4: Attack and sideslip angles

Figure 2 Figure 2 . 6 :

 226 Figure 2.5: BLDC Motor

Figure 3 . 1 : 3 . 2 :

 3132 Figure 3.1: Airflow around a flying-wing aircraft Figure 3.2: Wind tunnel test (source: NASA Langley CRGIS)

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Friction shear stress and pressure

Assumption 3 . 1 .

 31 For vehicles flying at low Reynold numbers and low airspeeds, the aerodynamic coefficients are indepedent of the Reynold and Mach numbers.Hence, the steady aerodynamic coefficients vary only with the direction of the air-velocity h a . Equivalently, they are only functions of the attack angle α and the sideslip angle β, and we can write C(.) = C(α, β) and C M (.) = C M (α, β).

Figure 3 . 5 :

 35 Figure 3.5: Aerodynamic forces on an airplane

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Airfoil nomenclature

Figure 3 Figure 3 . 9 :

 339 Figure 3.8: NACA2412 airfoil

Figure 3 . 11 :

 311 Figure 3.10: Drag Coefficient for NACA2412

Figure 3 . 12 :

 312 Figure 3.12: Change of body axis

20 )CFigure 3 . 14 :

 20314 Figure 3.13: Drag approximation for NACA2412

  surface, this flow establishes a circular motion called a vortex that trails downstream of the wing (see figure3.15). Locally at the wing surface, the flow induces a downward component of the wind velocity called downwash, this decreases the effective angle of attack at the leading edge and therefore decreases the global lift on the wing, it also adds an additional drag component called the induced drag, which is induced by the slight deviation of the lift force from the direction of v a , for more details see[3, chapter 5].

Figure 3 . 15 :

 315 Figure 3.15: Vortices on a finite wing

Figure 3 . 16 :

 316 Figure 3.16: Action of ailerons Figure 3.17: Action of elevator Figure 3.18: Action of rudder

35 )

 35 C m,0 represents a pure torque, and xcpC l,α c α takes into account the torque resulting from the displacement x cp of the center of pressure w.r.t. the center of mass (i.e. the abscissa of the point of application of F a ), which depends among others on the design of the horizontal tail. A design with a negative value of x cp , i.e. when the center of pressure is rearward to center of mass, creates a restoring moment on α, and a stabilizing configuration. The elevator effect on the pitching moment is represented with the term C m,δe . And the negative term C m,ω 2 corresponds to a damping torque. The term C m, α represents the transitory effects of the downwash on the rear components of the aircraft.

37 )

 37 The expressions of the vectors C M,0 , C M,α , C M, α and C M,β and the matrices A and B can be easily deduced from equations (3.34),(3.35) and(3.36).

Figure 4 . 1 :

 41 Figure 4.1: Manual control

Figure 4 . 2 : 4 . 3 :

 4243 Figure 4.2: Speed control system Figure 4.3: Attitude control system

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Path following objective

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Trajectory tracking objective

Figure 4 . 8 :

 48 Figure 4.8: SISO linear control

  has many advantages over NDI techniques such as increased robustness and flexibility in the design. It is common to design an hierarchical architecture for these nonlinear controller similarly to what we already presented in the SISO linear control section, except that here the equations are not linearized and nonlinear feedback linearizing control terms are employed. Indeed, the architecture of the flight controller is as follows, where the symbol → signify "controlled by": • Sideslip angle β → Yawing angular rate ω 3 → Rudder action δ r • Airspeed |v a | → Thrust T • Course angle χ → Roll angle φ → Rolling angular rate ω 1 → Ailerons action δ a • Path angle γ → Pitch angle θ or Attack angle α → pitching angular rate ω 2 → Elevator action δ e

  Step 1, equation 4.11: A desired value for the pitch angle θ d is considered as an intermediate control variable for the stabilization of the path angle γ at its desired value. • Step 2, equation 4.12: A desired pitch angular velocity ω 2d plays the role of a control variable for the stabilization of θ -θ d at zero. • Step3, equation 4.13: Finally, the elevator deflection δ e stabilizes the error ω 2 -ω 2d at zero.

5 Control

 5 Model for Control Design F lying vehicles of different classes (see chapter 1) led to the development of different control strategies often based on different control tools (see chapter 4).

First, it is

  assumed that the actuator's proper dynamics are sufficiently fast so that they can be neglected in the first approximation. Small fixed-wing airplanes, for instance, are equipped with electrical brushless motors for the generation of thrust, and servomotors for the actuation of angular positions of the control surfaces. These motors have their own electronic controller modules, which should ensure sufficiently fast response of the actuators with respect to the mechanical dynamics of the aircraft. Hence desired values of T and δ are supposed to be reached instantaneously.Second, and in view of equation(5.4), ω can be modified at will via the choice of the active control torque (Γ(δ)) produced by the control surfaces. One can thus consider the angular velocity ω as an intermediate control input. This corresponds to postponing the study of the inner-loop of the hierarchical design corresponding to equation(5.4).

7 )C

 7 L (α, 0) = c 1 sin(2α) (5.8) For small attack angles, the drag coefficient C D is thus approximately equal to c 0 and the lift coefficient C L is approximately proportional to the attack angle with the coefficient of proportionality given by 2c 1 , and the expression of C Y for small values of β can also be written as C Y c0 β. Those are the classical first-order approximations of the static aerodynamic coefficients at low attack and sideslip angles, that were discussed in chapter 3. For the two models (the classical and the proposed one) to correspond at small attack and sideslip angles, one can choose to set c 1 = C L,α 2 , c 0 = C D0 and c0 = C C,β . A clear advantage of this model with respect to the classical linear approximations

Figure 5 . 1 :

 51 Figure 5.1: Aerodynamic forces on a sphere

Figure 6 . 1 :

 61 Figure 6.1: Attitude control block diagram

Figure 6 . 2 :

 62 Figure 6.2: Flying-wing actuation configuration

17

 17

  ) to converge to zero to obtain the closed loop equation,ã = ξ( p, ṽ, I p) + o(t), lim t→∞ o(t) = 0 (7.20)

Remark 7 . 3 .Remark 7 . 4 .

 7374 The conditions pointed out in corollary 7.1 prevent us from stating a more global stability result. However, the practical stability domain can be quite large because the set where | F | and |v × F | F | | are equal to zero is very "thin". No condition has so far been put on the sign of the thrust intensity T , whereas only positive thrust can be produced for many aircraft. To take this limitation into account and obtain a result similar to Theorem 1, one has to add it as a constraint in the definition of an admissible trajectory, i.e. by further requiring that T r = (c 0 |v r (t)|v r (t) + m(a r (t) -g)) • ı r

30 )

 30 Indeed, let p = p + k I I p applying the previous feedback control yields, p = -k p s at ∆p ( p) -k d s at ∆v ( ṗ)(7.31)

Figure 7 . 1 :

 71 Figure 7.1: Trajectory

Figure 7

 7 Figure 7.2: Horizontal trajectory

Figure 7 . 3 :Figure 7 . 4 :

 7374 Figure 7.3: Vertical trajectory4 To avoid these errors, one must design a three-times differentiable trajectory.

Figure 7 . 6 :Figure 7 . 8 :

 7678 Figure 7.6: Angle of Attack 5

Figure 8 . 1 :

 81 Figure 8.1: Desired path and position errors

Figure 8 . 2 :

 82 Figure 8.2: Control block diagram

From which results the. 12 ) 8 . 1 .

 1281 ultimate boundedness of | İev |. This Proportional-Integral (PI) feedback controller, complemented with a pre-compensation term, is well defined if ı • v = 0 (a nominally satisfied condition) and yields the closedloop equation, d dt e v = -k T,1 e v -k T,2 α e I ev + v a,2 Ō(v a ) • h (8Proposition Assume that ı • v is always larger than some positive number. If the sideslip speed v a,2 converges exponentially to zero, then the application of the thrust control(8.8)

  t. the tangent to the desired path at the point Q. The usefulness of choosing d 1 = d 2 is related to the possibility of imposing different upper bounds upon | ẏ1 | and | ẏ2 |, a feature which may be useful to separate the rates of convergence along the directions ū and ū, and limit the rate of descent or of climb of the aircraft when, for instance, ū is chosen as a vertical vector.

25 )

 25 Note that ´t 0 | ō(s)|ds ≤ sup(|v|) ´t 0 |o(s)|ds and that this integral is thus bounded. Because ṗF = ẏ1 ū + ẏ2 ū, it comes that ẏ = -|v|ȳ + ō = -k 1 D s at ∆ h (y) + ō, with ō denoting the vector of coordinates of ō along the unit vectors ū and ū. The announced ultimate upper bounds of | ẏ(t)| and | ẏi (t)|(i = 1, 2) follow directly from

2

 2 ds is also uniformly bounded. Therefore the integrals ´t 0 | h(s)| 2 ds and ´t 0 (1 -α 1 (s))|z(s)| 2 ds are uniformly bounded. Because the timederivative of | h| is bounded, the first of these integrals would diverge if | h| did not converge to zero. The same reasoning applies to the second integral and leads to the convergence of (1 -α 1 )|z| 2 to zero. This term also converges to 1 -α ∆z (|z|) |z| 2 = |z| -s at ∆z (|z|) |z| when h converges to zero. Because s at ∆z (|z|) < |z| when z = 0, the convergence of this latter term to zero in turn implies the convergence of |z| to zero.

Theorem 8 . 1 (

 81 Convergence). Consider an aircraft whose motion equations satisfy the kinematic equations (5.1)-(5.3) and the Newton dynamic equation(5.2), complemented with the model (5.6) of aerodynamic forces. Given a desired heading vector h * , apply to this system the attitude angular velocity control(8.38), combined with the thrust control defined by (8.5) (8.8) (8.9). Assume that during the flight 1 -γ 1 y 1 -γ 2 y 2 , |v|, |a -ḡ|, |a * -ḡ|, |h • (a -ḡ)|, |ī × v a |,and |ı • v| are always larger than some positive number. Then |v| converges to v * . Provided that θ(0), i.e. the initial angle between the aircraft frame and B, is different from π, the aircraft frame converges to B and the sideslip angle converges to zero. Furthermore, if (h, z) does not converge to the unstable point (-h * , 0), then (h, z) converges to (h * , 0). In this latter case, if h * is given by (8.22)-(8.23),then the path following error p converges to zero and h converges to the desired direction, i.e. sign vu u. Rates of convergence to the desired equilibria are ultimately exponential. Proof. The convergence of |v| to v * when |ı • v| remains larger than some positive number was established in Proposition 8.1. The convergence of the aircraft frame to B when θ(0) = π was established in Proposition 8.4. Therefore, since |v|, and thus |v a |, are assumed to be bounded, and since  converges to , v a,2 = v a •  = v a • ( -) and the sideslip angle converge to zero. From the convergence of ı to ī one deduces from relations (8.34) and (8.35) that ξ := a-ḡ |a-ḡ| -a * -ḡ |a * -ḡ| converges to zero. The convergence of |v| -v * to zero entails the convergence of d dt |v| -v * to zero. Since aa * = ( d dt |v| -v * )h + |v|( ḣωh × h) one then deduces that (a -a * ) • h converges to zero. This latter property combined with the convergence of ξ to zero in turn implies that ((a -ḡ) • h)( 1 |a-ḡ| -1 |a * -ḡ| ) converges to zero. Since |(a -ḡ)

  1 and k T,2 denoting positive gains, and I ev some bounded integral of e v . The simple PI controller obtained by omitting the pre-compensation term T * may be sufficient in practice to bring and maintain e v close to zero. Stabilizing v a,1 instead of |v| implicitly means that the current aircraft's speed that results from airspeed control, coincides with the desired speed, i.e. v * (t) = |v(t)|. This in turn yields to setting v * = d dt |v| in the relation (8.33) that defines the desired acceleration a * . This supposes that d dt |v| is either estimated or measured.8.6.3 Calculation of PWM control inputsThe controller produces four control inputs, namely the desired thrust T and the three deflection angles composing the vector δ. The calculated values of these inputs are communicated to the motors of the actuators as Pulse-Width-Modulated (PWM) signals, that are encoded in most software as setpoint values comprised between 0 and 1 in the case of T (0 for a non-rotating motor and 1 for maximum rotating speed), and between -1 and +1 for the deflection angles. Concerning the transformation between the generated thrust and the encoded setpoint u T of the squared angular velocity of the propeller's motor, although it is in reality not linear according to section 2.4, we have considered a simple linear approximation T k T u T with k T a positive constant gain. The integral action incorporated in the thrust control law allows for the compensation of this modeling error. As for the electrical servo motors acting on the control surfaces, which interpret the PWM inputs as angular setpoints, a scaling has to be done to bring back calculated angle values within the interval [-1, +1], with ±1 corresponding to chosen saturated values of these angles.

Figure 8 . 3 :

 83 Figure 8.3: Switching scenarios

Figure 8 .

 8 Figure 8.4: HITL tools

Figure 8 . 5 :Figure 8 . 6 :

 8586 Figure 8.5: Simulation 1, 3-axis aircraft, trajectory and reference path (with no wind)

Fig. 8 .

 8 Fig. 8.7 shows the time evolution of various variables. From these figures one can pinpoint the descent phase starting at t ≈ 65s and yielding the zeroing of the thrust during the time-period when the requested computed value is negative and the air-velocity component v a,1 exceeds the desired speed.

Figure 8 . 7 :

 87 Figure 8.7: Simulation 1, 3-axis aircraft: (a) thrust setpoint (b) airspeed v a,1 (c) angles of control surfaces

Figure 8 . 8 :Figure 8 . 9 :Fig. 8 .

 88898 Figure 8.8: Simulation 1, 2-axis pitch-roll aircraft, trajectory and reference path (with no wind)

Figure 8 . 10 :Figure 8 . 11 :

 810811 Figure 8.10: Simulation 2, trajectory and reference path (with wind)

Figure 8 . 12 :

 812 Figure 8.12: Simulation 2, 3-axis aircraft, airspeed v a,1 (with wind).

Figure 8 . 15 :Figure 8 . 16 :

 815816 Figure 8.14: Simulation 3, altitude

Figure 8 . 17 :

 817 Figure 8.17: Simulation 4, trajectory and reference path

• 1 |va| 2 17.5 21 17.5 21 ω 1 -ω * 1 ω 2 -ω * 2 .

 1221211122 Attitude control: k ω = 7.0 Control surfaces angles are calculated according to δ ev,r δ ev,l = -Simulation 4: flying-wing and aggressive maneuvers

Figure 8 . 18 :

 818 Figure 8.18: Simulation 4, position error: components of y. The vertical dashed lines indicate the instances when the autopilot switches to another reference path.

Figure 8 . 19 :

 819 Figure 8.19: Simulation 4, airspeed v a,1

Figure 9 . 1 :

 91 Figure 9.1: Timber airplane

Figure 9 . 2 :

 92 Figure 9.2: Avionics

Figure 9 . 3 :

 93 Figure 9.3: Flight operation

First 4 • 15 • 10 •

 41510 flight: waypoint followingThe model parameters used for the control calculations are η a c 0 = 0.006, η a c 1 = 0.52, m = 2 and g 0 = 9.81. As in the case of HITL simulations, the values of c 0 and c 1 were determined using the empirical formulas in sections 3.2.1 and 3.2.2 that relate the aerodynamic coefficients to the aircraft wingspan and wing surface. The control parameters are:• Guidance: k 1 = 0.75, µ = 0.5, d 1 = 1, d 2 = 0.Airspeed control: k T,1 = 1, k T,2 = k T,1 /2, k T,3 = 10, ∆ ev = Heading stabilization: k h,1 = 1, k h,2 = 0.25, ∆ z = 0.6, k z = Attitude control: k ω = 7.0Desired deflection angles were calculated according to δ * = -1 |va| 2 diag([START_REF] Hua | Control of vtol vehicles with thrust-tilting augmentation[END_REF][START_REF] Kai | Design and experimental validation of a new guidance and flight control system for scale-model airplanes[END_REF][START_REF] Hua | Introduction to feedback control of underactuated vtol vehicles: A review of basic control design ideas and principles[END_REF])(ωω * ). The thrust command was calculated according to 8.6.3, with k T = 14.

Figure 9 . 4 :

 94 Figure 9.4: flight 1, horizontal trajectory and reference path. The waypoints are numbered from 1 to 6. The last waypoint (6) is the center of a reference circular path.

Figure 9 . 6 :

 96 Figure 9.6: flight 1, position error: components of y. The vertical dashed lines indicate the instances when the autopilot switches to another waypoint. The numbers of the waypoints are also indicated and correspond to those shown in figure 9.4

Figure 9 . 7 : 15 • 10 •

 971510 Figure 9.7: flight 1, (a) Thrust setpoint (b) Airspeed

Figure 9 . 8 :

 98 Figure 9.8: Flight 2, trajectory and reference path

9. 9 )

 9 . A position error that remains smaller than the wingspan of the vehicle represents a good result for this kind of small and light airplane, particularly sensitive to wind and equipped with low cost sensors. The airspeed was also maintained close to the desired setpoint, with an error of about ±1m/s and RMS value of 0.5m/s (Fig.9.10). These results highlight the robustness of the control w.r.t. model errors and uncertainties.

Figure 9 . 13 :

 913 Figure 9.13: Flight 3, trajectory and reference path

Figure 9 . 14 :

 914 Figure 9.14: Flight 3, position error: components of y

Figure 9 . 16 :

 916 Figure 9.16: Flight 3, (a) Thrust setpoint (b) Airspeed

  sideslip velocity, i.e. ∀t : v r,2 (t) = 0, 2. strictly positive angles of attack, i.e. ∃ 1 > 0, ∀t : v r,3 (t) > 1

	3. The inequality,	∃ 2 > 0, ∀t : |g -	c0 m	|v r (t)|v r (t) -a r (t)| > 2	(7.15)

Definition 7.1. A trajectory p r (t) such that 0 < < |v r (t)| < v max < +∞ is admissible if assuming zero wind velocity, the equilibrium equation g + Fa m + T m ıa r = 0 is satisfied with 1. zero

  Proof. Let r ∈ R 3 denote a local parametrization of the orientation error between the frames B and B, and h e ∈ R 3 the vector of coordinates of hh * in the inertial 8.6 Practical issues8.6.1 Estimation of the air-velocityThe air-velocity vector v a appears in the control expressions and thus needs to be either measured or estimated. For small UAVs, which fly at low speeds, one cannot assume that the inertial velocity v is a good approximation of v a . Also, since these vehicles are usually not equipped with angle of attack and sideslip angle sensors, the air-velocity is typically not measured directly. Nevertheless, the use of a single Pitot tube allows for the measurement of the component v a,1 in the direction of the aircraft fuselage. Then one can build a model based estimate of v a,3 . Indeed, using equation 8.3, one gets, v a,3 = m c0 |va| (g -a) • k. By further assuming that v a,1 is the main component of v a , an estimate of v a,3 is,

	va,3 =	m η a c0 |v a |	(g -â) • k	(8.41)
	Theorem 8.2 (Local exponential stability). Given the model (5.6) of aerodynamic forces,
	consider the control system composed of the kinematic equations (5.1)-(5.3) and the Newton
	dynamic equation (5.2), augmented with the integrators (8.9) and (8.28). In the case of zero
	wind velocity and a constant v * ( = 0), if the assumptions A1-A2 are satisfied and h * is
	a desired heading vector defined by (8.22)-(8.23), then the control inputs (T, ω) defined by
	(7.19),(8.8),(8.9) and (8.38) locally exponentially stabilize the equilibrium (p, v, B, z, I ev ) =
	(q, sign vu v * u, B, 0, 0).

Thanks as well for providing me with this cool template.

Thanks for having corrected english mistakes in the first chapter.

The term navigation refers to the state estimation techniques used to localize the vehicle in its environment.

[START_REF] Anderson | Fundamentals of Aerodynamics[END_REF] A gyrocompass is a non-magnetic compass whose operation is based on a fast-spinning disc, and uses the principle of gyroscopic precession.

Other secondary flight control surfaces may exist on some airplanes. Among others we mention the wing flaps whose role is to increase the lift force when flying at lower airspeeds, the trim system for the horizontal stabilizer, and the spoilers. The use of these actuators is not considered in this thesis.

The aerodynamic forces are commonly expressed in the wind frame.

A rigid body is a solid body in which the distance between any two given points remains constant.

The engines are the main source of energy, however in the presence of wind, the vehicle is capable of harvesting energy from its surrounding air, an example is the glider airplane which is able to achieve sustainable flights in certain conditions without being equipped with any propulsive system.

Mainly oxygen and nitrogen

The freestream quantities are the characteristics of air far from the body.

At standard sea-level the density of air is approximately 1.225kg/m 3 .

At standard sea-level temperature, µ = 1.7894 × 10 -5 (kg) (m)(s)

At least for flying regimes below the stall region as is explained in later sections.

The definition of the normal force coefficient may correspond to -CZ in some references[START_REF] Stevens | Aircraft control and simulation[END_REF]

This is the direction of the air-velocity that results in a pure drag force, i.e. the resulting aerodynamic force is parallel to the air-velocity.

If the wing has a varying chord, it is common to use the mean chord length along the wing as a reference length for the pitching torque.

The dimension of the position vector p ∈ R 3 is equal to three, and since there is a one-to-one correspondence between the orientation of the vehicle and the Special Orthogonal group (SO(3)), the dimension of the orientation is three as well.

We refer roughly by joystick to all elements used to control the surface deflections, as joystick, yokes and their combination with pedals and even radio commands in the case of radio controlled airplanes.

Note that determining a desired heading vector h * is also equivalent to determining desired path (γ * ) and course (χ * ) angles.

It should be noted that in some works on path following, it is considered that the kinematical guidance also involves second order kinematics, so that the output of the kinematical control block is the desired angular velocity ω h * associated with the unitary vector h * , i.e., ω h * = h * × ḣ * , or equivalently the desired values for the variations of the path ( γ * ) and course ( χ * ) angles.

This is the component of the angular velocity along the direction of the speed.

the derivatives Ẋ * are generally computed as functions of the desired trim speed, path angle and turning radius.

In some references, hybrid vehicles are referred to as VTOL vehicles for vertical takeoff and landing vehicles, while in other references VTOLs refer to multicopters and helicopters.

Precise modeling of aerodynamic forces and torques is important for designing an airplane, for simulation and for evaluation of flight performance, which are not the subject of this work.

This is equivalent in practice to running this discrete controller at a high rate, in accordance with its name as a "fast" inner-loop.

we use directly the notation T (t) instead of Tr(t), because the desired thrust is considered to be applied instantaneously.

In more advanced studies, it can be shown that first-order aerodynamic effects due to the physics of the propellers can be taken into account in control design to enhance the performance of the closed-loop system, see for instance[START_REF] Kai | Nonlinear feedback control of quadrotors exploiting first-order drag effects[END_REF].

As long as r is chosen orthogonal to ır.

This figure shows a high angle of attack at the end of the flight. This is because the reference trajectory ends with a decelerating phase and a slow reference velocity. The situation is of course unrealistic because this simulation doesn't take into account the stall phenomenon.

Acknowledgements

NACA National Advisory Committee for Aeronautics.

NASA National Aeronautics and Space Administration.

NDI Nonlinear Dynamic Inversion.

NED North-East-Down.

With this formalism any smooth curve is characterized by an initial point at s = 0, the choice of a parallel transport frame at this point, and the functions γ 1 (s) and γ 2 (s).

These functions are themselves related to the curve curvature κ and torsion τ according to κ = γ 2 1 + γ 2 2 and τ = d ds (arctan( γ 2 γ 1 )). For instance, γ 1 = 0, γ 2 = 0 in the case of a straight line, and γ 1 = 1 r , γ 2 = 0 in the case of a circle with radius r. By definition of Q (point on the curve closest to G), the vector p is perpendicular to the tangent to the curve at Q (see figure 8.1) . It thus belongs to the plane {Q; ū, ū}. Let y = [y 1 , y 2 ] ∈ R 2 denote the vector of non-zero coordinates of p expressed in the basis of the parallel transport frame F = {Q; u, ū, ū}, i.e. p = y 1 ū + y 2 ū with y 1 = p • ū and y 2 = p • ū. The convergence of p to zero is equivalent to the convergence of y to zero, and one can make p converge to zero by considering its variations w.r.t. the reference frame F. We will use the notation pF and ṗF when differentiating the vector p in the reference frame F. We have, ṗ = d dt

and, ṗF = ẏ1 ū + ẏ2 ū (8.17)

. This relation may also be written as,

We also need an expression of ṡ. Recall that by definition of Q, the error p is orthogonal this latter equality and the convergence of |ō| to zero imposed by the boundedness of the integral of this term. The convergence of y and ẏ to zero then follows from re-writing the previous equality as ẏi = -k

The domain of stability is not global because h = -h * is also an equilibrium. The instability of this latter equilibrium, and the convergence to the former equilibrium when h is initially different from -h * , comes from examining the non-increasing cost function 1 -h • h * which has its maximum value (equal to 2)

A more complete solution involves a complementary integral action in charge of compensating for stationary effects of errors in the modeling of the aircraft dynamics that could prevent the convergence of h to h * . Considering that the error vector h would typically converge to a constant vector w.r.t. a frame rotating with the angular velocity frame. Recall that (0, y ) is the vector of coordinates of the position error pq in the frame F. Define x 1 = r, x 2 = (e v , I ev ) , x 3 = (h e , z ) , x 4 = (0, y ) , and x = (x 1 , x 2 , x 3 , x 4 ) ∈ R 14 . Consider the error system

whose origin x = 0 is an equilibrium. It is clear that the (local) exponential stability of this equilibrium is equivalent to the (local) exponential stability property stated in Theorem 8.2. In view of Propositions 8.1-8.4 and their proofs, in the neighborhood of x = 0 this system is an interconnection of sub-systems in the form

for some non-negative constants γ i,j , and the origin of each (isolated) subsystem ẋi = f i (t, x i ) being exponentially stable due to the existence of a positive function

for some positive constants α i and β i . Define the "interconnection" matrix as

Being lower triangular with positive diagonal terms, there exists a diagonal weightmatrix D such that DS + S T D is symmetric positive definite. Then, by application of Theorem 5.4 (page 233) in [START_REF] Hassan | Nonlinear systems[END_REF], x = 0 is locally exponentially stable.

Remark 8.1. Despite the no-wind and constant desired speed assumptions, Assumptions A1-A2 are less conservative than the convergence conditions of Theorem 8.1 from which they derive, because they bear only upon the desired aircraft trajectory and speed. On the other hand, Theorem 8.1 shows that convergence is possible even when the aircraft starts far away from the desired equilibrium. Non-satisfaction of the convergence conditions over a long period of time is a remote possibility in practice, even in the absence of a good control. Nevertheless, this possibility cannot be discarded a priori. It thus matters to take practical precautions and implement control expressions which, besides from being efficient during the time-periods when these conditions are met, yield control inputs that are well defined and bounded in all circumstances. In particular, of course, no division by zero should be allowed.

Application to particular curves 8.5.1 Straight line

C is a straight line passing through the point p c and with constant unit direction vector u. Then γ 1 = γ 2 = 0. Assumption A1 is verified provided that the path is not vertical, i.e. not parallel to the gravitational acceleration. As for Assumption A2, it gives a condition relating the desired speed to the path slope. More precisely it is verified when v * 2 = mg 0 ηac 0 sin(ν), with ν = arcsin(k 0 • u) = -γ (γ is the path angle), i.e. the angle between the path and the horizontal plane. Moreover, the point Q is always unique and its position can be directly calculated from the aircraft position p and the curve characteristics (p c , u). More precisely, q = p c + u.(pp c ) u, p = u × (p -p c ) × u and any pair ( ū, ū) of constant orthonormal vectors perpendicular to u can be used for the control calculations.

Circular path

C is a circle centered at p c with radius r and constant unit vector ū orthogonal to the circle's plane. Note that this plane does not have to be horizontal. Then γ 1 = 1 r , γ 2 = 0. As for the straight line case, the point Q on the curve and the unit vectors ( ū, u) associated with the parallel transport frame (which because of our choice of ū, coincides in this particular case with the Frénet frame) can be directly calculated from the aircraft position p and the curve characteristics (p c , r, ū). More precisely, ū = The condition of positivity of (1 -γ 1 y 1 -γ 2 y 2 ) ensuring the good conditioning of the projection of the aircraft CoM on the circle may also be written as (p -p c ). ū = 0. It is not satisfied only when the aircraft is located on the circle's axis, which corresponds to the case where the aircraft is equidistant to all points on the circle (loss of uniqueness of the closest point). Assumption A1 is always verified, except in the particular case when the circle is vertical and g 0 = v * 2 r . As for Assumption A2, it is verified when v * 2 > mg 0 ηac 0 sin(ν), with ν denoting the angle between the circle's plane and the horizontal plane.

Which clearly coincides with the first unitary vector defining the parallel transport frame.

Hardware-In-The-Loop simulations

The object of this section is to test the path following control approach in Hardware-In-The-Loop (HITL) simulations, which consist in implementing the control algorithms on actual hardware that could equip a physical UAV aircraft. The difference with a real experimentation is that this hardware is connected to a computer simulated model of an aircraft that closely mimics the dynamics of a physical aircraft. This technique is convenient to validate both the control algorithm and the functioning of the embedded system, before moving on to experiments.

Aircraft simulator

We use the X-plane 10 software, a Laminar Research product, as a flight simulator. Xplane implements an aerodynamic model based on the so-called blade element theory (www.x-plane.com/desktop/how-x-plane-works/). This method takes into consideration the geometry of the plane and the different airfoil shapes. It decomposes the wings, and the horizontal and vertical stabilizers into a finite number of elements. Aerodynamic forces are then determined for each element depending on its orientation and the air-velocity at its location on the aircraft. Force calculations take downwash and propwash effects into consideration, with finite-wing corrections depending on the wings geometry. Compressible flow effects are also taken into account. This approach differs from others that traditionally calculate aerodynamic forces from stability derivatives and lookup tables, and whose precision is tied to data acquired from wind tunnels measurements and/or advanced computational fluid dynamics (CFD) simulations. Another asset of this simulation technique, particularly interesting for scale-model UAVs, is that the obtained flight model is not limited to a small flight envelope. One can also design a custom aircraft model by using the "Plane-Maker" tool included in X-plane and create a custom airfoil force model by using the "Airfoil-Maker" tool.

For the present simulations, we have used an available model of a small scale UAV created with the previously cited tools. Some specifications of the model are:

Simulation 3: waypoint following For this simulation, a typical path derived from a sequence of waypoints is designed with the ground control station. Two consecutive waypoints determine a joining line segment and switching to the next segment is done via an acceptance radius of 20m.

This rule applies to all waypoints except for the last one that is interpreted as the center of a horizontal circle of radius equal to 30m. The same wind speed and direction that were used for the second simulation are added to this scenario. Simulation results are reported in Figs. 8.14-8.15.

Results

The first test is a classical waypoint follow-up. The last part of this mission consists in following a circle of radius 40m in the clockwise direction (as seen from above).

Once the airplane reaches the penultimate waypoint, it is at a distance of 125m from the reference circle (and 165m from its center), it is also flying in the opposite direction and has to achieve a tight turning maneuver to track the desired heading direction. The commanded airspeed was 11m/s and the weather was relatively calm. The results in figures 9.4-9.7 show that the airplane managed to approach every part of the path with a root mean square (RMS) position error (|y|) of 0.9m at proximity of the path and RMS airspeed error (|e v |) of 0.5m/s. The airplane also managed to turn back and approach the circular reference path. It kept following the circle until the pilot decided to switch back to manual mode and perform the landing. Note that discontinuities and large transient errors on the position error correspond to switching between different parts of the reference path.

The instants where the airplane was still far away from the reference path and approaching were considered as outliers. The computation of this RMS value only considers the situations where the position error is less than 3m

As in HITL simulations, this switching is performed automatically according to a specified acceptance radius.