Geometric Reduced-Attitude Control of Fixed-Wing UAVs

Featured Application: Although the focus in this article is on unmanned aerial vehicles, the geometric reduced-attitude controllers presented apply to all ﬁxed-wing aircraft with fully actuated rotational dynamics. The proposed approach could also be applied to other bank-to-turn vehicles such as missiles. The method can be particularly useful for situations where the vehicle experiences large deviations from the attitude reference. Abstract: This paper presents nonlinear, singularity-free autopilot designs for multivariable reduced-attitude control of ﬁxed-wing aircraft. To control roll and pitch angles, we employ vector coordinates constrained to the unit two-sphere and that are independent of the yaw/heading angle. The angular velocity projected onto this vector is enforced to satisfy the coordinated-turn equation. We exploit model structure in the design and prove almost global asymptotic stability using Lyapunov-based tools. Slowly-varying aerodynamic disturbances are compensated for using adaptive backstepping. To emphasize the practical application of our result, we also establish the ultimate boundedness of the solutions under a simpliﬁed controller that only depends on rough estimates of the control-effectiveness matrix. The controller design can be used with state-of-the-art guidance systems for ﬁxed-wing unmanned aerial vehicles (UAVs) and is implemented in the open-source autopilot ArduPilot for validation through realistic software-in-the-loop (SITL) simulations.


Background and Motivation
In recent years, technology advancements have led to increased use of small unmanned aerial vehicles (UAVs) in civil, commercial, and scientific applications. Fixed-wing UAVs [1], as illustrated in Figure 1, have superior range and endurance when compared to rotary-wing UAVs, which enable applications such as environmental monitoring, search and rescue, aerial surveillance and mapping, and medical transportation [2]. To further develop the field, and enable safe and efficient autonomous operation of UAVs, requires robust autopilots that can handle a range of environmental conditions, including turbulent wind conditions, and operate in the presence of highly uncertain aerodynamics [3].
As underactuated vehicles, conventional fixed-wing aircraft have fewer control inputs than the dimension of their configuration space. One or more propellers provide a thrustforce in the longitudinal direction, but the forces orthogonal to the thrust axis (lift, sideforce) are not directly controllable. Therefore, fixed-wing UAVs have to resort to using guidance schemes [4], where the UAV's geometric path in 3-D space is controlled by specifying course and flight path angle commands to lower-level autopilots [5]. Due to the fact that small fixed-wing UAVs experience winds that are large relative to their operating airspeeds [1], path-following methods [6] are usually preferred over trajectory tracking control [7]. In path following, the goal is to reach and follow a geometric path, but without any temporal constraints. This also deals with performance limitations of trajectory tracking for systems with nonminimum phase characteristics, such as aircraft [8]. See [9,10] for a comparison of different path-following algorithms for fixed-wing UAVs, in two and three dimensions, respectively. For a recent survey with a focus on quadrotor UAVs, see [11].
Guidance and control systems for unmanned vehicles can be integrated, or separated [12]. For integrated guidance and control (IGC) systems, the guidance system and inner-loop autopilot are designed simultaneously, taking cross-coupling effects into account. On the other hand, in separated guidance and control (SGC), inner and outer loops are designed separately, with modularity and cross-platform use in mind [13]. Examples of separate guidance algorithms for fixed-wing UAVs include nonlinear guidance laws [14,15], vector-field path following [16,17] and a guidance law based on nested saturations [18]. In [19], path following is achieved by using an existing commercial inner-loop autopilot but augmented with an L 1 adaptive controller to deal with modeling uncertainty and environmental disturbances. While most guidance algorithms use only kinematic models, an integrated approach is presented in [20] that uses a simple model of the aerodynamic forces acting on the aircraft. Common to all the mentioned approaches, both IGC and SGC, is the reliance on attitude control in the inner-most loop. The rotational dynamics is not considered but rather assumed to be stabilized by some low-level controller. This motivates further research on attitude controllers, specifically tailored towards fixed-wing UAVs. Several different attitude representations have been employed for fixed-wing UAV path following, including Euler angles [21], rotation matrices [22] and unit quaternions [23]. Minimal representations such as Euler angles are often used because of their intuitive interpretation but suffer from "gimbal-lock" singularities [24]. Unit quaternions [25] are singularity-free, but provide a double cover of SO(3), the space of 3-D rotations. This might lead to unwinding, where the UAV unnecessarily makes a full rotation, even when arbitrarily close to the target attitude [26,27]. Rotation matrices, on the other hand, provide a global and unique representation. This has led to a significant research effort into socalled geometric attitude control, where singularity-free controllers are designed directly on SO(3), using rotation matrices, that avoid the unwinding phenomenon and often controls the system along geodesics, i.e., paths of minimum length in rotation space [28][29][30][31][32][33]. These advantages are desirable when the controlled vehicle is subject to large angle rotations, e.g., a fixed-wing UAV recovering from large attitude errors resulting from severe wind gusts [34].
Fixed-wing UAVs use one of two main mechanisms for turning: bank-to-turn, where a lateral acceleration is generated by reorienting the lift-force by rolling/banking the UAV, or skid-to-turn, where turning is achieved by generating a sideslip angle, which in turn generates a lateral force that turns the vehicle [35]. In [36], these methods are combined to reduce lateral distortion of camera images gathered by a fixed-wing UAV. In general, bankto-turn is often preferred over skid-to-turn because for most aircraft the lift force is of orders of magnitude greater than thrust forces [37]. Thus, the course angle, yaw angle, and turn rate of aircraft are not controlled directly, but rather through banked-turn maneuvers. For aircraft in coordinated turns, i.e., with zero sideslip angle, the coordinated-turn equation provides a simple relationship between roll angle and resulting turn rate, and is for this reason often used in autopilot design [1,[38][39][40][41][42], including those used in state-of-the-art open-source autopilots [43,44].
Controllers designed using rotation matrices or quaternions control the full attitude, and therefore cannot be directly applied to fixed-wing aircraft using banked turn maneuvers. One approach could be to feedback the true yaw angle into the desired rotation matrix and as such use a rotation error representation for roll and pitch only. However, this representation is highly redundant, as 9 parameters are used to parametrize a twodimensional subspace. A simpler approach, that does not require the full machinery of working in SO (3), is to consider a reduced-attitude representation, evolving on the twosphere, S 2 ⊂ R 3 [45]. In this space of reduced attitude, all rotations that are related by a rotation about some fixed axis, are considered the same [27]. Control systems with reduced attitude evolving on S 2 have previously been studied in the context of spin-axis [45] and boresight-axis [46] control for satellites, pendulum stabilization [31], path-following control of underwater vehicles [47], thrust-vector control for multirotor UAVs [48,49] and for general rigid bodies [50][51][52]. Controllers developed on S 2 are relatively simple compared to those developed using rotation matrices and require fewer matrix operations.
It is well known that a desired attitude (full or reduced) cannot be globally stabilized using continuous state-feedback control laws [26]. This stems from the topological properties of SO(3) and S 2 , which are compact, boundaryless manifolds that are not diffeomorphic to any Euclidean space. The largest possible attraction basins under continuous feedback are almost global, i.e., excluding a zero-measure set, which corresponds to the stable manifolds of additional unstable equilibrium points [53]. However, global asymptotic stability can be achieved by using tools from hybrid dynamical systems, where hysteresis-based switching ensures that all trajectories converge to the desired equilibrium [49,50,52,[54][55][56][57].

Scope and Contributions
In this paper, we present smooth, nonlinear reduced-attitude controllers for fixed-wing UAVs, in a coordinate-free manner, using a global, singularity-free attitude representation on S 2 . The method applies to UAVs with fully actuated rotational dynamics, e.g., those that are equipped with a full set of control surfaces, such as ailerons, elevator, and rudder. The chosen reduced-attitude representation is independent of the yaw angle and thus enables traditional banked-turn maneuvers. A consequence of this is that the presented approach can be deployed in conjunction with state-of-the-art hierarchical flight control architectures that rely on roll and pitch control in the inner loop, such as [1], and those implemented in open-source autopilots such as ArduPilot [43] and PX4 [44]. Furthermore, no lateral/longitudinal decoupling assumptions are used in the design, allowing the attitude controller to compensate for coupling effects that arise when such assumptions are violated.
The reduced-attitude representation allows for a convenient decomposition of the dynamics and a natural corresponding decoupling of the control objective into two parts: (1) reduced-attitude (roll/pitch) control, and (2) control of the angular velocity about the inertial z-axis (turn rate control). Using Lyapunov theory, almost global asymptotic stability is established for three controllers: one constructed based on an energy-like Lyapunov function, a variation of this based on a backstepping procedure, and lastly an adaptive version of the latter that estimates the net aerodynamic moment caused by the translational dynamics (flow angles). This alleviates the need for expensive flow angle measurement equipment, as well as the knowledge of an accurate aerodynamic model. Furthermore, we show that only a rough estimate of the input matrix is needed to achieve ultimate boundedness. The suitability of the proposed attitude control algorithm is demonstrated in realistic software-in-the-loop simulations.

Related Work
The existing work in the literature that shares the most similarities with this paper can be found in [58][59][60][61], where nonlinear attitude controllers for fixed-wing UAVs are devel-oped using quaternions, and that also use a model of the rotational dynamics. In [58,59], the translational and rotational subsystems are decoupled by estimating the higher-order derivatives of the angle of attack and sideslip angle. This enables controllers for the two subsystems to be designed separately. In [60], a nonlinear PID controller for fixed-wing UAV (full) attitude control is presented. The control law is based on unit quaternions and compensates for aerodynamic coupling effects using integral action. This approach is extended in [61] to apply also to rudderless (i.e., underactuated in attitude) fixed-wing UAVs by using a projection of the quaternion error to a yaw-free subspace. In [62], a gainscheduled attitude controller based on Euler angles is given. An algorithm for automatic tuning is provided, and the control system is verified experimentally in a wind tunnel.
Reduced-attitude control has been extensively applied for thrust-vector control of multirotor UAVs, e.g., [48]. Fixed-wing UAVs on the other hand, are subject to additional aerodynamic forces and moments that make control of such vehicles fundamentally different. Besides, the reduced-attitude representation used in this paper (gravity direction represented in body-fixed frame) is different than the thrust-direction of multirotors (bodyfixed axis represented in the inertial frame). The representation used here is similar to that used to stabilize the inverted equilibrium manifold of the 3-D pendulum in [31,63].
The idea of separately controlling reduced attitude and another variable that is decoupled from the reduced-attitude vector is not new. In [64], the reduced attitude is steered along a geodesic path, while the full attitude is stabilized. In [65,66], the attitude control of a quadrotor is decoupled into thrust-vector control on S 2 , and control of the angle of rotation about the thrust vector. A similar approach is taken in [67] with a control allocation strategy that allows to prioritize reduced-attitude correction over yaw errors. Different rotational error metrics for quadrotor control, defined in terms of both full and reduced attitude are compared in [68]. In [69], a vector-projection algorithm is used for trajectory tracking for an agile fixed-wing UAV (where aerodynamics are dominated by the propeller). The roll angle is decoupled from the reference attitude such that thrust and lift forces can be pointed such that position tracking is achieved. Compared to these works, we simultaneously control reduced attitude and an angular velocity around the reduced-attitude vector.
While the present work employs Lyapunov-based methods to develop lightweight control laws with stability guarantees, other approaches using optimal control algorithms have also been proposed, using deep reinforcement learning [70] and nonlinear model predictive control [71].
Preliminary results of the work presented in this paper have previously been reported in [72], and some initial work towards extending this by applying tools from hybrid control can be found in [73].

Organization of the Paper
The rest of the paper is organized as follows: Section 2 presents some notation and preliminaries on the reduced-attitude representation. The UAV equations of motion are given in Section 3, and the control objective is stated in Section 4 along with the definitions of the error functions used. In Section 5, the controllers for the nominal model are presented. Some robustness considerations are stated in Section 6, where we also give an adaptive version of the backstepping-based control law. The simulation results are presented in Section 7, and some concluding remarks are given in Section 8. All lengthy proofs have been relegated to the appendices.

Preliminaries
In this section, we establish some notation and useful mathematical relations that are used throughout the text, before presenting the reduced-attitude representation.

Notation and Definitions
For a ∈ R and x ∈ R n , let |a| and x = √ x x denote the absolute value and the Euclidean norm, respectively. Positive (resp. non-negative) real numbers are denoted R + (R ≥0 ), and the maximum and minimum eigenvalues of a square matrix A is denoted λ max (A), λ min (A), respectively. The induced 2-norm of a matrix A is A = σ max (A), where σ max (A) is the largest singular value of A. For square, real symmetric positive semidefinite matrices A, λ max (A) = σ max (A). Any square matrix A can be written as the sum of a symmetric and skew symmetric part, A = sym(A) + skew(A), where sym(A) = (A + A )/2 and skew(A) = (A − A )/2. For a symmetric matrix A = A , we have the following inequality for quadratic forms: For any u, v ∈ R 3 , the matrix S(u) = −S (u) ∈ so(3) is the skew-symmetric matrix such that S(u)v = u × v. From properties of the cross product we have S(u)v = −S(v)u, S(u)u = 0 and u S(v)u = 0, which implies that u Au = u sym(A)u for any square matrix A.
We make use of standard right-handed coordinate frames: {n}, a local north-eastdown tangent frame (assumed inertial), and {b}, a body-fixed frame centered at the center of gravity of the UAV, with the x-axis in the longitudinal direction and the y-axis pointing towards the right wing.
The three-dimensional special orthogonal group is the set of three-dimensional rotation matrices, given by where I 3 ∈ R 3×3 is the identity matrix. The two-sphere S 2 ⊂ R 3 is defined by The tangent space at a point x ∈ S 2 can be identified with the vectors that are orthogonal to x: and the normal space N x S 2 is the set of vectors parallel to x: Define the orthogonal and parallel projections Π ⊥ x : R 3 → T x S 2 and Π x : Then, any vector v ∈ R 3 can be written as the sum v = Π ⊥ x v + Π x v.

Reduced-Attitude Representation
Let R ∈ SO(3) be the rotation matrix transforming vectors from {b} to {n}, and let e 3 = [0 0 1] represent the inertial z-direction (direction of gravitational acceleration). We employ the following reduced-attitude representation: which is interpreted as the inertial z-axis, expressed in {b}. By expanding (2) using the roll-pitch-yaw Euler-angle parametrization of R [1], the reduced-attitude vector η can be expressed in terms of the roll angle φ ∈ [−π, π] and pitch angle θ ∈ (−π/2, π/2) as follows: Observe that this particular choice of attitude representation is invariant to changes in the heading/yaw angle ψ. The reduced attitude representation is illustrated in Figure 2, where a section of the sphere corresponding to θ = 0 is shown. Figure 3 shows another section where the aircraft is shown from the side with a possible nonzero roll angle. As shown, the vector η is expressed in the body-fixed frame and points towards the ground. The reduced-attitude representation (2) is the same as the one considered for the 3-D pendulum in [31], but different compared to the one used for thrust-vector control for multirotor UAVs, which is the thrust direction in the inertial-frame [48].
Let ω ∈ R 3 be the angular velocity of the body-fixed frame relative to the inertial frame, expressed in the body-fixed frame. The reduced-attitude vector η satisfieṡ which can be derived from (2) and the relationṘ = RS(ω) [74].
Using (1), we can perform an orthogonal decomposition of the angular velocity ω with respect to η such that Applying this decomposition of ω in combination with (4) giveṡ The parallel component ω is the angular velocity about the inertial z-axis (expressed in the body-fixed frame) and does not influenceη.

Remark 1.
Note that since the two-sphere S 2 is a two-dimensional manifold, in principle two degrees of freedom (DOFs) are sufficient to control reduced attitude. However, since η is fixed in the inertial frame, the two required DOFs (control directions) vary with the orientation of the vehicle and are thus not fixed in {b}. Therefore, we need three actuators to make the reduced attitude fully controllable throughout the configuration space. In this paper, we consider only UAVs with fully actuated rotational dynamics, and at each time instant, we use the remaining DOF to control ω .

UAV Rotational Dynamics
A standard dynamic model of the rigid-body rotational dynamics is given by the Euler equations [1] Jω where J = J > 0 is the inertia matrix and M ∈ R 3 is a vector of applied torques, typically a sum of aerodynamic and propulsion effects. In this respect, we write M = M a + M p , where M a denotes the aerodynamic torque, while M p is caused by a rotating propeller.

Aerodynamics
Aerodynamic forces and moments are in general nonlinear functions that are difficult to model accurately. Identification of parameters for even simple linear models from flight data remains a challenging problem [75,76]. Following [1,77] we define the aerodynamic torque as a function of the angular velocity ω, the body-fixed relative velocity v r ∈ R 3 of the UAV (with respect to the surrounding air mass), and a vector u ∈ R 3 of control surface deflections used to control the attitude of the UAV: Reynolds and Mach number effects are usually ignored for small UAVs moving at airspeeds well below the speed of sound [1].
The airspeed V a ∈ R ≥0 , angle of attack α ∈ [−π, π] and sideslip angle β ∈ [−π, π] are defined by where atan2(y, x) is the four-quadrant inverse tangent. Let ρ, b, c, S ∈ R + be the air density, wingspan, mean wing chord, and wing planform area of the UAV, respectively. A first approximation of the aerodynamic moments that is commonly used in the literature [1,77], and can be useful for control design is given by the control-affine model where the parameters C (·) are dimensionless aerodynamic coefficients.

Propulsion Effects
Let Ω p ∈ R be the rotational speed of the propeller, given in radians per second, and without loss of generality, assume that the propeller thrust axis is aligned with the body-frame x axis. Following [1], for some constant k Ω ∈ R, we write This is a reaction torque caused by the motor of the UAV. Since the motor torque is bounded, we can write M p ≤ c Ω . If the propeller axis is not properly aligned with the x-axis of the body frame, we will get additional small non-zero elements in M p , but the bound still holds for some c Ω . If we also consider a gyroscopic torque (typically small, but sometimes actually used to control aircraft attitude, see Lomcevak maneuver [78]), the subsequent analysis must be adjusted slightly, since the gyroscopic moment also depends on the angular velocity of the UAV. Instead of considering M p as a bounded time-varying exogenous signal, we could then write M p ≤ a + b ω for suitable constants a and b. Additional modeling of complex phenomena generated by the interplay between the main body of the UAV and its propeller (slipstream effects) can be found in [79].

Control-Oriented Model
To summarize, the UAV rotational dynamics can be written In horizontal, level flight, the angular velocity ω is zero. To ensure equilibrium flight ("trim conditions"), define where V * a , v * r , M * p is the trim airspeed, trim relative velocity and trim propeller moment (corresponding to the trim throttle setting), respectively. If u = u trim and ω = 0, thenω = 0 during trimmed flight. Now define which represents the deviation from trimmed flight. We can now combine (8), the rotational dynamics (7) and the reduced-attitude kinematics (4) to obtain the following model that is the basis for control design: The state is represented by (η, ω) ∈ S 2 × R 3 , the control input is u ∈ R 3 , and we consider v r (and thus V a ) as an exogenous bounded input.
To fascilitate control design, we will assume the following: The airspeed V a is strictly positive and bounded with bounded derivative

Assumption 2.
The moment vector ∆(v r , t) and its derivative∆(v r , t) are bounded.

Assumption 3. The control effectiveness matrix B is invertible.
Assumption 4. The damping matrix D satisfies x Dx ≤ 0, ∀x ∈ R 3 .

Remark 2.
Assumption 2 is an assumption on the translational dynamics, which is assumed to affect the rotational dynamics through the exogenous signal v r (may also be considered as "internal dynamics"). In practice, since we are dealing with a physical system, ∆(v r , t) and∆(v r , t) will always be of bounded magnitude. However, since the control input is bounded, we would want these bounds to be relatively small. In particular, during nominal flight, the angle of attack α is usually small, and the lift coefficient is such that a perturbation in α tends to be restored [1]. However, if the stall angle of attack is reached, the slope of the lift coefficient changes such that the α-dynamics might go unstable, which in turn results in a high aerodynamic moment ∆(v r , t).

Remark 3.
A square matrix B corresponds to a fixed-wing UAV that has fully actuated rotational dynamics, i.e., three independent actuators. Now B is invertible if it has full rank. It can be shown that the full rank condition corresponds to primary control coefficients being larger than the coefficients associated with secondary roll-yaw coupling effects. The full rank assumption is therefore reasonable for most common fully actuated control surface configurations.

Remark 4. Assumption 4 is a dissipation assumption and is equivalent to requiring that sym(D)
has nonpositive eigenvalues. In nominal flight conditions, this will be true for most airframes [77] but can be relaxed by using a higher derivative gain (adding damping to the system). See Remark 8.

Error Functions
The goal is to design a state-feedback control law u ∈ R 3 to make the reduced attitude η ∈ S 2 asymptotically track a smooth, time-varying reference η d ∈ S 2 and at the same time drive ω to ω d , where ω d ∈ N η S 2 denotes the desired value of ω , yet to be specified. Furthermore, let the desired reduced-attitude vector η d satisfy the reference model where Assumption 5. The angular velocity references ω ⊥ d , ω d and their derivativesω ⊥ where c Let a smooth configuration error function Ψ : where ν is the angle between η and η d . The function Ψ measures the "distance" between two points η and η d on S 2 , and is clearly positive definite with respect to η = η d . There are two critical points: A minimum when η = η d , and a maximum when η = −η d . In subsequent Lyapunov analysis, Ψ is used as a pseudo-potential energy term in Lyapunov functions.
We proceed by defining the following error vectors: where The error vector e η can be viewed as a gradient vector field on S 2 induced by the potential function Ψ [53]. As e η = |sin ν|, e η vanishes at the critical points of Ψ. The error terms e η and e ω are also compatible in the sense thatΨ = e ω e η , which will cancel with the proportional feedback term defined later when calculating the derivative of a Lyapunov function. The error vector e η is geodesic in the sense that its direction defines an axis of rotation which connects η and η d with the shortest possible curve on S 2 .
Differentiating e η giveṡ where we have used (15), the fact that η × ω = 0 and the identity for any a, b ∈ R (which can be derived using the Jacobi identity of vector cross products). From (10), the derivative of e ω satisfies

Control Objective
From our definition of e ω , Equation (15), note that e ω can be decomposed into two orthogonal parts: This means that as e ω converges to zero, The reduced-attitude error vector e η is zero when η = η d . However, this is also the case when η = −η d . Naturally, this choice of configuration error leads to an additional undesired equilibrium point at η = −η d , but due to the topology of the sphere (it is a compact manifold), this is unavoidable when using continuous feedback [26]. The presence of more than one equilibrium prevents us from designing globally stabilizing feedback laws. A suitable notion of stability in this context is the concept of almost global asymptotic stability.

Definition 1.
An equilibrium solution of a dynamical system is said to be almost globally asymptotically stable if it is asymptotically stable with an almost global domain of attraction, i.e., the domain of attraction is the entire state space excluding a set of Lebesgue measure zero [56].
As we consider continuous feedback on a compact configuration manifold, almost global asymptotic stability is the best possible achievable result [27]. In our setting, if the equilibrium point (η, e ω ) = (η d , 0) is almost globally asymptotically stable, then almost all trajectories converge to it, except for those with initial velocity (depending on the initial configuration error) that are exactly such that ω . This set of initial conditions has a dimension lower than the dimension of the state space, and therefore has measure zero.
We now presicely state the control objective as follows:

Almost Global Reduced-Attitude Tracking
Design a state-feedback control law u such that for almost all e η (t 0 ), e ω (t 0 ), η(t) → η d (t) and e ω (t) → 0 as t → ∞.

Remark 5.
Other configuration error vectors (with corresponding potential functions) on S 2 could be used in place of (14), without changing the general approach considered in this paper. The advantage of using (14) for proportional feedback is that it is simple, smooth, and globally defined. However, there are some performance issues, since for initial reduced-attitudes arbitrarily close to −η d , the control action will be close to zero, and the reduced attitude will stay there for an extended period before converging to the desired reduced attitude. Some alternative error vectors that do not vanish when approaching −η d , but are not defined at this point, are given in [45,51,63].
Before continuing with the controller design, we proceed with a discussion on different design choices for ω d .

Coordinated-Turn Equation
The coordinated-turn equation provides an approximation of the relationship between heading rate and the roll angle during banked-turn maneuvers, and is given by [1], where g is the acceleration of gravity, V a > 0, and the roll angle φ has to satisfy |φ| = π/2.
With ω = [p q r] , the heading rate can be also be written as a function of ω and the Euler angles roll and pitch asψ From (3) and (21) we can relate ω toψ as follows: Although the heading rateψ, given by (21), is not well defined for |θ| = π/2, ω is globally defined. Furthermore, if |θ| = π/2, the body-fixed roll rate p satisfies p =φ −ψ sin θ. Therefore, for |θ| < π/2, Equations (21) and (22) can be combined to obtain Motivated by (23) and the coordinated-turn Equation (20), we propose the following design for ω d that satisfies Assumption 5: where θ d and φ d are consistent with η d (in the sense that (3) is satisfied). Clearly, since (23) is only valid for |θ| = π/2, and (24) contains tan φ d , we restrict the desired reduced-attitude as follows:

Remark 6.
We stress that the mentioned singularities at φ = ±π/2 and θ = ±π/2 are only present for the reference angles. The allowed reference orientations cover most typical flight conditions, except for certain aerobatic maneuvers. The controller design, however, is globally defined, which enables recovery from large reduced-attitude errors, e.g., resulting from large wind gusts.
Alternative design choices for ω d : It is possible to consider some variations of the preceding design of ω d . We now present a few of these options, but leave it as an exercise to the reader to fully explore these possibilities.

•
An alternative to (24) is to define ω d in terms of η d and then project to N η S 2 : The extra term η d η = cos(ν) puts less emphasis on turn coordination when errors in reduced attitude are large. • Equation (24) only satisfies the coordinated-turn Equation (20) asymptotically, as η → η d (and φ → φ d ). One might consider to instead use the actual value of φ instead of φ d , but in the case, we cannot guarantee a priori that ω d and its derivative are bounded. This means that the subsequent stability analysis needs to be adjusted. A pragmatic solution could be to use a saturation function in combination with (24).
To summarize, the expression for the total desired angular velocity ω d in (15) is where ω d ∈ N η S 2 is given by (24), and ω ⊥ d ∈ T η d S 2 . An explicit expression forω d , which is needed in the control law, is given in Appendix A. Equation (24) satisfies Assumption 5 with Remark 7. The coordinated turn Equation (20) has an alternative formulation in terms of the course angle [1], which is often used to perform course control. Course control based on the coordinated-turn equation is thoroughly studied in [41].

Control Laws-Nominal Case
In this section, we present nominal state-feedback control laws assuming perfect knowledge of the rotational dynamics. Two different controllers are presented: one based on an energy-like Lyapunov function, and another based on the backstepping procedure. Although perfect model knowledge is assumed, we do not perform feedback linearization/dynamic inversion, but rather exploit model structure such as skew-symmetry and positive definiteness of matrices. This way, we avoid canceling "good" terms, while other terms are dominated in the stability proof.

Control Design Based on an Energy-Like Lyapunov Function
Proposition 1. Consider the tracking error dynamics (18), and for k p > 0, K d = K d > 0, define the control input as where and the matrix K d is chosen such that for some γ > 0. Then the following holds: (i) There are two closed-loop equilibria, given by (η, e ω ) = (±η d , 0). (ii) The equilibrium (η, e ω ) = (−η d , 0) is unstable.
Proof. See Appendix B.

Remark 9.
The region of exponential convergence to the desired equilibrium point can be made (almost) arbitrarily large by increasing k p ("semi-global" property). However, the region of convergence can never include the unstable equilibrium point and its corresponding unstable manifold [53]. Figure 4 shows a block diagram that illustrates how this controller integrates into a typical guidance, navigation and control (GNC) architecture for a fixed-wing UAV. The references for reduced-attitude, angular velocity, and angular acceleration are generated by some outer-loop guidance controller, and the reduced-attitude control law is combined with a control law for airspeed control, e.g., a PI-controller [1]. The controller uses estimates of the rotation matrix R, the angular velocity ω, as well as the relative velocity v r , which are all made available through a state estimation module. The use of v r is relaxed in Section 6.
The control law (26) is based on proportional action that is proportional to the error term (14), which defines an axis of revolution for the direct, shortest rotation connecting η and η d (forming a geodesic curve on the sphere). This is convenient when dealing with large rotation errors and is a property that is not shared with controllers based on Euler angles. A comparison between a geodesic controller like (26) and one based on Euler angles is presented in [72], which indicates that the geodesic controller spends less control energy than the controller based on Euler angles.

Backstepping Design
A disadvantage of the controller design in the previous section is that the scalar proportional gain k p is restrictive. In this section, we present a backstepping controller that allows for a matrix proportional gain, which gives the flexibility for the control to be more aggressive along certain body-fixed axes, which is important due to geometric and aerodynamic asymmetries of aircraft. In the previous section, the proportional action defines a torque that is aligned with the axis of shortest rotation. The backstepping controller, on the other hand, defines a desired angular velocity that generates a geodesic curve on the sphere.
To this end, define the virtual control signal where κ ∈ R + is a user specified parameter. We will show that ω ⊥ = ϕ(η, η d , ω ⊥ d ) solves the kinematic reduced-attitude tracking problem (see the proof of Proposition 2). Now, introduce the tracking-error signal Note thatω d = ω d − κe η and z can be written as z = e ω + κe η . Due to orthogonality properties, z defined as in (33) has the nice property that as z converges to zero, ω ⊥ → ϕ(η, η d , ω ⊥ d ), which stabilizes the desired reduced-attitude, and at the same time, ω converges to ω d .

Proposition 2.
Consider the tracking error dynamics (18), and for k 1 > 0, K 2 = K 2 > 0, define the control input as where and the matrix K 2 is chosen such that for some γ > 0. Then the following holds: (i) There are two closed-loop equilibria, given by (η, z) = (±η d , 0).

Proof. See Appendix C.
As for the previous design, a statement similar to Remark 8 holds true also here. The control laws (26) and (34) might seem similar at first glance, but by inserting z = e ω + κe η we can rewrite Equation (34) in terms of e ω and ω d : , and the time-dependence is implicit through V a , ω ⊥ d and ω d . Here, the feed-forward part is the same as (28), but the change of variables imposed by the backstepping procedure has introduced a time-varying matrix proportional gain K(t), a time-varying derivative gain, as well as a nonlinear feedback-term −κ 2 S(Je η )e η .

Robustness Considerations
There are a few drawbacks to the controller designs presented in Section 5. In particular, the control laws (26) and (34) require the knowledge of the inertia matrix J, the damping matrix D, the input-matrix B, and the aerodynamic moment ∆(v r , t). In this section, we focus on the control law (34) and state some properties regarding robustness to uncertainty in our model estimates. In addition, an adaptive version of (34) is presented, that provides integral action by estimating ∆(v r , t) under a slowly time-varying assumption. Assumption 7. The aerodynamic moment disturbance ∆(v r , t) is slowly varying, satisfyinġ ∆(v r , t) ≈ 0.

Integral Action
The assumption that ∆(v r , t) is known is particularly restrictive. The aerodynamics of aircraft is highly uncertain. Moreover, the explicit computation of ∆(v r , t) requires the knowledge of the surrounding flow field. Although the airspeed can be measured using a small pitot-static tube, equipment that measures the flow angles α and β is usually not readily available for small UAVs. There exists some available technologies [80], but such equipment can be expensive, too large or too heavy, or just impractical to install on small UAVs that often perform belly landings [81]. Several approaches for flow angle estimation have been proposed in the literature [81][82][83][84], but it remains a challenging problem. Therefore, we focus our attention to instead estimating the aerodynamic moments directly. The control input during trim, u trim can often be quite easily identified during manual flight, so we turn our attention to estimating ∆(v r , t) instead of h r (v r ). This also removes the need for an explicit estimate of M p . (18), and let∆ be an estimate of ∆(v r , t). Define the estimation error∆ ∆ − ∆(v r , t), let K 2 , K 3 be symmetric, positive definite matrices and define the control input as

Proof. See Appendix D.
While the controller in the previous section is of PD type, this is a PID controller with feedforward terms. Integral action removes any steady-state error between the desired and actual angular velocity.

Uncertain Model
Sometimes it is desirable not to include integral action in the inner loops of cascaded control systems since this introduces limitations on achievable bandwidth in the inner loop [1]. Therefore, we focus on a version of the controller that uses a fixed-possibly time-and state-varying, but bounded-disturbance estimate. This estimate does not necessarily equal the true value of ∆. Besides, for some hierarchical flight control loops, the angular velocity reference is not made available for the inner loop. We thus remove the assumption that ω ⊥ d is known in the control design, and redefine the backstepping variable asẑ = ω −ω d , andω d = ω d − κe η . The backstepping procedure in the previous section has provided us with a strict Lyapunov function that can be used to show uniform ultimate boundedness of the solutions of the closed-loop system. To account for model uncertainties, consider again the control law (34) where and∆(v r , t),B,Ĵ,D are estimates of ∆(v r , t), B, J, D, respectively. The termω d is defined as the parts ofω d that do not require reference angular velocities or accelerations:ω The next proposition states that, for sufficiently small model uncertainties, and sufficiently small ω ⊥ d , the solutions are ultimately bounded. This is essentially a local input-tostate stability (ISS) property [85].
To parametrize the model uncertainty, let δB For compactness, we define c J = J + E Ĵ and c D = D + E D .

Proposition 4.
Consider the tracking error dynamics (18) and the perturbed controller (45). Assume that δB satisifes x δBx > 0, ∀x = 0. Then, there exists some gain matrix K 2 = K 2 > 0 such that the matrix sym(δBK 2 ) is positive definite. If the matrix K 2 is chosen such that and if c as defined by Equation (A11) is sufficiently small, then the solutions of the closed-loop system are uniformly ultimately bounded, with an ultimate bound that depends on the controller parameters, the model estimation errors, and the reference velocity bounds.

Proof. See Appendix E.
In essence, the matrix K 2 can be chosen such that the controller is robust to model uncertainties, even when the derivatives of the reduced-attitude reference are not available. However, a necessary condition is that the uncertainty in the input matrix is not too large.
The condition x δBx > 0 implies that the control direction is known up to an error of 90 degrees.

Remark 10.
Global stability of the nominal system is a necessary condition for (global) ISS. Due to the topological obstruction to global stabilization on compact manifolds such as S 2 , a relaxed property of almost (global) ISS has been proposed in [86], and sufficient conditions based on dual Lyapunov techniques (density functions) [87] are given. In [88], a combination of Lyapunov and density functions are used to show almost ISS for systems with rotational degrees of freedom, illustrated using a perturbed nonlinear observer. In [89], a complementary set of tools is given, based on Lyapunov functions and the theory of stable and unstable manifolds of dynamical systems. It is shown that the downward equilibrium of a perturbed pendulum with friction is almost ISS. This is a system that is very much of a similar nature to the one considered in this paper. In [90], robustness on SO(3) is considered in the context of nonlinear complementary filters in the presence of measurement errors. "Divergence" of trajectories on SO(3) is defined as trajectories that converge to the manifold of maximum distance, i.e., the manifold of all rotations of angle 180 degrees. In contrast to [89], only kinematic systems are considered. While almost ISS could probably be shown in our case, the result of [89] only considers a perturbation that is independent of the state. Since in our case, the perturbation is state-dependent, we settle for a local property.

Simulation Results
In this section, simulation results are presented. We show results from an ideal Matlabenvironment, as well as realistic software-in-the-loop simulations where discretization effects and simulated sensor noise are present. In both cases, reduced-attitude references are generated from roll and pitch angle references using Equation (3).

Matlab
Figures 5-10 show the simulation results for the adaptive backstepping controller (40) applied to a simulation model of the Aerosonde UAV [1]. The controller uses perfect estimates of the matrices B, J, and D but no information about ∆(v r , t). While the control surface deflections are controlled by the attitude controller, a PI controller is used to control airspeed using throttle [1]. The airspeed reference is constant and set to 35 m s −1 . The attitude controller parameters are set to κ = 1, k 1 = 1, K 2 = diag(7, 5, 7) and K 3 = diag (40,30,40). During the first 20 s, the reduced-attitude reference is constant, corresponding to φ d = 60 deg and θ d = 15 deg, which might correspond to a sharp, climbing turn. During the last 20 s, we use the time-varying reference A cos(2π f (t − 20)), with amplitude A equal to the initial 20 s, and f = 0.1 Hz for roll and f = 0.08 Hz for pitch. The initial conditions are set to ω(0) = 0, φ(0) = −40 deg and θ(0) = −20 deg. Figure 5 shows the tracking performance in terms of roll and pitch angles, while Figure 6 shows the vector coordinates η ∈ S 2 . The errors converge quickly from large initial values, and the velocity errors are kept close to zero throughout the maneuver. During the latter half of the simulated trajectory, a slightly deteriorated tracking performance is observed in pitch. This also applies to the turn rate, visualized in Figure 7. This is explained by looking at Figure 8: When stabilizing a constant reference, the assumption that aerodynamic moments are slowly-varying applies quite well, and the disturbance estimates converge towards their true values. When tracking a time-varying trajectory, however, this assumption seems to break down, which has a negative impact on tracking performance. This variation seems to be attributed to the variations in the angle of attack seen in Figure 9. Nevertheless, this simulated case study shows adequate tracking performance for both constant and time-varying reference trajectories. The effect of turn-coordination can be seen by observing the sideslip angle in Figure 9. During the first 20 s, the sideslip angle is reduced to zero. During the last 20 s, some variation is seen, but the sideslip angle is still kept at small values (less than 2 degrees). The control surface deflections are shown in the bottom half of Figure 9, and are smooth and well below the saturation limits, which in the simulation is set to ±20 deg. Finally, the start of the maneuver is illustrated as a path on the two-sphere in Figure 10, which is an alternative to showing roll and pitch angles when visualizing reduced-attitude trajectories.

Software-in-the-Loop Simulation
This section showcases the efficacy of the control design via realistic software-in-theloop (SITL) simulations. The controller is implemented in the ArduPilot [43] open-source autopilot framework for fixed-wing UAVs. We simulate our code using ArduPilot's SITL framework, using the JSBSim flight dynamics engine with a model of a SIG Rascal 110.
Roll and pitch reference angles are provided by ArduPilot's guidance system. Lateral guidance is performed using a nonlinear guidance law [14,15], often called L1 guidance, by commanding a lateral acceleration a s cmd using where V g is the ground speed of the UAV, L 1 is the distance to a reference point on the desired path, ahead of the UAV, K L 1 is a tuning parameter, and ϕ is the angle between the ground speed vector and the L 1 vector pointing from the UAV to the reference point on the path. The desired lateral acceleration is then converted into a desired roll angle φ d using a simplified version of the coordinated turn equation: The pitch angle reference is calculated using the total energy control system (TECS) [91]. TECS is based on energy principles, and accounts for dynamic coupling in the longitudinal dynamics of the aircraft by simultaneously controlling altitude and airspeed using pitch and throttle. The pitch angle is used to control the energy distribution E D , i.e., the difference between (specific/per mass) potential and kinetic energy, given by where h is the altitude. For a desired altitude h d ∈ R and desired airspeed V a,d ∈ R + , define the desired specific energy distribution E D,d = gh d − V 2 a,d /2 and the errorẼ D = E D,d − E D . Then, the pitch reference is prescribed as follows: where k i ∈ R + , i = 1 . . . 4 are tuning gains. Figures 11-14 shows the results of a simulation run of the adaptive backstepping controller (40). As the derivatives of the reduced-attitude reference are not available in the ArduPilot code, we use a version of (40) where no information about the angular velocity reference is used in the feedforward part of the controller. In addition, up to 20 percent uncertainty is added to all elements of the matrices J, B and D. This makes the controller more akin to (45), but with added integral action. The controller parameters are set to κ = 2, k 1 = 10, K 2 = diag (5,7,5) and The simulated UAV is tasked with following a square pattern, shown in Figure 11. The actual horizontal position and altitude are shown, from takeoff and until a few rounds have been completed. It is clear that the proposed reduced-attitude controller successfully integrates into the ArduPilot infrastructure. Roll and pitch responses are shown in Figure 12, while Figure 13 shows the vector coordinates η. The UAV tracks the reference well, except when there are large steps in roll angle going into a turn. This is where a feedforward from the angular velocity reference could help reduce the errors. Anyhow, the errors are relatively small and do not interfere with the overall control objective. Figure 11 shows that the altitude is kept approximately constant at 100 m, with only minor drops in altitude during sharp turns. The control input is shown in Figure 14. Except for some large spikes in the control surface deflections (due to large steps in roll reference when going into sharp turns), the control input is well behaved.

Conclusions
Nonlinear reduced-attitude controllers for fixed-wing UAVs have been proposed, using geometric methods on the unit two-sphere. The attitude representation is singularityfree, independent of the yaw/heading angle and allows the UAV to perform bankedturn maneuvers while simultaneously tracking a turn rate satisfying the coordinatedturn equation. Using an aerodynamic model of the rotational dynamics, almost global asymptotic stability is established for the proposed controllers.
The suitability of the presented approach has been verified using Matlab-simulations as well as more realistic SITL simulations, where the control law shows that it successfully completes the defined control objectives in the presence of uncertain aerodynamics and reference velocities, and integrates into a state-of-the-art open-source autopilot for fixedwing UAVs. The next step is to validate the controllers in flight experiments using a physical fixed-wing UAV.
The simulation results show some deterioration of tracking performance when the flow angles are not slowly varying. Therefore, future work could consider different robust and/or adaptive control techniques to better compensate for a wider class of exogenous disturbances, including harsh wind conditions. Also, the coordinated-turn equation, Equation (20), is only an approximation of the relationship between roll angle and turn rate. Future work could seek to relax the assumptions used to derive this relation, e.g., by considering the equations derived in [41], and do a comparative study of different turn coordination methods for fixed-wing UAVs. Another topic for future work is to investigate if a tailor-made guidance scheme can be designed that directly produces a reduced-attitude vector reference.

Appendix B. Proof of Proposition 1
The control law (26) in combination with (18) results in the closed-loop error dynamics Note that the system is time-varying due to the presence of ω d and V a .
To show convergence of e η to zero we make use of the following lemma: Lemma A1. Let x(t) denote a solution to the differential equationẋ = a(t) + b(t) with a(t) a uniformly continuous function. Assume that lim t→∞ x(t) = c and lim t→∞ b(t) = 0, with c a constant value. Then, lim t→∞ẋ (t) = 0 [48,92].
Since e ω converges to zero, we know that b(t) converges to zero. From (17), the derivative of a(t) is given byȧ which is bounded because e η , e ω and ω ⊥ d are bounded. Therefore, a(t) is uniformly continuous, and convergence of e η to zero follows from Lemma A1.
To summarize, all solutions converge to one of the two equilibria given by (η, e ω ) = (±η d , 0). Remark A1. The Lyapunov function (A2) is quadratic, and as we will show, leads to exponential stability, which in general leads to good performance and robustness to perturbations [85]. In [95], it is shown that non-quadratic Lyapunov functions could lead to better performance. Therefore, future work based on the general approach presented in this paper might explore whether different Lyapunov functions than (A2), possibly non-quadratic ones, could lead to better performance.

Appendix B.2. Instability of the Undesired Equilibrium Point
At the equilibrium point (η, e ω ) = (−η d , 0), the value of the Lyapunov function is V = 2k p . To show that this equilibrium is unstable, it suffices to show that for any neighbourhood U around this point, one can find η * , e * ω such that V < 2k p . Since V is non-increasing and all solution converge to either of the two equilibria, any solution starting at (η * , e * ω ) must converge to (η, e ω ) = (η, 0). Consider η * arbitrarily close to −η d , say at an angle away from −η d . Then, ψ(η, η d ) = 1 − cos(π − ) ≈ 2 − and V ≈ k p (2 − ) + e * ω Je * ω /2. This means that, if we choose e * ω small enough, then V < 2k p and we conclude that the equilibrium point is unstable.
Remark A2. This line of reasoning parallels that of Chetaev's Theorem, for which a version for time-invariant systems is given in Theorem 4.3 in [85].

Appendix B.3. Stability of the Desired Equilibrium
We proceed by studying the asymptotic stability of the equilibrium point where η = η d . To this end, consider again the Lyapunov function candidate (A2). From [52], we know that in some neighborhood of (η, e ω ) = (η d , 0), V can be lower and upper bounded by V ≤ 0 together with the positive definite bounds on V makes the equilibrium point (η, e ω ) = (η d , 0) uniformly stable ( [85], Theorem 4.8).
Convergence combined with Lyapunov stability leads to asymptotic stability of the desired equilibrium point (η, e ω ) = (η d , 0). The stable manifold of the unstable equilibrium is less than the dimension of the state space of the system and therefore has measure zero. The region of attraction to the stable equilibrium point excludes this manifold, so we conclude that the desired equilibrium is almost globally asymptotically stable.

Appendix B.4. Exponential Stability
To show exponential stability, let > 0 (arbitarily small) and consider the Lyapunov function (A2) augmented with a cross-term: which is positive definite for small . The time-derivative of V along the closed-loop trajectories satisfiesV =V + e η Jė ω + e ω Jė η .
Combining this with (A3) giveṡ Let x [ e η e ω ] . Then, we can writeV ≤ −x Mx, where the matrix M is given by The matrix M is positive definite if the following inequality is satisfied: Since can be chosen arbitrarily small, this inequality can always be satisfied. Witḣ V negative definite, and with quadratic bounds on V , we conclude that the desired equilibrium is exponentially stable. For estimation of the region of exponential convergence, see [72].

Appendix B.5. Convergence of Angular Velocities
Since e η converges to zero, Π ⊥ η ω ⊥ d → ω ⊥ d . Now, Equation (19) proves our point due to the orthogonality of the two parenthesized terms.
The rest of the proof follows the same arguments as in the proof of Proposition 1.
For some symmetric, positive definite gain matrix K 3 = K 3 > 0, let an augmented Lyapunov function candidate be given by If the update law for∆ is chosen as˙∆ = K 3 z, we are left withV 3 ≤ −κk 1 e η 2 − γ z 2 .
By Barbalat's lemma,V goes to zero asymptotically, and so does e η and z, and therefore also e ω . From (A6), we have Jż = a(t) + b(t), where a(t) −∆ b(t) −k 1 e η − [K 2 − V a D + S(ω d )J]z + S(J(z +ω d ))z Since e η and z converges to zero, we know that b(t) converges to zero. The derivative of a(t) is given byȧ (t) = −˙∆ = −K 3 z which is bounded by the boundedness of z. Therefore, a(t) is uniformly continuous. From Lemma A1, we get that∆ converges to zero as well. The rest of the proof follows closely that of Proposition 1.