A nonlinear global approach to
scale-model aircraft path following control

arXiv:1803.05184v1 [] 14 Mar 2018

Jean-Marie Kai, Tarek Hamel and Claude Samson
Abstract— This paper addresses the path following control
problem for scale-model fixed-wing aircraft. Kinematic guidance and dynamic control laws are developed within a single
coherent framework that exploits a simple generic model of
aerodynamics forces acting on the aircraft and applies to almost
all regular 3D paths. The proposed control solutions are derived
on the basis of theoretical stability and convergence analyses.
They are complemented by addressing several practical issues,
and validated via realistic hardware-in-the-loop simulations.

I. I NTRODUCTION
During the last two decades important progress has been
made in the development of advanced methods for motion
control of autonomous Unmanned Aerial Vehicles (UAVs).
The commercial landscape of these vehicles is characterised
by a plethora of small start-up companies proposing an
increasing variety of devices for various applications such as
power lines inspection, data collecting, crops monitoring, real
estate use of aerial photographs, etc. This rapidly growing
market has in turn boosted research studies in the field of
UAVs, particularly from the Systems and Control community, with specific control strategies developed in relation to
different classes of aerial vehicles, as exemplified by scalemodel airplanes, quadrotors and, more recently, tail-sitters.
A typical aerial vehicle is an underactuated dynamic system
with six main degrees of freedom (for the position and
orientation of the vehicle’s main body) and usually four
(sometimes three) independent control inputs, i.e. a thrust
force along the main direction of the vehicle, and three (or
only two) control torques for body orientation modification.
Associated motion control problems are commonly classified
into two sub-categories, namely trajectory tracking and path
following.
Trajectory tracking refers to the problem of stabilising
a time-parametrised reference position trajectory. It is
best adapted to highly manoeuvrable Vertical-Take-Off-andLanding (VTOL) vehicles (quadrotors, helicopters, ductedfans) [11], [14], [17], [23], convertible vehicles (tail-sitters
in particular) [28] and acrobatic aircraft [26] whose positions
have to be precisely and timely monitored near hovering.
Its applicability to fixed-wing aircraft, although possible
[12], [24], is more limited because the position-timing issue
for cruising vehicles is less essential than following a preplanned path with a given (possibly high) velocity needed
J.M. Kai and T. Hamel are with I3S, Université Côte d’Azur, CNRS,
Sophia Antipolis, France, kai(thamel)@i3s.unice.f r.
C. Samson (corresponding author) is with INRIA and I3S UCACNRS,
Sophia
Antipolis,
France,
claude.samson@inria.f r,
csamson@i3s.unice.f r.

for the production of strong air-lift forces on profiled wings
that dramatically reduce energy expenditure.
Path-following refers to the problem of stabilising a desired forward velocity and of zeroing the ”distance” between
the vehicle position and a desired geometric path. This
”distance” can be characterized by the mathematical distance
between a specific point located on the vehicle body and
the orthogonal projection of this point onto the path [1],
[18], [25], or by considering some variant characterization
to avoid the orthogonal projection non-uniqueness issue in
some particular situations [6], [21]. The path itself may
be attached to a frame that is moving with respect to
some fixed inertial frame [20]. Path following control has
become a central feature of any advanced automated flight
control system. It can be used for waypoints navigation, to
fly a programmed course [1], [6], [18], and for a landing
approach with a given glide slope [5]. It is commonly
decomposed into sub-problems, namely guidance (sometimes
termed as guidance-based control), velocity control, and
attitude control. Guidance essentially consists in determining
the desired heading direction for the aircraft given a path
to follow [18], and involves body kinematic equations only.
It is mostly independent of the aircraft characteristics and
is usually associated with the notion of outer-loop control.
By contrast, velocity and attitude controls take into account
the specificities of the aircraft and of aerodynamic reaction
forces that influence its dynamics. These lower level control
tasks are associated with the notion of inner-loop control.
Automatic monitoring of the aircraft velocity is commonly
done via thrust production and is sometimes called autothrottle control. Automatic monitoring of the aircraft attitude
(or orientation) is in charge, via the production of torques,
of making the actual aircraft heading direction converge to
the direction specified at the guidance level and is typically
performed by modifying the orientation angles of the aircraft
control surfaces. For this latter task three-axis autopilots
monitor angles of the aircraft ailerons (for roll motion),
elevator(s) (for pitch motion), and rudder(s) (for yaw motion). Less sophisticated two-axis autopilots monitor roll and
pitch motions only. The aforementioned separation between
kinematic guidance and dynamic velocity and attitude control
is conceptually attractive and convenient, all the more so that
imprecise knowledge of aerodynamic forces acting on the
aircraft dynamics constitutes the main source of difficulty for
the design of robust controllers. Inner-loop controllers have
historically, and to these days, been essentially designed on
the basis of linearised modelled dynamic equations about socalled trim trajectories (for which the aircraft longitudinal

and angular velocities are constant) with attack and bank
angles often taken as intermediary control variables [3],
[27]. They are also reported in all major flight dynamics
textbooks 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. Among the most
advanced ones let us cite, for instance, those based on
linear control techniques combined with gain scheduling [7],
feedback linearisation combined with LQ control [8], or
Linear Parameter Varying (LPV) modelling and control [9],
[16].
In the present paper we propose a nonlinear control
approach that encompasses within a single framework the
different steps involved in aircraft path following control,
i.e. guidance, velocity control and attitude control. This
approach offers an alternative to current state-of-the-art automated flight solutions, and also to the Total Energy Control
approach [2], [13] centred on the aircraft velocity/altitude
monitoring, which has gained popularity in the field of small
fixed-wing UAV control through its implementation on several open-source autopilot simulation packages. More specifically, the proposed control design methodology exploits a
simple model of aerodynamical forces applied to the vehicle
that is both representative of the physics underlying the
creation of the environmental forces and sufficiently simple
to lend itself to control design and analysis. It complements
our previous work devoted to the trajectory tracking problem
[12], [23] and involves three stages: i) thrust monitoring
of the aircraft velocity, ii) design of a kinematic guidance
vector yielding convergence to the desired path with the prespecified route direction, and iii) determination of a desired
body frame whose asymptotic stabilisation yields a balanced
flight, i.e. sideslip angle zeroing, and convergence of the
aircraft heading vector to the previously defined guidance
vector. This control methodology originally assumes that the
orientation dynamics is fully-actuated. However we show
that it also applies, modulo minor adaptations, to aircraft
possessing only two orientation control axes, i.e. either
pitch-and-roll or pitch-and-yaw, provided that the aircraft is
equipped with a rear vertical stabiliser and, for the pitchand-yaw configuration, with a main wing dihedral.
The paper is organized as follows. Section II recalls the
equations of motion of an aircraft, presents the model used
for the aerodynamic forces acting on an aircraft, and specifies
the control variables used for a generic control design that
does not depend on specific aircraft actuation means. Section
III first recalls the general control objectives associated with
the path following problem and some useful relations related
to the calculation of the distance between the aircraft and the
desired path. The three aforementioned control design stages
are subsequently detailed, yielding control laws that are
theoretically justified via stability and convergence analyses.
Assumptions under which exponential stability is achieved
are then particularized and discussed in the common cases of
straight and circular paths. Complementary practical issues,
namely i) the extension of the control design to curves
defined w.r.t. a translating frame, ii) its adaptation to two-axis

autopilots, iii) the determination of control surfaces angles
from previously determined desired body angular velocity,
and iv) the monitoring of thrust limitations and attack angle
for stall avoidance, are addressed in Section IV. Hardware-inthe-loop simulation results involving a scale-model aircraft
and challenging reference paths, with large initial tracking
errors, imposed positive thrust, and air-velocity measurement
approximations, purposely introduced to test the control
performance and its robustness, are reported in Section V.
Concluding remarks and a short list of possible extensions
of this work are given in the last Section VI.
II. C ONTROL MODEL
A. Notation
Throughout the paper, E 3 denotes the 3D Euclidean
vector space and vectors in E 3 are denoted with bold letters.
Unless specified otherwise the associated reference frame is
the inertial frame with respect to (w.r.t.) which the aircraft
position is defined. Inner and cross products in E 3 are
denoted by the symbols · and × respectively. The following
notation is used.
n
3
• When x ∈ R (resp. x ∈ E ) |x| (resp. |x|) denotes
the Euclidean norm of x (resp. x).
¯ ∆ (y) (∆ > 0, y ∈ Rn ) denotes a twice differ• sat
entiable adaptation, with bounded derivatives, of the
classical vector-valued saturation function sat∆ (y) =
∆
¯ ∆ (y) = α∆ (|y|)y with
)y. More precisely sat
min(1, |y|
α∆ : [0, +∞) 7→ (0, 1] a decreasing twice differentiable
d2 ∆
d ∆
α (0) = dx
function such that α∆ (0) = 1, dx
2 α (0) =
∆
0, α∆ (x) ≤ ∆
,
lim
(α
(x)x)
=
∆.
A
typical
exx→+∞
x
x
tanh(
).
From
these
definitions
ample is α∆ (x) = ∆
x
∆
¯ ∆ (y) ≈ y when |y| is small and |sat
¯ ∆ (y)| ≤ ∆, ∀y.
sat
• G denotes the aircraft center of mass (CoM);
• I = {O; ı0 , 0 , k0 } is an inertial frame;
• B = {G; ı, , k} is an aircraft-fixed frame;
• ω is the angular velocity of B w.r.t. I, i.e.

•
•
•
•
•
•

d
(ı, , k) = ω × (ı, , k)
(1)
dt
m is the body mass;
v is the CoM velocity w.r.t. the inertial frame;
a is the CoM acceleration w.r.t. the inertial frame;
g = g0 k0 is the gravitational acceleration;
vw is the ambient wind velocity w.r.t. I;
va = v − vw is the aircraft air-velocity. The direction
of va in body frame can be characterized by two angles
α and β such that
va = |va |(cos α(cos β ı + sin β ) + sin α k)

(2)

The coordinate vector of any ξ ∈ E 3 w.r.t. the bodyfixed frame B is denoted by the ordinary letter ξ, i.e.
ξ = ξ1 ı + ξ2  + ξ3 k with ξ = [ξ1 , ξ2 , ξ3 ]> .
In the case of an airplane, and with a proper choice of
the body frame unit vectors, α = arcsin(va,3 /|va |) and
β = arctan(va,2 /va,1 ) denote the attack angle and sideslip
angle respectively. This definition slightly differs from the
•

neglect here the dependence of the aerodynamic characteristics on the Reynolds and Mach numbers. This all the more
justified in the case of scale-model aircraft evolving at nearly
constant altitudes and at speeds much lower than the speed
of sound.

Fig. 1.

Frames and forces

one classically used in aeronautics but the two definitions
are locally equivalent for small angles. Note also that, like
any parametrization of the unit sphere by two angles, this
representation is not uniquely defined everywhere. More
precisely β cannot be defined by continuity at α = ±π/2.
B. Aerodynamic forces
The resultant aerodynamic force Fa applied to a rigid body
moving with air-velocity va is traditionally decomposed into
the sum of a drag force FD along the direction of va and a
lift force FL perpendicular to this direction, i.e.
Fa = FD + FL

(3)

The intensities of drag and lift forces are essentially proportional to |va |2 modulo variations characterized by two
dimensionless functions CD and CL , which depend in the
first place on the orientation of va w.r.t. the body, but also
on the Reynolds number Re and Mach number M . These
dimensionless functions are called the aerodynamic characteristics of the body, or drag coefficient and lift coefficient
respectively. More precisely
FD = −ηa |va |CD va , FL = ηa |va |CL va⊥

(4)

with
va⊥ some vector perpendicular to va and such that
|va⊥ | = |va |,
ρΣ
• ηa := 2 with ρ the free stream air density, and Σ an
area germane to the body shape.
Throughout the paper, we assume that the torque produced
by the resultant aerodynamic force Fa can, in practice, be
overcome by a chosen control torque produced, for instance,
by control surfaces. A common assumption in the flight
control literature, when α, β are well defined, is that the
aerodynamic characteristics depend essentially on α and β
(i.e. CD = CD (α, β) and CL = CL (α, β)). Accordingly, we
•

C. Body dynamics
We assume that the control inputs consist of a thrust force
T = T ı applied at G, with T the thrust intensity, and a
torque vector Γ that is essentially produced by the aircraft
moving control surfaces. In other words we assume that
the torque produced by these surfaces easily dominates the
torque produced by fixed surfaces so that this latter torque
can, in the first approximation, be omitted in the modelling
equations of the aircraft. On the other hand, since moving
surfaces are typically much smaller than fixed surfaces,
the relative variation of the aerodynamic forces about the
neutral positions of these surfaces is limited and can also
be neglected in the first approximation. For the sake of
simplification we further assume that the thrust direction
ı is parallel to the aircraft zero-lift plane and that the
corresponding propulsive force does not produce a significant
torque. This direction thus represents the longitudinal motion
direction of the aircraft. The unit vector  can also be
chosen parallel to the zero-lift plane. Being orthogonal to
ı it represents the aircraft lateral motion direction. Under
these assumptions and approximations, the aircraft dynamic
equations are given by (1) complemented with the classical
Newton-Euler equations:
ma = mg + Fa + T ı

(5)

J ω̇ = −S(ω)Jω + Γ

(6)

with J the body inertia matrix, and Γ the vector of coordinates of Γ in the body frame. In view of Eq. (6), ω can be
modified at will via the choice of the control torque Γ so that
one can consider the angular velocity ω as an intermediate
control input. The corresponding physical assumption is that
“almost” any desired angular velocity can be obtained after
a short transient time. This is a standard “backstepping”
assumption. Once it is made, the vehicle actuation consists
in four input variables, namely, the thrust intensity T and
the three components of ω. This assumption allows one to
”eliminate” or, more precisely, postpone the complementary
issue of producing a desired angular velocity via a torque
Γ that can be produced via various and specific physical
means (control surfaces, in the case of a standard airplane).
The body dynamics equations are then reduced to (1) and
(5).
D. Modelling of aerodynamic forces
To be usable for control design, the Newton equation
(5) has to be complemented with a model of the resultant
aerodynamic force Fa . It is well known that the norm of
aerodynamic forces are commensurate with the squared norm
of the air-velocity so that one can safely assume that there
exists two positive numbers c and d such that
|Fa | < c + d|va |2

(7)

sudden orientation changes when the curves straightens out
[4] [10]. 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 ū and
¯ are (arbitrarily) chosen at some point on the curve. The
ū
corresponding variational frame equations are

Fig. 2.

Desired path and position errors

The model that we propose to use is:
Fa = −(c0 va,1 ı + c̄0 va,3 k)|va | + va,2 O(va )

(8)

with c0 and c1 denoting positive numbers, c̄0 = c0 +2c1 , and
O(va ) any Euclidean vector-valued function such that the
a )|
ratio |O(v
is bounded. For instance, in [12] we have used
|va |
O(va ) = −c0 |va | for a model of a disc-shaped aircraft. Relation (8) is compatible with the assumed aircraft symmetry
about the plane (G; ı, k). Note that if the drag coefficient c0
were equal to zero then, in the case of zero sideslip angle
(va,2 = 0), the resultant aerodynamic force would be orthogonal to the zero-lift plane with an amplitude proportional to
sin α|va |2 . This model is also compatible with the previous
relations (3) and (4). Indeed, one easily verifies that, in this
|va |
2
case, va⊥ = cos
α k−tan αva , CD (α) = (c0 +2c1 sin α)/ηa ,
and CL (α) = c1 sin 2α/ηa . For small attack angles |α| the
drag coefficient CD is thus approximately equal to ηc0a and
the lift coefficient CL is approximately proportional to the
attack angle with the coefficient of proportionality given by
2c1
ηa . This is coherent with experimental data performed on
a variety of wing profiles and axisymmetric bodies [23].
However, for lift-optimized wing profiles this model fails to
account for stall phenomena occurring at large attack angles.
A way to tentatively overcome this shortcoming consists in
changing the coefficients c0 and c1 beyond the stall angle.
III. C ONTROL
A. Objectives and useful relations
Consider a three-times differentiable curve C in 3D-space
parametrized by its curvilinear abscissa s and, at the point
Q(s) on this path, an associated parallel transport frame
¯ } with u the vector tangent to the path
F = {Q; u, ū, ū
C at the point Q. An advantage of a parallel transport
frame over the more conventional Frénet frame is that
it is (continuously) well-defined at points where the path
curvature vanishes. It does not suffer from ambiguity and

d
d
d
¯ = −γ2 u ;
¯ (9)
ū
ū = −γ1 u ;
u = γ1 ū + γ2 ū
ds
ds
ds
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
pthe curve curvature
γ12 + γ22 and τ =
κ and torsion τ according to κ =
γ2
d
ds (arctan( γ1 )). For instance, (γ1 = 0, γ2 = 0) in the case
of a straight line, (γ1 = 1r , γ2 = 0) in the case of a circle with
radius r, and (γ1 (s) = κ cos(τ s + c), γ2 (s) = κ sin(τ s + c))
–with c denoting a constant depending on the choice of the
initial frame– in the case of a helix with constant curvature
and torsion.
Given a curve C, the path following control objective
consists in i) stabilising the aircraft speed at a desired value,
and ii) having the aircraft converge to the curve and then
follow it. Let p denote the aircraft CoM position and q
denote the position of the point Q on the curve closest to
G. 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 aircraft 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). Define the position error vector p̃ =: p − q, and
let v ∗ ∈ R+ denote the desired magnitude of the aircraft
velocity. A way to achieve the previously evoked control
objectives consists in stabilising |v| − v ∗ and p̃ at zero.
By definition of Q (point on the curve closest to G) the
vector p̃ is perpendicular to the tangent to the curve at Q. It
¯ }. Let y = [y1 , y2 ]> ∈ R2
thus belongs to the plane {Q; ū, ū
denote the vector of non-zero coordinates of p̃ expressed in
¯ },
the basis of the parallel transport frame F = {Q; u, ū, ū
¯ with y1 =: p̃.ū and y2 =: p̃.ū
¯ . The
i.e. p̃ = y1 ū + y2 ū
convergence of p̃ to zero is equivalent to the convergence
of y to zero. With p̃ taken as a Euclidean vector w.r.t. the
reference frame F one can make this vector converge to zero
by considering its variations w.r.t. to this frame. From now
¨ when the vector
on we will use the notation p̃F , p̃˙ F and p̃
F
p̃ and its derivatives are taken as Euclidean vectors in the
reference frame F. One has
d ~
d ~
− OQ
= v − ṡu
(10)
p̃˙ = OG
dt
dt
and
¯
p̃˙ F = ẏ1 ū + ẏ2 ū
with
ẏ1

˙ ū + p̃.ū˙
= p̃.
= (v − ṡu).ū + p̃.(−γ1 u)
= v.ū

ẏ2

=
=
=

˙ ū
¯ + p̃.ū
¯˙
p̃.
¯ + p̃.(−γ2 u)
(v − ṡu).ū
¯
v.ū

¯ ]> . This relation may also be written
Therefore ẏ = [v.ū, v.ū
as
¯ )ū
¯
p̃˙ F = (v.ū)ū + (v.ū
(11)
= Πu v
with Πu denoting the operator of projection on the plane
orthogonal to u. Forming the second time-derivative of p̃F
yields
¨
p̃
F

=
=

d
¯ )ū
¯
+ dt
(v.ū
¯˙ )ū
¯
˙
¯
(a.ū + v.ū)ū + (a.ū + v.ū
d
dt (v.ū)ū

¯˙ = −γ2 ṡu. Therefore ÿ = [a.ū −
with ū˙ = −γ1 ṡu and ū
¯ − ṡ(u.v)γ2 ū
¯ ]> , and this relation may also
ṡ(u.v)γ1 ū, a.ū
be written as
¨
¯ )ū
¯ − ṡ(u.v)(γ1 ū + γ2 ū
¯)
p̃
= (a.ū)ū + (a.ū
F
(12)
¯)
= Πu a − ṡ(u.v)(γ1 ū + γ2 ū
This latter relation calls for the calculation of ṡ. We have
already established that
ṡu

= v − p̃˙
¯˙ )
= v − (p̃˙ F + y1 ū˙ + y2 ū
˙
= v − p̃F + (γ1 y1 + γ2 y2 )ṡu

Control block diagram

|v|2 (the term −c0 va,1 |va | in the model contributes to this
dissipation), that |v| is itself bounded.
The proposed control design involves the interconnection of
three stages see Fig. 3, namely i) velocity monitoring via
thrust adaptation, ii) desired heading vector determination
(guidance), and iii) attitude control for the balanced flight
stabilisation of the desired heading vector. These three stages
are detailed next.
1) Velocity control via thrust adaptation: The model of
aerodynamic forces (8) may also be written as
Fa = −c̄0 |va |va + 2c1 va,1 |va |ı + mva,2 Ō(va )

T̄
ı + va,2 Ō(va )
m

(16)

c̄0
|va |va
m

(17)

with
ḡ := g −
and

T̄ := T + 2c1 va,1 |va |

Using this relation in (12) yields
(14)

B. Control design
A central quantity involved in the control of any aircraft
v
is the velocity vector v = |v|h with h := |v|
denoting the
heading vector that specifies the direction of displacement
of the aircraft and |v| the inertial speed. In the case of the
path following problem here addressed the decomposition
of v into the product of |v| by h is all the more justified
that convergence to the desired path can be performed at
various speeds with the same heading policy and, vice versa,
that changing the heading policy does not imply modifying
the aircraft speed. Decoupling the |v| control problem from
the control of h is also natural because these problems
essentially involve separate physical inputs: monitoring of
|v| is done via thrust adaptation, i.e. via the control input
T , whereas monitoring of h is performed via the control of
the aircraft attitude, i.e. via the control input ω. The task
of a human pilot is facilitated by this decoupling, just as in
the case of car-driving. In order to avoid useless theoretical
complications, we assume from now on that the thrust T
applied to the aircraft is bounded. Assuming also that the
wind velocity is bounded, this in turn implies, by virtue of
energy dissipation in the air, which increases roughly like

(15)

1
(O(va ) + c̄0 |va |). Combining this relation
with Ō(va ) = m
with (5) yields

a = ḡ +

The scalar product of both members of the previous equality
with u yields
(u.v)
(13)
ṡ =
1 − γ1 y1 − γ2 y2
¨ = Πu a + η
p̃
F
2
¯)
(γ1 ū + γ2 ū
η =: − 1−γ(u.v)
1 y1 −γ2 y2

Fig. 3.

(18)

Let v ∗ denote the desired value of |v| and define err :=
|v| − v ∗ . By using (16) one finds that
d
T̄
err = ḡ.h + (ı.h) − v̇ ∗ + va,2 Ō(va ).h
dt
m
This relation in turn suggests to set

T̄ = m − ḡ.h + v̇ ∗ − kT,1 err − kT,2 Ierr /ı.h

(19)

(20)

with kT,1 and kT,2 denoting positive gains and Ierr some
bounded integral of the velocity error err. This ProportionalIntegral (PI) feedback controller complemented with a precompensation term is well defined if ı.v 6= 0 (a nominally
satisfied condition on which we will come back further) and
yields the closed-loop equation
d
err = −kT,1 err − kT,2 Ierr + va,2 Ō(va ).h
dt
from which the exponential convergence of err to zero
follows provided that the lateral velocity component va,2
itself converges to zero exponentially, i.e. provided that a
”balanced” flight with zero sideslip angle is achieved. How
to control the aircraft attitude so as to achieve this latter
property is addressed later on via the control of the aircraft
attitude. In practice, since no division by zero is allowed,
division by ı.h can be replaced by the multiplication by
|v|(ı.v)/((ı.v)2 +) with  denoting a small positive number.

This rule applies to any division by a term x that may
cross zero. Note also that the gains kT,1 and kT,2 do not
have to be constant. In particular they may be of the form
0
0
kT,j = |ı.v|kT,j
with kT,j
> 0 so as to avoid numerical
problems in the eventuality where |ı.v| becomes very small.
Note finally that an approximation of the above controller
consists in exploiting the fact that the aircraft angle of attack
is nominally small, which is equivalent to assuming that the
vectors ı and h are nearly aligned. With this assumption, and
neglecting the wind velocity one has va,1 = v.ı ≈ v.h = |v|
and ı.h ≈ 1 so that, after simple calculations, relation (20)
yields
T ≈ T ∗ − m(kT,1 err + kT,2 Ierr )
(21)
with T ∗ := −mg.h + c0 |v|2 + mv̇ ∗ . A pure PI controller,
i.e. without the pre-compensation term T ∗ and with large
enough gains, may also in practice be sufficient to monitor
the aircraft speed satisfactorily. Note that the integral action
is not only useful to compensate for the absence of the precompensation term, but also to compensate for the imperfect
knowledge of the physically applied thrust.
Alternatively, instead of |v| one may wish to monitor the airvelocity in the direction of ı, i.e. the component va,1 = va .ı,
which can be measured, for instance, with a pitot tube. Define
now the velocity error err := va,1 −v ∗ . Using (16) the timederivative of err satisfies the relation
d
d
∗
dt err =
dt ((v − vw ).ı) − v̇
= (a − v̇w ).ı + ω.(ı × va ) − v̇ ∗
= (g − v̇w ).ı + ω.(ı × va ) − cm0 |va |va,1
T
+ va,2 Ō(va ).ı
−v̇ ∗ + m
Assuming that va,2 converges exponentially to zero (balanced
flight), exponential convergence of va,1 − v ∗ to zero is then
obtained by setting
T = T ∗ − m(kT,1 err + kT,2 Ierr )

(22)

with T ∗ := −m((g − v̇w ).ı + ω.(ı × va ) − cm0 |va |va,1 − v̇ ∗ ),
kT,1 and kT,2 denoting positive gains, and Ierr some
bounded integral of err. As in the previous case the simple
PI controller obtained by omitting the pre-compensation
term T ∗ may be sufficient in practice to bring and maintain
err close to zero.
2) Desired heading vector determination (guidance): By
viewing the aircraft as a point moving in 3D-space with
a given velocity |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 signvu u with signvu chosen in
advance and equal to either 1 or −1. There are obviously a
multitude of solutions to this problem, a certain number of
which have been reported in the literature [6] [22] [20]. The
solution that we propose here consists in setting
h∗ = sin(θh )l + cos(θh )signvu u

(23)

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
∆h

¯
ȳ := k1 Dsat

(y)/|v|
p
θh : arctan(|ȳ|/ 1 − |ȳ|2 )
¯
ȳ1 ū + ȳ2 ū
l := −
|ȳ|

(24a)
(24b)
(24c)

with D = diag{d1 , d2 } a diagonal matrix with di ∈ (0, 1]
µ|v|
(i = 1, 2), ∆h := k1 max(d
, µ ∈ (0, 1), and k1 a positive
1 ,d2 )
gain. Note that |ȳ| < 1 so that θh is well defined and belongs
to [0, π/2). Note also that sin(θh ) = |ȳ| ≤ µ so that the term
sin(θh )l entering the expression of h∗ is always well defined.
Moreover, if d1 = d2 then sin(θh ) tends to µ when |y| tends
to infinity. Therefore, in this case arcsin(µ) characterizes the
angle of incidence of the desired heading direction w.r.t. the
tangent to the desired path at the point Q when the aircraft
is far from the path. The usefulness of choosing d1 6= d2 is
related to the possibility of imposing different upper bounds
upon |ẏ1 | and |ẏ2 |, a feature which may be useful to separate
¯ and
the rates of convergence along the directions ū and ū
limit the rate of descent or of climb of the aircraft when, for
¯ is a vertical vector as in the case of a linear or
instance, ū
circular horizontal path.
The following proposition summarizes the stability and
convergence properties associated with this desired heading
direction.
Proposition 1 If h = h∗ with h∗ given by (23) and
(24), then (y, ẏ) = (0, 0) is globally asymptotically (locally
exponentially) stable and h tends to signvu u. Moreover
|ẏ(t)| is upper bounded by µ|v| and |ẏi (t)| is upper bounded
di
µ|v| (i = 1, 2).
by max(d
1 ,d2 )
Proof: When h = h∗ the variation of p̃F is given by
p̃˙ F

= Πu v = |v|Πu h∗
= |v|Πu sin(θh )l
= |v||ȳ|l
¯)
= −|v|(ȳ1 ū + ȳ2 ū

¯ , it comes that ẏ = −|v|ȳ =
and, since p̃˙ F = ẏ1 ū + ẏ2 ū
∆h
¯
−k1 Dsat (y). The announced upper bounds of |ẏ(t)| and
|ẏi (t)| (i = 1, 2) follow directly from this latter equality.
The convergence of y and ẏ to zero then follows from
re-writing the previous equality as ẏi = −k1 di α∆h (|y|)yi
d 2
yi = −k1 di α∆h (|y|)yi2 (≤ 0).
(i = 1, 2), so that 0.5 dt
Finally the exponential stability of (y, ẏ) = (0, 0) in the
sense of Lyapunov is simply inherited from the non-saturated
equation ẏ = −k1 Dy that holds when |y| is small. Finally,
since |y| and thus |ȳ| tend to zero, θh converges to zero and,
in view of (23), h (= h∗ ) converges to signvu u.
(end of proof).
Any other guidance law h∗ yielding the same properties
when h = h∗ , i.e. exponential stability of (y, ẏ, ) = (0, 0)
and convergence of h∗ to signvu u, can be used in combination with the attitude controller derived next.

3) Attitude control: Attitude control, as we define it here,
is in charge of making the aircraft heading direction h
converge to the desired one h∗ and of ensuring a balanced
flight, i.e. of zeroing the side-slip angle by zeroing the lateral
velocity component va,2 = va .. We show next that these
two objectives can be achieved via the determination of a
desired mobile frame B̄ := {G; ı̄, ̄, k̄} and the convergence
of the aircraft frame B = {G; ı, , k} to this desired frame.
We assume in this section that |v.ı| is always larger than
some positive number and that thrust control is used to
stabilise |v| at the desired speed v ∗ .
Let us first define the desired unit vector ı̄. From relation
(16) one deduces that
ı=

a − ḡ − va,2 Ō(va )
|a − ḡ − va,2 Ō(va )|

(26)

with ∆z > 0 the chosen upperbound for |z(t)| and kz
denoting a typically large positive number. Note that z is
the integral of the error h̃F as long as |z| does not reach the
saturation value ∆z > 0. A desired angular velocity for h is
then
(28)

¯ , kh,1 a positive gain (not
with zF := z1 u + z2 ū + z3 ū
necessarily constant), kh,2 ≥ 0 a constant gain, and
α1 (h, h∗ , z) := α∆z (|z + h̃F /kz |) (∈ (0, 1])

ωh = ωh∗ + kh,1 h̃ + kh,2 Ih
Forming the time-derivative of (1 − h.h∗ ) yields
− h.h∗ )

d
dt (1

−ḣ.h∗ − h.ḣ∗
−(ωh × h).h∗ − h.(ωh∗ × h∗ )
−(ωh − ωh∗ ).h̃
−(ωh − ωh∗ − kh,2 Ih ).h̃
−kh,2 Ih .h̃
= −kh,1 |h̃|2 − kh,2 Ih .h̃
=
=
=
=

Note that

with ωh∗ := h∗ × ḣ∗ denoting the angular velocity of h∗ ,
h̃ := h×h∗ , and kh a positive gain (not necessarily constant)
whose value determines the rate of convergence of h to h∗ .
The (almost global) exponential stability of h−h∗ = 0 when
d
(1 − h.h∗ ) = −kh |h̃|2 ≤
ωh = ω̄h then results from that dt
0.
A more complete solution involves a complementary integral action in charge of compensating for stationary effects
of errors in the modelling of the aircraft dynamics that could
prevent the convergence of h to h∗ . Define h̃F as the vector
of coordinates of h̃ in the moving frame F, i.e. h̃F =
¯ .h̃]> . Without a complementary integral action,
[u.h̃, ū.h̃, ū
in the absence of wind, and along simple desired paths
corresponding to classical trim trajectories (straight lines,
horizontal circles,...), this vector would typically converge
to a constant vector. To compensate for this type of bias, let
z denote the (bounded) integral of h̃F that we propose to
calculate as follows

¯ ∆z (z + h̃F /kz )
ż = kz − z + sat
(27)
z(0) = 03×1

ω̄h := ωh∗ + kh,1 h̃ + kh,2 α1 (h, h∗ , z)zF

Proof: Define Ih := α1 (h, h∗ , z)zF so that, in view of (28)

(25)

Let ωh denote the angular velocity of h, i.e. ωh = h× ḣ, and
ω̄h denote a ”desired” angular velocity for the heading vector
h that exponentially stabilises h − h∗ = 0 when ωh = ω̄h .
Take, for instance
ω̄h := ωh∗ + kh h̃

Proposition 2 If ωh = ω̄h with ω̄h given by (27)-(28)
then the equilibrium (h, z) = (h∗ , 03×1 ) is almost globally
asymptotically (locally exponentially) stable.

Ih .h̃ = α1 (h, h∗ , z)(z.h̃F )
and also that (27) may be written as
ż = kz − z + α1 (h, h∗ , z)(z + h̃F /kz )



Therefore
d
|z|2
0.5 dt

= −kz (1 − α1 )|z|2 + α1 (z.h̃F )
= −kz (1 − α1 )|z|2 + Ih .h̃

with the arguments of the scalar function α1 omitted for the
sake of legibility. Therefore, by forming the time-derivative
of the positive function V := (1 − h.h∗ ) + 0.5kh,2 |z|2 one
obtains
V̇ = −kh,1 |h̃|2 − 2kh,2 kz (1 − α1 )|z|2 (≤ 0)
This latter relation yields the convergence of |h̃| and |z| to
zero. The exponential stability of the equilibrium (h, z) =
(h∗ , 03×1 ) then follows from the existence (easily verified)
of a positive number c such that V̇ ≤ −cV when |h̃| and |z|
are close to zero.
(end of proof).
Let us then define the ”desired” acceleration
a∗ := v̇ ∗ h + |v|(ω̄h × h)
that we in turn use to define ı̄ as follows
a∗ − ḡ
ı̄ := ∗
|a − ḡ|

(29)

(30)

For the vector ̄ we set
̄ :=

va × ı̄
va × (a∗ − g)
=
|va × ı̄|
|va × (a∗ − g)|

(31)

and the third vector k̄ is just calculated as the cross product
of ı̄ and ̄, i.e.
k̄ := ı̄ × ̄
(32)
The above relations point out that the frame B̄ is well defined
as long as |a∗ − ḡ| and va × ı̄ are different from zero. An
important property is that, like p̃, v, va , 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 ωı̄ := ı̄ × ı̄˙ 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 stabilising B̄ = B is thus wellposed.
Proposition 3 Assuming that the frame B̄ and its angular
velocity ω̄ are well-defined, an angular velocity control
that almost globally asymptotically (locally exponentially)
stabilises B̄ = B is

ω = ω̄ + kω (t) (ı × ı̄) + ( × ̄) + (k × k̄)
(33)
with kω (t) >  > 0.
This is a known result. Let us just mention, for the sake of
completeness, that the proof involves the candidate Lyapunov
function V = 0.5 tan2 ( θ̃2 )), with |θ̃| denoting the angle
between the frames B at B̄, whose time-derivative is found
equal to V̇ = −2kω (t)V (≤ 0) when applying the control
(33).
Proposition 4 If some conditions (specified and discussed
next) are satisfied, then the attitude angular velocity control
(33), combined with the thrust control (18)-(20), ensures the
exponential convergence of the side-slip angle to zero and of
h to h∗ .
Proof: By definition of ̄: va .̄ = 0. Therefore va,2 =
v. = va .( − ̄). The exponential convergence of  to ̄
thus yields the exponential convergence of va,2 (and of the
side-slip angle) to zero. This also means that the objective
of balanced flight is achieved with the angular velocity
control (33). Now, since va,2 exponentially converges to
zero and ı exponentially converges to ı̄, ∗one deduces from
a−ḡ
a −ḡ
relations (25) and (30) that ξ := |a−ḡ|
− |a
∗ −ḡ| exponentially
converges to zero (provided, of course, that |a − ḡ| is always
larger than some positive number). In the section devoted
to thrust control we have also shown that the exponential
convergence of va,2 to zero entailed the convergence of
d
|v| − v ∗ , and subsequently of dt
|v| − v̇ ∗ , to zero when
an adequate thrust control, like (18)-(20), is applied. Since
d
a − a∗ = ( dt
|v| − v̇ ∗ )h + |v|(ḣ − ω̄h × h) one deduces
∗
that (a − a ).h exponentially converges to zero. This latter
property combined with the exponential convergence of ξ
1
to zero in turn implies that ((a − ḡ).h)( |a−ḡ|
− |a∗1−ḡ| )
converges exponentially to zero. Therefore, provided that
|(a − ḡ).h| is always larger than some positive number,
(|a − ḡ| − |a∗ − ḡ|) converges exponentially to zero, and so
does (a−a∗ ). Therefore (ḣ− ω̄h ×h) also converges to zero
exponentially. Since ωh := h × ḣ and ωh .h = 0 one then
infers that ωh converges exponentially to ω̄h and, by a minor
extension of Proposition 4, that h converges exponentially to
h∗ .
(end of proof).
The conditions evoked in Proposition 4 are that |a − ḡ|,
|a∗ − ḡ|, |va × ı̄|, |v.ı|, and |h.(a − ḡ)| must always
be larger than some positive number. The vector u must

also be uniquely defined and its time-derivative u̇ must be
bounded to ensure that ḣ∗ and ωh∗ are themselves welldefined and bounded. In view of the guidance stability
result of Proposition 1 and of the attitude control stability
result of Proposition 3 one can also conclude that, under
the same conditions, the angular velocity (33) ensures the
convergence of the (modelled) aircraft to the desired path,
with the desired direction. The conditions evoked previously
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
linearised equations of the Newton equation (5) complemented with the model (8) of aerodynamic forces are not
controllable. Although they are not satisfied in only very
specific situations, they nonetheless rule out the possibility
of global stability results. However, it remains possible to
state local 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
the first five conditions are satisfied on the desired path if
|u×ı| and |u.ı| are positive (and larger than a small number)
a−ḡ
with
on the path. Since, for a balanced flight, ı = |a−ḡ|
2
c̄
∗2
∗
¯ ) and ḡ = g − m0 v u, these conditions
a = v (γ1 ū + γ2 ū
are themselves satisfied if
¯ − γ2 ū)| > 1 > 0
A1: |g × u − v ∗ 2 (γ1 ū

(34)

and

c̄0 ∗ 2
v − g.u| > 2 > 0
(35)
m
As for the conditions concerning the boundedness of |u̇| and
of |ωh∗ | they are satisfied if
A2: |

A3: |γ1 | and |γ2 | are bounded, and 1−γ1 y1 −γ2 y2 > 3 > 0
(36)
The satisfaction of latter condition also ensures the uniqueness of the projection of the aircraft position on the path. We
can then state the following local stability result
Proposition 5 Consider the Newton equation (5) complemented with the model (8) of aerodynamic forces. In the
case of zero wind velocity and v ∗ constant (6= 0), if
the assumptions A1-A3 are satisfied and h∗ is a desired
heading vector that verifies the properties of Proposition
1 then the control inputs (T, ω) given by (18)-(20) and
(33) locally exponentially stabilise (p̃F , p̃˙ F , z, B, v) =
(0, 0, 03×1 , B̄, signvu v ∗ u).
C. Application to particular curves
1) Straight line: C is a straight line passing through the
point pc and with (constant) unit direction vector u. Then
γ1 = γ2 = 0 so that Assumption A3 in Proposition 5
is verified. Assumption A1 is also 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
0
it is verified when v ∗ 2 6= mg
c̄0 sin(ν), with ν denoting the
so-called 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
(pc , u). More precisely,

q = pc + u.(p − pc ) u, p̃ = u × (p − pc ) × u and any
¯ ) of constant orthonormal vectors perpendicular to
pair (ū, ū
u can be used for the control calculations.
2) Circle: C is a circle centered at pc with radius r and
¯ orthogonal to the circle’s plane. Note
(constant) unit vector ū
that this plane does not have to be horizontal. Then γ1 = 1r ,
γ2 = 0. Alike the straight line case, the point Q on the
curve and the unit vectors (ū, u) associated with the parallel
transport frame (that coincides in this particular case with the
Frénet frame) can be directly calculated from the aircraft
¯ ). More
position p and the curve characteristics (pc , r, ū
¯
¯ ×ū
(p−pc )×ū
¯ , q = pc − rū,

, u = ū × ū
precisely, ū =
¯|
¯ ×ū
| (p−pc )×ū

2

(u.v)
and p̃ = p − q. One also verifies that η = (p−p
ū.
c ).ū

In this case Assumption A3 becomes − 1r (p − pc ).ū >
3 > 0, and 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
∗2
is vertical and g0 = vr . As for Assumption A2, it is verified
0
when v ∗ 2 > mg
c̄0 sin(ν), with ν denoting the angle between
the circle’s plane and the horizontal plane.

IV. C OMPLEMENTARY ISSUES
A. Extension to a curve defined on a translating frame
An extension of the path following problem addressed
previously consists in considering a curve that is fixed w.r.t.
a translating frame whose origin pc moves with velocity vc ,
d
d
i.e. dt
pc = vc , and acceleration ac , i.e. dt
vc = ac . For
the sake of simplification we assume here that this frame
does not rotate, but the extension to a rotating frame is also
possible. A practical application would, for instance, consist
in having an aircraft fly in circles over a moving ground
target. This problem is also referred to as Moving Path
Following (MPF) in [20] [19]. We leave the interested reader
to verify that this extension just involves the modification of
equations (10), (11), (13) and (14) according to
d ~
d ~
d ~
p̃˙ = OG
− OC
− CQ
= (v − vc ) − ṡu
dt
dt
dt
p̃˙ F

¯ )ū
¯
= ((v − vc ).ū)ū + ((v − vc ).ū
= Πu (v − vc )
ṡ =

((v − vc ).u)
1 − γ1 y1 − γ2 y2

¨ = Πu a + η
p̃
F
¯)
η =: −Πu ac − ((v − vc ).u)ṡ(γ1 ū + γ2 ū

(37)
(38)

(39)
(40)

v
The desired heading direction h∗c (= |v|
) must satisfy the
equation
vc
h∗c − |v|
v − vc
= h∗ (= ∗ vc )
|v − vc |
|hc − |v| |

with h∗ given by (23). One verifies that the solution to this
equation is
p
h∗c = Πh∗ vc /|v| + 1 − |Πh∗ vc |2 /|v|2 h∗
(41)
with Πh∗ denoting the operator of projection on the plane
orthogonal to h∗ . For this desired direction to be well defined
it suffices that |v| > |vc |. As expected, h∗c = h∗ when vc =
0.
B. Adaptation to two-axis autopilots
The three-axis (pitch-roll-yaw) control system proposed so
far does not specifically rely on the existence of a stabilising
tail equipped with pitch and yaw control surfaces. All that
matters is the possibility of modifying the aircraft attitude
at will, like multirotor drones that use differential blade
rotation to this purpose. By contrast, two-axis (either pitchroll or pitch-yaw) control systems rely on the existence of
a stabilising tail that compensates, via passive control and
torque coupling, for the lack of an active control torque
about one of the aircraft rotation axes. Typically tail vertical
and horizontal surfaces serve to stabilise the aircraft heading
along the air-velocity direction by creating passive torques
that are opposed to angular variations w.r.t. the air-velocity
direction.
1) Pitch-roll control: In particular, provided that an adequate roll angle is created, the tail vertical surface is
very efficient at maintaining the side-slip angle small. This
explains why active yaw control (often termed as yaw
damping in relation to the comfort it brings to passengers
by allowing for precise side-slip angle zeroing) via the use
of a rotating rudder surface is of secondary importance for
most common airplanes. This also explains the common
use of two-axis pitch-roll autopilots that control airplanes
automatically. With respect to the three-axis desired angular
velocity ω ? given by (33) it thus essentially suffices to create,
via elevator and ailerons actions, pitch and roll torques that
asymptotically stabilise ω1 − ω1? and ω2 − ω2? at zero. This
can be achieved without creating a yaw torque with the tail
rudder.
2) Pitch-Yaw control: This type of two-axis attitude control, involving tail elevator and rudder action only, is more
involved due to the absence of wing ailerons that directly
produce the desired roll torque. In this case roll is indirectly,
and less efficiently, produced by yaw via the action of the
tail rudder in combination with the existence of a main-wing
dihedral and/or a main wing located significantly above the
aircraft CoM. This configuration also requires the existence
of a tail vertical surface that creates yaw stabilisation and
resistance to slideslip.
To be more specific, and give some insight of a possible yawroll control solution, let us consider the case of an airplane
with a (small) dihedral angle θ. By applying the model (8)
of aerodynamic forces to both sides of the main wing with
c0  c̄0 one shows that the aerodynamic roll torque is
approximately given by

Γa,roll ≈ lc̄0 − (θva,2 + lω1 )va,1 + lva3 ω3 ı

with l denoting the lateral distance between the aircraft CoM
and the points of application, on both sides of the wing,
of resultant aerodynamic forces. In the absence of actuated
ailerons this is the main roll torque acting on the aircraft, and
the Euler equation thus becomes in the first approximation


Γa,roll
J ω̇ ≈ −S(ω)Jω +
Γ2,3
with Γ2,3 the vector of directly controlled pitch and yaw
torques. Assuming that J = diag{j1 , j2 , j3 } with j2 ≈ j3 ,
the first line of this vectorial equation, combined with the
expression of Γa,roll , yields

lc̄0
− (θva,2 + lω1 )va,1 + lva3 ω3
(42)
ω̇1 ≈
j1
This relation already indicates that a balanced flight (i.e.
with va,2 = 0) is no longer possible when the aircraft turns.
Indeed, consider for instance the case of a horizontal circular
trajectory followed with a constant speed and a constant roll
angle, i.e. such that ω1 = ω̇1 = 0 and ω3 is constant. Then
(42) yields va,2 ≈ lα
θ ω3 , with α = va,3 /va,1 approximately
equal to the attack angle. This relation also indicates that
va,2 is inversely proportional to the dihedral angle θ.
Another relation of interest concerns the dependence of va,2
upon the other velocity components. In the model (8) of
aerodynamic forces, setting va,2 O(va ) ≈ −c¯0 va,2 |va |, with
c¯0 > 0 denoting the coefficient of resistance to lateral
motion, and using this model in the Newton equation (5),
yields
c¯0
v̇2 ≈ −ω3 v1 + ω1 v3 + g2 − va,2 |va |
m
with g2 = g.. By neglecting the wind velocity and using
|va | ≈ va,1 this relation becomes
c¯0
v̇a,2 ≈ −ω3 va,1 + ω1 va,3 + g2 − va,1 va,2
m
¯
A large resistance coefficient c0 , whose value essentially
depends on the aircraft body profile and the size of the
vertical stabiliser, in turn implies that
m
va,3
g2
va,2 ≈ (−ω3 +
ω1 +
)
c¯0
va,1
va,1
Using this relation in (42) then yields
ω̇1 ≈

lc̄0
mθ
mθ
mθ 
−(
va,3 +lva,1 )ω1 +(
va,1 +lva,3 )ω3 −
g2
j1
c¯0
c¯0
c¯0

and, using the fact that the attack angle (≈ va,3 /va,1 ) is
nominally very small
2

ω̇1 ≈

l c̄0 va,1
mθ
g2 
− ω1 +
(ω3 −
)
¯
j1
c0 l
va,1

This latter relation illustrates the known fact that yaw authority over roll increases with the dihedral angle θ. In view
l2 c̄0 va,1
of the nominally large gain
, if the dynamics of ω1?
j1
is ”enough slow” (a feature that can be enforced via filtering
and saturation), the previous relation suggests to choose
g2
ω3? =
+ kω3 ω1? (kω3 > 0)
(43)
va,1

for the desired value of ω3 , so as to have ω1 approximately
k mθ
converge to ωc3¯0 l ω1? when ω3 converges to ω3? . Accordingly,
c¯0 l
choosing kω3 = mθ
would make the yaw induced roll
velocity ω1 tend to the desired roll angular velocity ω1? ,
except that this choice is only theoretical due to the imprecise
knowledge of the parameters l and (more particularly) c¯0 .
This latter issue is fortunately not very sensitive because ω1? ,
as given by (33), is essentially a feedback control term whose
gain can be multiplied by a positive term without much
affecting the control performance. From there it remains to
produce (via tail rudder action) a yaw torque that stabilises
ω3 − ω3? = 0 asymptotically. The complementary production
of a pitch torque that stabilises ω2 − ω2? at zero then yields
a pitch-yaw control that mimics (up to some limitations) the
pitch-roll control evoked previously.
C. From desired body angular velocity to control-surfaces
angles velocities
Consider a fixed-wing aircraft equipped with standard rollpitch-yaw control surfaces. Let δ = [δ1 , δ2 , δ3 ]> denote
the vector of control-surfaces angles that have to be calculated from the desired body angular velocity ω ∗ determined
previously. The preliminary step relies on the determination/estimation of the matrix-valued function A(va , ω) =:
∂Γ
∂δ |δ=0 that relates angles deviations to torque production for
small angles. This function is aircraft specific. It depends in
particular on the placement of the control surfaces with respect to the aircraft CoM and on their dimensions. Typically,
along classical trim trajectories and nominal air velocities
involving small attack and side-slip angles, A(va , ω) can
be approximated by |va |2 Ā, with Ā a constant matrix. We
assume from now on that A, or at least an estimate of this
matrix-valued function, is known. Since J ω̇ ≈ −S(ω)Jω+Γ
one deduces that a torque equal to Γ? = S(ω)Jω ∗ −kγ J(ω−
ω ∗ ), with kγ a large positive gain, is effective to make ω track
ω ∗ . This in turn suggests to make δ track δ ? =: A−1 Γ? by
setting δ̇ = uδ with uδ =: −kδ (δ − δ ? ) and kδ a positive
gain. However, control surfaces angular velocities, as well
as angles, are usually bounded. Let ūδ,i (i = 1, 2, 3) denote
the saturated value of uδ,i , and ∆δi the maximal value of δi .
A possible driving law for each control surface angle that
respects these bounds is then
∆δi

¯
δ̇i = −k̄δ δi + k̄δ sat

(δi + ūδ,i /k̄δ ) , i = 1, 2, 3

with k̄δ denoting a large positive gain.
D. Thrust bounds and attack angle monitoring
We have so far assumed that the aircraft could produce the
desired thrust T calculated according to (18)-(20). In practice
this desired value of T may leave the physical thrust interval
[Tmin , Tmax ]. When this happens at least one of the control
objectives –i.e. convergence of the aircraft heading direction
h to the desired one h∗ (or h∗c ), or stabilisation 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 circular path with the convergence
dynamics specified by h∗ may require a negative thrust (to

slow down the aircraft) that a common aircraft, for which
Tmin ≈ 0, cannot produce. Similarly, ascending with these
convergence dynamics and velocity may require a thrust that
exceeds Tmax .
To avoid this situation a possibility consists in reducing the
rate of convergence to the desired path in order to have
the control action focused on the stabilisation of the aircraft
velocity about the desired speed v ∗ . This can be done, for
instance, by choosing the parameter µ in the expression of
µ|v|
small enough. Provided
the saturation ∆h := k1 max(d
1 ,d2 )
∗
that v can be maintained on the desired path with the
available thrust (meaning for instance, in the case of a
straight path and assuming that (−signvu mg.u + c0 v ∗ 2 ) ∈
[Tmin , Tmax ]), this strategy can ensure that the calculated
thrust T will ultimately remain in the interval [Tmin , Tmax ]
whatever the initial position and orientation of the aircraft.
In the case of a horizontal path one may alternatively set the
parameter d2 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.
The aforementioned possibility aims at keeping the calculated thrust within the range [Tmin , Tmax ] so that both
control objectives can simultaneously be achieved. However,
one may also temporarily accept an increase of the aircraft
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
Tmin = 0). In this case only the objective of convergence of
the heading direction h to the desired one h∗ is maintained.
This implicitly means that the current aircraft speed obtained
with zero thrust (the saturated value closest to the calculated
negative thrust) coincides with the desired speed, i.e. v ∗ (t) =
d
|v| in the relation
|v(t)|. This in turn yields to setting v̇ ∗ = dt
(29) that defines the desired acceleration a∗ . This supposes
d
that dt
|v| is either estimated or measured. This choice can
also be opted for when thrust control is used to stabilise the
air-velocity component va,1 instead of |v|.
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 Tmax automatically limits the aircraft speed (even
in a deep dive situation). Moreover, when the maximal
thrust is significantly smaller than the aircraft 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) and the crash of the aircraft. The Air France
Airbus tragedy in 2009, between Rio de Janeiro and Paris, is
a sad illustration of this. 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
upperbound for the attack angle. Whenever the predicted
attack angle α∗ := arcsin( |vvaa | .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
va
, ̄ = ̄∗ , k̄ = ı̄ × ̄
(44)
ı̄ = rot(αmax ̄∗ )
|va |
with rot(αmax ̄∗ ) |vvaa | denoting the vector |vvaa | rotated by
the angle αmax about the unit vector ̄∗ . Note that this
choice corresponds to the minimization of |ı̄ − ı̄∗ | under
the constraint arcsin( |vvaa | .k̄) = αmax . In other words the
proposed procedure minimizes the modifications brought
to the calculated desired directions under the constraint
of keeping the attack angle at most equal to the chosen
threshold. Note also that the new attitude directions {ı̄, ̄, k̄},
alike {ı̄∗ , ̄∗ , k̄∗ }, do not depend on the aircraft orientation
so that the stabilisation of these directions via the angular
velocity input ω specified in relation (33) remains a wellposed problem. In this case the pre-compensation velocity
ω̄ is calculated with ωı̄ := ı̄ × ı̄˙ and, in view of (44)
ı̄˙ = αmax ω̄∗ × ı̄ + rot(αmax ̄∗ )Π |vva |
a

v̇a
|va |

The convergence of {ı, , k} to {ı̄, ̄, k̄} then ensures that
the aircraft attack angle is ultimately bounded by αmax .
When implementing this safety procedure (again imposed
by a limited maximal thrust) one must further pay attention
to choosing the desired speed v ∗ larger (with some security
margin) than the horizontal speed that can be sustained when
the attack angle is kept equal to αmax . The reason is that a
smaller horizontal speed calls for a larger attack angle that
is no longer allowed for by the controller. Imposing such a
speed will thus result in the aircraft going down.
V. S IMULATIONS
A. Hardware-in-the-loop
The object of this section is to test and validate via realistic
simulation the path following control approach reported in
the previous sections. A particularity of the proposed simulations is that they are hardware-in-the-loop implemented in
the sense that the control algorithms are implemented on a
piece of hardware that could equip a physical UAV aircraft.
The only difference with a true 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 type of simulation is widely recognized for its
convenience and soundness.
1) Aircraft simulator: We use the X-plane 10 software, a Laminar Research product, as a flight simulator (http://www.x-plane.com/desktop/how-x-plane-works/).
X-plane implements an aerodynamic model based on the socalled blade element theory. This method takes into consideration the geometry of the plane and the different airfoil
shapes. It decomposes the wings, horizontal and vertical
stabilisers 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 applied to 2D
airfoil data 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 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 ”AirfoilMaker” tool. For the present simulations, we have used a
publicly available model of a scale UAV created with the
previously cited tools and based on a PT-60 RC model. Some
specifications of the model are:
•
•
•
•

Wingspan: 1.8m
Wing surface: 0.6m2
Fuselage length: 1.72m
Weight: 3.6 kg

2) Controller hardware and software: A commercially
available ”Pixhawk”, equipped with a 168 MHz ARM Cortex
M4F CPU and 256 KB of RAM, is used as the autopilot
hardware. Our code implementation is based on the opensource PX4 flight stack that runs on top of a NuttX real-time
operating system and uses the PX4 middleware [15]. 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 control algorithms
presented in the paper. However we kept the other preexisting functionalities and, in particular, the EKF state
observer. Flight data is logged on a SD card, allowing for
post-processing and analysis of this data.
3) Ground control station: We use the ”Qgroundcontrol”
software1 as a control station installed on a base computer,
to design missions and change parameters online during
flights. In our simulation setup, Qgroundcontrol connects
the base computer to the ”PIXHAWK” controller via USB
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 measurements are processed by the included EKF state observer
and are used to calculate PWM (Pulse-width-Modulation)
control values. The controller sends these control values to
the ground station, which in turn sends them to the simulator
1 http://qgroundcontrol.com/

Fig. 4.

HITL tools

that uses them as setpoints for the airplane’s propeller angular
velocity and for the control surfaces angles.
4) Calculation of PWM control inputs: Thrust is produced by the airplane’s propeller. A relation classically used
between the propeller angular velocity δt , whose PWM
value is encoded between 0 and 1, and the produced thrust
force is T = kprop ((kmotor δt )2 − |va,1 |2 ), with kprop and
kmotor denoting positive constant gains. For the simulations
we have purposefully neglected the term |va,1 |2 and used
2
= 28. The integral action incorporated in the
kprop kmotor
thrust control law is there to compensate for this type of
modelling error and our imprecise estimation of the gains
kprop and kmotor . As for the servo motors acting on the
control surfaces they receive PWM setpoint signals whose
encoded values ∆i (i = 1, 2, 3) are comprised between −1
and 1. These latter values correspond to maximum deflection
angles of the control surface on both sides of the neutral zero
angle value. Knowing that the maximum deflection angle
for the simulated airplane is about 0.25rad we have used
δi = 0.25∆i to relate the control surface angle δi to the
corresponding value of ∆i .
B. Simulation Results
For these simulations the aircraft inertial position p,
velocity v, attitude (ı, , k), and angular velocity ω are
supposed to be measured (by using on-board GPS and
gyrometers, for instance) and are thus available for control
computations. The air-velocity component va,1 , as measured
by a pitot tube along the aircraft longitudinal direction ı,
is also available and the demanded airspeed for va,1 is
15m/s. For the reported simulations we have assumed that
an a priori unknown wind may be present so that we can
never presume that the air-velocity va is equal to the aircraft
inertial velocity v. Besides va,1 the other components of va
involved in the control expressions can be measured and/or
estimated in various manners by fusing the information
provided by complementary sensors. Now, it is also possible
to derive and use an approximation of va based on the

120
100

setpoint
3 axis
2 axis (pitch−roll)
2 axis (pitch−yaw)

−Z

80
60
40
20
−200
−100

Y

0
100
200
−250

−200

−150

−50

−100

0

50

100

150

200

250

X

Fig. 5.

Trajectory and reference path (with no wind)

40

3 axis
2 axis (pitch−roll)
2 axis (pitch−yaw)

35

position error (m)

aircraft model equations (5) and (8) and on the sole direct
measurement of the component va,1 . Let us further elaborate
on this possibility. From (16), and considering that the nonspecified transversal component va,2 Ō(va ) of the resultant
aerodynamic force is parallel to , one gets a.k = ḡ.k.
Therefore, using the expression (17) of ḡ one deduces that
va,3 = va .k = c̄0m
|va | (g −a).k. Approximating |va | by |va,1 |
and using an estimation â of the aircraft acceleration, an
approximation of va,3 is v̂a,3 = c̄0 |vma,1 | (g−â).k. To estimate
â one may use on-board accelerometers that measure the
gravity acceleration minus the body acceleration, i.e. acc =
g−a so that a = g−acc (= â). Another possibility consists
in estimating a from the measurement of the inertial velocity
v obtained, for instance, with a GPS. Further assuming that
v is almost constant in the body frame yields the crude
estimation â = ω ×v. An even cruder estimation is obtained
by assuming that the aircraft acceleration is small, i.e. â ' 0.
Concerning the sideslip velocity component va,2 , we assume
that it is nominally small by virtue of the strong energy
dissipation along the aircraft transversal direction that results
from the existence of a tail rudder and, to a smaller extent,
from the oblong shape of the body aircraft (this dissipation
is actively re-enforced by the balanced flight control policy).
The resulting approximation of va is then v̂a = va,1 ı+v̂a,3 k.
This solution seems convenient for scale-model UAVs that
are not equipped with dedicated sensors that measure va in
all three directions. We have used it with â = 0 to test the
control robustness.
Besides the three-axis controller developed in the paper
we have tested, for validation and comparison purposes, the
two-axis (pitch-roll and pitch-yaw) adaptations proposed in
Section IV-B. In the pitch-yaw control case we found it
useful to increase the main wing dihedral angle from 3◦
to 5◦ in order to improve the control authority and reach
a comparatively acceptable performance.
The aerodynamic coefficients used for the control calculations are: c0 = 0.014, and c1 = 0.8. These coefficients
correspond to a glide rate equal to 5.4 and a gliding speed
equal to 15.3m/s (see the Appendix).
Two sets of control parameters have been used, depending on
the controller’s type. For the three-axis and two-axis pitchroll controllers the parameters are:
• Guidance: k1 = 1.5, µ = 0.5, d1 = 1, d2 = 0.4
2
• Velocity/airspeed control: kT,1 = 2.5, kT,2 = kT,1 /4
• Heading stabilisation: kh,1 = 3, kh,2 = 2.25, ∆z = 0.3,
kz = 10
• Attitude control: kω = 8.0
For the two-axis pitch-yaw controller the parameters are:
• Guidance: k1 = 0.375, µ = 0.5, d1 = 1, d2 = 1
2
• Velocity/airspeed control: kT,1 = 2.5, kT,2 = kT,1 /4
• Heading stabilisation: kh,1 = 3, kh,2 = 2.25, ∆z = 0.3,
kz = 10
• Attitude control: kω = 8.0, kω3 = 0.75
¯ ∆ (x) = α∆ (|x|)x entering
For the saturation functions sat
the control laws we have implemented, in discrete time,
the classical vector-valued saturation function, i.e. with

30
25
20
15
10
5
0
0

20

Fig. 6.

40

60

time (s)

80

100

120

Position error |y| (with no wind)

∆
α∆ (|x|) = max(1, |x|
). As explained in section IV-C,
control surfaces angles are calculated according to δ ∗ =
k
− |vaγ|2 Ā−1 J(w − w∗ ) = − |va1|2 diag([85, 100, 75])(w − w∗ ).
The term S(w)Jw∗ is purposefully neglected because its
contribution is not important compared to other terms. The
actuators angular velocities are also limited to 1rad/s by
simply monitoring the difference between computed consecutive angle values.
The chosen closed reference curve (see Fig. 5) consists,
for the first part, of a horizontal segment connected to a
horizontal half-circle of radius equal to 50m, 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 is carried out with no wind so
that one can appreciate the controllers performance in this
”ideal” case. Comparative simulations results are reported
in Fig. 5 and Fig. 6. They show that the three controllers
allow the aircraft to approach the desired path and then
follow it closely. As expected, a performance degradation
is nonetheless observable in the case of the pitch-yaw controller. Despite the approximations involved in the estimation
of the air-velocity va , the imperfections of the model used
for the control design, and the noisy state estimates produced
by the embedded EKF observers, the other two controllers

120
100

setpoint
3 axis
2 axis (pitch−roll)
2 axis (pitch−yaw)

80

−Z

nominally maintain 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.
Fig. 7 shows the time evolution of various variables in
the case of the three-axis controller. From these figures one
can pinpoint the descent phase starting at t ≈ 70s and
yielding the zeroing of the thrust during the time-period
when the requested computed value is negative and the
air-velocity component va,1 exceeds the desired speed. For

60
40

20
−200
−100

Y

0
100
200
−250

−200

−150

0

−50

−100

50

100

150

200

250

X

Fig. 8.

Trajectory and reference path (with wind)

20
40

10
0
0

20

40

60

80

100

120

time (s)
30

airspeed (m/s)

3 axis
2 axis (pitch−roll)
2 axis (pitch−yaw)

35

State
Setpoint

20

30
25
20
15
10
5

10
0
0

position error (m)

thrust sp (N)

30

0
0
20

40

60

80

100

20

40

120

60

80

time (s)

100

120

140

time (s)
Fig. 9.

15

Position error |y| (with wind)

angles (deg)

10
5
0
−5

ailerons
elevator
rudder

−10
−15
0

20

40

60

time (s)

80

100

120

Fig. 7. 3-axis aircraft: (a) thrust setpoint (b) airspeed va,1 (c) angles of
control surfaces

the second set of simulations a steady (a priori unknown)
wind of intensity |vw | = 5.14m/s blowing from the North
(corresponding to the X-axis in Fig. 5) has been added to the
X-plane scenario. Despite an inevitable (mild) performance
degradation, Fig. 8 and Fig. 9 show that the proposed
controllers, implemented with the basic air-velocity estimator
evoked previously, continue to operate properly in this case.
VI. C ONCLUDING REMARKS
The present paper provides a detailed development of
a new control approach to the path following problem of
scale-model airplanes. It departs from existing approaches

by proposing a single coherent framework for kinematic
guidance and dynamic control. The approach builds upon
a simple analytic model of aerodynamic forces acting on the
vehicle, and the resulting nonlinear controllers are designed
to operate in a large spectrum of operating conditions.
In particular they overcome the limitations associated with
classical methods based on linearisation of aircraft dynamics
equations along so-called trim trajectories. Reported realistic
hardware-in-the-loop simulations illustrate the proposed control methodology, its applicability to three-axis and two-axis
autopilots, and its robustness w.r.t modelling, measurements,
and estimation approximations. The next logical stage of
this study will involve intensive experimentation on different
scale-model airplanes.
A PPENDIX : AERODYNAMIC FORCES COEFFICIENTS IN
RELATION TO GLIDE RATE AND VELOCITY

The glide rate and corresponding glide velocity of an aircraft can be determined experimentally in various manners.
We here recall the relations between these quantities and the
coefficients c0 and c̄0 = c0 + 2c1 used in the control design
model. Comparison of experimentally measured and modelbased values can in turn be used to estimate coefficients c0
and c1 yielding the best fit.
The glide rate gr is commonly defined as the largest possible
ratio between constant horizontal and vertical speeds of the

aircraft in the vertical plane (i.e. the smallest rate of descent)
when no thrust is applied. In the case of the dynamics model
here considered one easily verifies that
gr =

1 − c0 /c̄0
p
2 c0 /c̄0

In the nominal case where c0 /c̄0  1 this relation simplifies
to
p
gr ≈ 0.5 c̄0 /c0
With the same approximation the corresponding gliding
speed is given by
p
|vgr | ≈ m|g|/(c0 c̄0 )0.25
R EFERENCES
[1] A. Pedro Aguiar, Joo P. Hespanha, and Petar V. Kokotovi. Performance
limitations in reference tracking and path following for nonlinear
systems. Automatica, 44(3):598 – 610, 2008.
[2] Matthew E Argyle and Randal W Beard. Nonlinear total energy control
for the longitudinal dynamics of an aircraft. In American Control
Conference (ACC), 2016, pages 6741–6746. IEEE, 2016.
[3] Randal W. Beard and Timothy W. McLain. Small Unmanned Aircraft:
Theory and Practice. Princeton University Press, Princeton, NJ, USA,
2012.
[4] Richard Bishop. There is more than one way to frame a curve.
82:Mathematical Association of America, 03 1975.
[5] F. Le Bras, T. Hamel, R. Mahony, C. Barat, and J. Thadasack. Approach maneuvers for autonomous landing using visual servo control.
IEEE Transactions on Aerospace and Electronic Systems, 50(2):1051–
1065, April 2014.
[6] Morten Breivik and Thor I Fossen. Principles of guidance-based path
following in 2d and 3d. In Decision and Control, 2005 and 2005
European Control Conference. CDC-ECC’05. 44th IEEE Conference
on, pages 627–634. IEEE, 2005.
[7] L. Carter and J. Shamma. Gain-scheduled bank-to-turn autopilot
design using linear parameter varying transformations. Journal of
Guidance, Control and Dynamics, 27(5):1056–1063, 1996.
[8] B. Escande. Nonlinear dynamic inversion and LQ techniques. In:
Magni JF., Bennani S., Terlouw J. (eds) Robust Flight Control, Lecture
Notes in Control and Information Sciences, vol 224. Springer, Berlin,
Heidelberg, 1997.
[9] G. Duan G. Cai and C. Hu. A velocity-based lpv modeling and control
framework for airbreathing hypersonic vehicle. International Journal
of Innovative Computing, Information and Control ICIC International,
7(5):2269–2281, 2011.
[10] Andrew J. Hanson and Hui Ma. Parallel transport approach to curve
framing. Technical report, 1995.
[11] Minh-Duc Hua, Tarek Hamel, Pascal Morin, and Claude Samson.
Introduction to feedback control of underactuated vtol vehicles: A
review of basic control design ideas and principles. IEEE Control
Systems, 33(1):61–75, 2013.
[12] Jean-Marie Kai, Tarek Hamel, and Claude Samson. A nonlinear
approach to the control of a disc-shaped aircraft. In IEEE CDC 2017
56th IEEE Conference on Decision and Control, 2017.
[13] Anthony A Lambregts. Vertical flight path and speed control autopilot
design using total energy principles. 1983.
[14] Taeyoung Lee, Melvin Leoky, and N Harris McClamroch. Geometric
tracking control of a quadrotor uav on se (3). In Decision and Control
(CDC), 2010 49th IEEE Conference on, pages 5420–5425. IEEE,
2010.
[15] Meier Lorenz, Honegger Dominik, and Pollefeys Marc. Px4: A
node-based multithreaded open source robotics framework for deeply
embedded platforms. In IEEE International Conference on Robotics
and Automation, ICRA’15, pages 6235–6240, 06 2015.
[16] A. Marcos and G. Balas. Development of linear parameter varying
models for aircraft. Journal of Guidance, Control and Dynamics,
27:218–228, 2004.
[17] Daniel Mellinger and Vijay Kumar. Minimum snap trajectory generation and control for quadrotors. In Robotics and Automation (ICRA),
2011 IEEE International Conference on, pages 2520–2525. IEEE,
2011.

[18] Derek R Nelson, D Blake Barber, Timothy W McLain, and Randal W
Beard. Vector field path following for miniature air vehicles. IEEE
Transactions on Robotics, 23(3):519–529, 2007.
[19] T Oliveira. Moving Path Following Control System for Fixed-Wing
Unmanned Aerial Vehicles. PhD thesis, U. of Porto, Faculdade de
Engenharia, April 2017.
[20] Tiago Oliveira, A Pedro Aguiar, and Pedro Encarnação. Moving path
following for unmanned aerial vehicles with applications to single and
multiple target tracking problems. IEEE Transactions on Robotics,
32(5):1062–1078, 2016.
[21] Sanghyuk Park, John Deyst, and Jonathan P How. A new nonlinear
guidance logic for trajectory tracking. In AIAA guidance, navigation,
and control conference and exhibit, pages 16–19, 2004.
[22] Guilherme V Pelizer, Natássya BF da Silva, and Kalinka RLJ Branco.
Comparison of 3d path-following algorithms for unmanned aerial
vehicles. In Unmanned Aircraft Systems (ICUAS), 2017 International
Conference on, pages 498–505. IEEE, 2017.
[23] Daniele Pucci, Tarek Hamel, Pascal Morin, and Claude Samson. Nonlinear feedback control of axisymmetric aerial vehicles. Automatica,
53:72–78, 2015.
[24] Wei Ren and Randy W Beard. Trajectory tracking for unmanned air
vehicles with velocity and heading rate constraints. IEEE Transactions
on Control Systems Technology, 12(5):706–716, 2004.
[25] Claude Samson. Path following and time-varying feedback stabilization of a wheeled mobile robot. In Proceedings of the international
conference on advanced robotics and computer vision, volume 13,
pages 1–1, 1992.
[26] Lars Sonneveldt, ER Van Oort, QP Chu, and JA Mulder. Nonlinear
adaptive trajectory control applied to an f-16 model. Journal of
Guidance, control, and Dynamics, 32(1):25, 2009.
[27] Brian L. Stevens and Frank L. Lewis. Aircraft control and simulation. Wiley-Interscience. J. Wiley and sons, New York, Chichester,
Brisbane, 1992.
[28] Jinni Zhou, Ximin Lyu, Zexiang Li, Shaojie Shen, and Fu Zhang. A
unified control method for quadrotor tail-sitter uavs in all flight modes:
Hover, transition, and level flight. In Proc. of the IEEE/RSJ Intl. Conf.
on Intell. Robots and Syst. IEEE, 2017.

