An Annular Wing VTOL UAV: Flight Dynamics and Control

A vertical takeoff and landing, unmanned aerial vehicle is presented that features a quadrotor design for propulsion and attitude stabilization, and an annular wing that provides lift in forward flight. The annular wing enhances human safety by enshrouding the propeller blades. Both the annular wing and the propulsion units are fully characterized in forward flight via wind tunnel experiments. An autonomous control system is synthesized that is based on model inversion, and accounts for the aerodynamics of the wing. It also accounts for the dominant aerodynamics of the propellers in forward flight, specifically the thrust and rotor torques when subject to oblique flow conditions. The attitude controller employed is tilt-prioritized, as the aerodynamics are invariant to the twist angle of the vehicle. Outdoor experiments are performed, resulting in accurate tracking of the reference position trajectories at high speeds.


Introduction
Fixed-wing vertical takeoff and landing (VTOL) aircraft (also known as hybrid VTOL) combine both the hovering capability of traditional VTOL aircraft (i.e., helicopters), and efficient high-speed flight of traditional fixed-wing aircraft [1]. Since the 1950s, such vehicles have been of interest primarily for defense applications, where special environmental requirements for takeoff and landing are alleviated, while still allowing for efficient high-speed flight [2]. With the advent of unmanned aerial vehicles (UAVs), there has been renewed interest in hybrid VTOL vehicles, particularly for various civilian applications such as surveillance [3], disaster management [4], and package delivery [5].
Most of the literature has revolved around so-called convertiplanes [6], which use additional mechanisms to transition between hover and fixed-wing mode, such as tilt-rotors/wings. However, such an approach increases the mechanical complexity, thereby increasing cost, maintenance, and risk of failure [7,8]. As a result, there is growing interest in tail-sitter designs, where either control surface or differential thrust based designs are currently under investigation [9]. Although differential thrust tail-sitters typically require more rotors than control surface based tail-sitters, they are more robust due to their larger control authority [9,10]. However, it has been noted that for these types of vehicles, systematic development of flight dynamics models and control system design still need to be explored [4,9,11].
In this paper, a differential thrust tail-sitter is considered based on an annular wing design, as shown in Figure 1. The annular wing acts as a lifting surface in forward flight and encloses the propeller blades to enhance human safety (over a solution without any guard or shroud). Regarding the latter point, there are numerous examples of propeller blades causing serious injuries to the eye or lacerations to the body [12][13][14][15][16][17], and it is claimed in [16][17][18][19]] that a propeller guard or rotor shroud will help mitigate such injuries. Furthermore, the annular wing may lead to reduced noise emission by Drones 2020, 4, 14 2 of 34 the rotors [20], and may provide better cross-wind resistance at hover than planar wing VTOL UAVs. The vehicle builds upon a previous prototype [21], which augmented the Flying Machine Arena (FMA) quadrotor [22] with a foam ring. A similar closed wing design has been adopted by Amazon in their most recent package delivery drone prototype [23].
A literature review is provided in the remainder of this section. The remainder of this paper is organized as follows: • Section 2 presents the specifications for the developed annular wing vehicle. The dynamic model of the vehicle is presented and is used both for simulations and subsequent control system design. Wind tunnel tests are performed to obtain accurate models for both the wing and the propellers under oblique flow conditions. The data presented here can be used as a basis for a further optimized vehicle design, or as a basis for comparison to other VTOL UAVs. • Section 3 presents a novel control system that accounts for the propeller aerodynamics in forward flight, and can be applied to any quadrotor vehicle. The attitude controller presented prioritizes the tilt of the vehicle, as the aerodynamics for the annular wing vehicle are invariant to the vehicle's twist angle (the angle about the vehicle's axis of symmetry), and does not rely on feedback linearizing the angular dynamics. Furthermore, unlike most works in the literature, the control system presented is globally defined and requires no switching depending on the operating region of the VTOL UAV. • Section 4 briefly discusses the software architecture used to achieve high-fidelity software-in-the-loop simulations. • Section 5 presents the outdoor experimental results, where unlike most works, tracking performance on the outermost loop (i.e., position) is shown. A video of the outdoor experiments can be found at: https://youtu.be/ULc63HrgPro.

Literature Review
An early work that investigated transitioning a fixed-wing UAV to a hover position can be found in [24], where the transition control law is based on an approximate dynamic model of the vehicle and a neural network for additional compensation in the dynamic inversion scheme. The transitions were performed not as a consequence of an outer loop position reference, but rather by defining a pitch angle trajectory given directly to an attitude controller with position regulation disabled. A similar approach is also used in [11,[25][26][27][28][29][30][31][32][33]. This approach is related to the work done in [27], where the control system is decomposed into four controllers depending on the operating region of the vehicle: vertical flight, horizontal flight, transition and braking flight. It is noted in [26] that using such a discontinuous control system can cause issues in the form of undesirable transients due to the transitions between each mode. A similar approach is also used in [8,10,28,31,[34][35][36][37][38]. One of the earliest works to forgo the switching approach and use a unified control system is [39], however the experiments for the transition flights were performed in open loop.
Other early tail-sitter concepts can be found in [40,41], which describe the design process and theoretical modeling. Another interesting concept was shown in [42], which aims to reduce wing tip vortices using propellers, and mentions that transition flight maneuvers are difficult to perform as the rotors are at high inflow angles/oblique flow conditions, thereby increasing vibration and other complicated aerodynamic effects. A recent work [10] presents a novel tail-sitter vehicle that contains a single large main rotor as the main propulsion device for both hovering and forward flight. The control approach also uses a switching based approach, where during the transition phase the position is not regulated and a pitch trajectory is tracked, resulting in altitude deviations as large as 20 m. The authors state that this will be improved upon in a future work.
With regards to conventional fixed-wing tail-sitters that are based on a control surface design, ref. [25] considers a transition controller that utilizes a pre-planned reference attitude trajectory, which is given to an attitude controller that utilizes gain scheduling. In [43], a hovering controller was presented, and is one of the earliest works to decompose the attitude error into tilt and twist components. In [44], Drones 2020, 4, 14 3 of 34 this tail-sitter was made to hover and move in an indoor environment, while avoiding obstacles detected with an ultrasonic range sensor. Other related works include [45], where the transition is performed by using a gain scheduling controller to improve upon the traditional "stall-and-tumble" maneuver to achieve the transition. In [46], a fixed-wing vehicle that uses contra-rotating propellers is presented, and a nonlinear attitude control based on Euler angles is derived and experimentally shown to work under cruise conditions. One of the earliest work on VTOL control that takes into account the propeller aerodynamics under oblique flow conditions is [29], however the model was made for a fixed angle of attack. In [47], a unified control system based on an adaptive control to account for the aerodynamics of the wing is presented and experimentally shown to perform well indoors. In [48], an incremental control structure is used that relies on the control derivatives of the aircraft at certain operating points, as opposed to knowledge of a detailed aerodynamic model, with outdoor experiments showing its effectiveness.
In the context of quadrotor based tail-sitters, early simulation work can be found in [30], where the pitch is commanded directly to perform the transition. In [18], various VTOL design concepts are presented, and it is mentioned that some type of rotor shroud is advantageous to minimize damage after a casual impact, or for human safety when flying at low altitudes. In [31], a switching control system is used on a quadrotor tail-sitter that additionally utilizes control surfaces, where the transition is performed by adjusting the pitch angle reference directly. In [32], it is claimed that fixed-wing tail-sitters have the disadvantage of relying on control surfaces and result in less accurate flight performance. Two different wing arrangements are also tested on a quadrotor. The transition is performed by providing a pre-planned attitude trajectory to an attitude controller that uses a tilt-twist approach similar to [43]. A similar approach is investigated in [35] and in [34], where a switch is implemented for forward flight mode. In a later work [33], it is shown that the standard transition control approach used in their previous work [32] results in significant altitude loss of the vehicle. This was improved upon by computing an optimal trajectory for the transition maneuver and storing it onboard. An interesting design that combines both a quadrotor for attitude stabilization and contra-rotating propellers as the main propulsion source can be found in [37]. In [8], a package delivery drone is designed and a switching based control system is synthesized to perform the various modes of flight, although fully autonomous flight was not implemented. In [11], a novel bi-plane quadrotor is introduced with corresponding theoretical model and controller synthesis. This idea was extended in [49] by incorporating a passive rotary joint in between the two wings, allowing the rotors to swivel and thereby increase yaw torque authority. The recent work [20] improves upon the switching control system approach by designing a unified controller, and is another early work to account for the propeller aerodynamics in forward flight, specifically in axial flow. In contrast to this paper, their approach relies on training based on simulation data. A unified control system approach is also used in [50], capable of continuous transition between hover and fixed-wing flight. Their approach is based on a nonlinear optimization and is shown to perform well in simulations.
Closed wings have also been researched by the aerodynamics community since the 1920s [51], and have the advantage of eliminating wing tip vortices and reducing lift-induced drag [52,53]. There is renewed interest, e.g., from NASA and Lockheed Martin, on using such wing designs for new aircraft [54]. In [55], it is noted that studies on annular wings are lacking in the literature, and a wind tunnel characterization was undertaken for two different aspect ratio wings for small angles of attack. In [56], wind tunnel characterization was also performed for an annular wing, with additional characterization of a ducted propeller mounted in both pusher and tractor configurations.

Overview
Two vehicles were developed, both shown in Figure 1. The annular wing of the white vehicle is based on a flat plate airfoil. The annular wing of the blue vehicle is based on the E169-il airfoil [57], and was manufactured in collaboration with the Product Development Group at ETH Zurich [58]. The dimensions of the vehicle can be found in Figure 2. The dimensions were determined by aiming at a target operating condition as follows: A first principles based model of the wing lift is derived based on integrating the airfoil aerodynamic coefficients around the wing, as shown in Appendix A. The wing diameter D was fixed to 70 cm as a design parameter. The wing chord c was determined to be that which achieves 7 N (about the target weight of the vehicle) of lift at the target operating condition of 10 m/s and 12 degrees angle of attack (approximate location of the optimal lift-to-drag ratio using the theoretical model in Appendix A, and in agreement with [53]), using the aforementioned theoretical model. The vehicle's center of mass lies on the quarter chord point of the wing, and is the point at which the carbon tubes that connect to the wing intersect.   Table 1 details additional specifications, such as the sensors used and measured power consumption of the vehicles. The blue vehicle (0.75 kg) is more efficient than the white one (0.71 kg), as expected due to the more optimized airfoil despite being heavier; it is able to sustain level flight at 10 m/s at a specific power [2] of 100 W/kg, whereas the white vehicle consumes a slightly larger value of 139 W/kg. For reference, the FMA quadrotor [22], which is of similar size and uses similar flexible propellers, consumes 275 W/kg at 10 m/s (determined from the tethered flight experiments performed in [59]).
These numbers can be improved with further design enhancements. The blue vehicle's wing is very sensitive to heat: the experiments were performed in cold temperatures, around 0 • C, and this would cause the rib structure to contract and the foil to expand, thereby inducing wrinkles in the surface of the airfoil. Both wings are quite flexible, and vibration due to aeroelastic effects was noticed in high-speed forward flight, see Section 2.2.1. The fuselage and the ribs for the blue wing were also printed using standard 3D printed plastic (i.e., PA 2200). Thus, lighter, rigid, and less heat sensitive materials, such as carbon fiber, should be used for a more efficient vehicle design. The tail of the fuselage should also be made less flat to reduce pressure drag and thereby improve aerodynamic efficiency, however the flat tail was used to yield a stable landing surface as the wings were quite fragile. Stiff, high-pitch propellers should also be used to improve the efficiency at high forward flight speeds.
The flight computer and Global Positioning System (GPS) printed circuit boards were custom built, based on the STM32 F4 microcontroller and the Ublox M8P module, respectively.

Modeling
In the following, the vector p B/I represents the position of the body frame B with respect to the inertial frame I, where the orientation of a frame is defined by the three orthonormal unit vectors {i, j, k}. A superscript will be used to realize this vector with respect to a particular frame, e.g., p I B/I ∈ R 3 represents the vector expressed in frame I. If a dot over a vector is present, the superscript/frame assignment is performed first, then the time derivative is taken. The same is true for the velocity vector v B/I , and the angular velocity ω B/I . An index subscript, e.g., p I B/I l , represents the lth element, l ∈ {0, 1, 2}, of that vector. The coordinate systems used are shown in Figure 3, where B refers to the body-fixed frame attached to the center of mass of the vehicle, R refers to a similar frame, but is rotated about the k axis and will be used in modeling the rotor force and moments, A denotes the air frame fixed to the incoming airflow to the vehicle, and P i , i = 0, . . . , 3 represents the propeller frame that is fixed to the ith propeller. Furthermore, the forward frame rotation from frame I to frame B is denoted The wind velocity vector, or the velocity of the air with respect to the inertial frame, is denoted as v A/I . Thus the velocity of the air with respect to the body is given by Let n(·) denote the normalization operator, that is n(x) = x/ x for some vector x. The forward frame rotation from the body frame to the air frame is given by where V ∞ is the air speed with respect to the body. The other common frame rotation encountered is between the body and the rotor frame: Frame A is fixed to the incoming air which forms an angle of attack α with the vehicle's k-axis. The rotor frame R has its k axis aligned with that of the body frame, but is rotated such that the incoming air lies solely on its i, k-plane.
The vehicle dynamics are given bẏ where ω B B/I × is the skew symmetric matrix form of ω B B/I , Ω i is the rotation rate of the ith propeller, ω B P i /B = (0, 0, Ω i ), g I = (0, 0, 9.81)m/s 2 is the acceleration due to gravity, v A/P i = v A/B − ω B/I × p P i /B is the velocity of the air with respect to the ith rotor, p P i /B is the position vector of the ith rotor with Drones 2020, 4, 14 7 of 34 respect to the center of mass, and J B B , J B P are the inertia tensors for the entire vehicle and propeller, respectively, represented in the body frame. (9) represents the motor dynamics, where u i is the input voltage to the ith motor, and R m , K m are the motor resistance and torque constants, respectively. Henceforth, the approximation will be used, under the reasonable assumption that the body rates and the rotor arm lengths are small. The maps F B W (·), M B W (·) represent aerodynamic wing loads, and F B P (·), M B P (·) represent aerodynamic propeller loads, and will be discussed in the subsequent subsections.

Wing
By construction of the air frame A as shown in Figure 3, the annular wing's aerodynamic force and moment can be written in the following way: where F WD (·), F W L (·) and M WP (·) map the free-stream velocity expressed in the body frame to the aerodynamic drag, lift, and pitching moment of the wing, respectively, and can be expressed in the following form [61]: where is the angle of attack (see Figure 3), ρ is the air density, and S := Dc is the reference planform area [61]. To account for some of the aeroelastic vibration effects observed during the wind tunnel tests, the coefficients C W L (·), C WD (·), and C WP (·) will be written as the sum of a function dependent solely on the angle of attack and the output of a colored noise process as follows (minding the slight abuse of notation by including the time t as an additional argument to the coefficients): where η C W L (t) is the output of a first-order stochastic procesṡ where η u (t) is drawn from a to-be determined distribution and τ is a constant to be determined. A similar approach is used for the remaining coefficients C WD (·) and C WP (·). The motivation for including this is to add an additional level of realism to the simulation, and to explore the robustness of the control system. The aerodynamic coefficient maps C W L (·), C WD (·), and C WP (·) are determined using a wind tunnel using the experimental set-up shown in Figure 4. Note that the maps will also incorporate the aerodynamic loads of the vehicle's body, in addition to the wing. The wind speeds V ∞ were set to 3, 5, 6.5, 8, 10, 12 and 15 m/s. For each fixed wind speed, the angle of attack α was swept at a rate of 0.2 deg/s from −10 to 90 deg (hardware limitations prevented sweeping to higher angles of attack). Note that aerodynamic hysteresis effects [62,63] are ignored. The load data was collected at a rate of 1 kHz using a six degree of freedom load cell sensor, and were translated from the point of measurement to the center of mass of the vehicle/quarter chord point of the wing. The load cell used is an ATI Mini40 with the SI-20-1 calibration, which has a negligible noise floor of at most 0.02 N and 0.0005 Nm, peak-to-peak. The load cell communicates the load data to the host computer via a direct Ethernet connection. Both the wind speed V ∞ , which is measured using a pitot tube and pressure sensor, and servo angle α, are sampled via a LabJack U12 data acquisition board and sent to the host PC over a USB connection. Note that the measurements were taken using a different fuselage tail, so the results are expected to differ from the vehicle used in the actual flight tests.
The measured aerodynamic coefficients are shown in Figure 5 as function of α. From the figure, it can be seen that the stall behavior is less pronounced than for planar airfoils. This can partly be explained using the theoretical model in Appendix A. Also, the drag curves do not exhibit the usual curves of planar airfoils, which would have maximum drag at α = π/2. In the case of the annular wing, the drag at α = π/2 is small because of the curved surface that the incoming air can glide around, similar to a cylinder in cross-flow [64]. These behaviors were also observed in the early wind tunnel testing done in [53]. This could suggest that annular wing vehicles have better cross-wind resistance at hover than more traditional planar wing VTOLs [65]. The mean lift and drag curves are parametrized by a piece-wise linear fit, each assumed to be an odd and even function [66], respectively, where for positive angles of attack have the form: where the constants c l,α , c l,0 , α l , c d,α , c d,0 , and α d are to be determined. For the lift coefficient, this is done by solving the following optimization s.t. c l,α,0 α l,0 =c l,α,1 α l,0 +c l,0,1 c l,α,1 α l,1 +c l,0,1 =c l,α,2 α l,1 +c l,0,2 where S(·) is the Huber loss function for robustness against outliers [67] (alternatively, the sum of squares function or sum of absolute values can be used), M is the set of all measurements, and the superscript "meas" is used to denote measured quantities. The constraints on the inner optimization are to ensure C W L (·) is continuous in its argument α. The inner optimization is convex and can be computed efficiently using CVXPY [68,69], while the outer optimization is not convex and is computed using the minimize routine in Scipy's optimize package [70]. A similar procedure is performed for fitting the drag coefficient. For the pitching moment coefficient, a simple sinusoid term suffices: which can be determined easily from the measured data using a similar convex optimization routine. The next step is to characterize the colored noise process (18). This is done by first discretizing (18) using an Euler approximation to yield where T s = 0.001s is the sampling period and a = 1 − T s /τ. Then, the measured noise sequence is defined as Using this measured noise sequence, the parameter a and the variance of η u [k] can be computed using the autocorrelation (Yule-Walker) method as described in [71]. The distribution of η u [k] can be determined by plotting a histogram of as shown in Figure 6, from which it can be observed that the distribution can be approximated by a Gaussian. With the distribution of η u [k] determined, and the parameter a determined earlier, the stochastic process (18) can be simulated and is shown in Figure 7. The simulated noisy sequence resembles the measured η meas C W L [k], though a higher fidelity model can be realized by using a higher order autoregressive model. The fitted aerodynamic parameters can be found in Table 2.  (25), where a is determined using the Youle-Walker method [71] for the first-order autoregressive model (23).

Propulsion
The rotor force and moment maps are modeled using the parametric models provided in [72]: where F PH (·) is the propeller H-force, F PT (·) is the propeller thrust, F PS (·) is the lateral/side force, and M PR (·), M PP (·), and M PQ (·) are the rolling, pitching and rotor torque, respectively. It is shown in (Section V, [72]), that these forces and moments can be written as are the rotor climb and advance ratios for the ith rotor, respectively, R is the rotor radius, and the twelve c P constants are to be determined. Note, however, as mentioned in (Section V, [72]), the pitching moment model (32) may be inaccurate. The constants are again determined by measuring the loads in a wind tunnel. The experimental set-up, shown in Figure 8, is the same as before with the same sweeping performed on V ∞ and α, but with two propellers mounted onto the vehicle. Two propellers were used because if all four were placed, and if all were sent the same rotation rate commands, the rotor torque M PQ and rolling moments M PR would cancel each other out. Furthermore, if only one was used, then this might cause a biased reading at that location. The propellers are commanded to sweep between 150 and 600 rad/s, sinusoidally with a period of 60s. The measured rotation rates Ω i are computed by the onboard STM32F4 microcontroller (see Section 3.1), which is then communicated to the host computer via the Laird RM024 radio.
The mean wing force and moment are subtracted from the measured loads using the identified model (19)- (22), to yield the measured rotor forces and moments from which the constants in the propeller model (28)-(33) are identified using a convex Huber regression [67]. The results are shown in Figure 9.
The thrust and rotor torque models fit very well, outperforming the standard hover model typically used in the drone literature [72], which is also shown for reference in Figure 9. The propeller H-force is very small for small angles of attack. The pitching moment model is inaccurate, however it can be shown that this is about an order of magnitude smaller than the pitching moment of the wing, and thus can be neglected. The slightly non-zero side force, and slightly inaccurate rolling and H-force models, may be the result of an asymmetrical induced flow observed by the rotor disk [73], resonance from the aeroelastic and vibration effects, or due to a slight misalignment of the load cell/vehicle. The numerical values of the aerodynamic coefficients can be found in Table 3.

Two Propellers
12V Power cables

Overview
The overall control architecture is based on a cascaded scheme as shown in Figure 10. A state estimator is used to fuse the sensor readings and is based on an extended Kalman filter (EKF) using an attitude reset scheme as described in [74], fusing the inertial measurement unit (IMU), magnetometer, pressure sensor and GPS. The state estimator provides an estimate of the body frame, denotedB, in the form of a position estimate p IB /I , velocity estimate v IB /I , and attitude estimate R IB . The estimate of the angular velocity, ωB B/I , is read directly from the inertial measurement unit.
The estimate of the propeller rotation rates,Ω i , is measured using a timer on the STM32F4 and a pulse representing the zero-phase crossings that is provided by the SimonK firmware running on the ESCs [75,76]. GPS latency is also taken into account, by capturing the EKF estimate at the time of a measurement epoch to generate the innovation.

Position Controller
Consider the position trajectory of the reference frame, where k p,p , k p,i , and k p,d are scalars representing the PID gains. This describes the acceleration of the frame B cmd . The orientation of B cmd will be determined in the Outer Control Allocation block, described in the following section.

Outer Control Allocation
The commanded acceleration a I B cmd /I is then fed to a control allocation scheme which determines the corresponding vehicle attitude, R I B cmd , and total collective vehicle thrust, T cmd , to achieve the commanded acceleration. To do so, the approximation F PH ≈ 0 is made, as the resulting tracking performance is adequate (see Section 5), although taking this into account may be subject for future work. Furthermore, as shown in Section 2.2.2, although the propeller H-force can become significant at high flight speeds at large angles of attack, this is not an operating condition in which the vehicle will spend much time.
Additionally, the approximation v A/I ≈ 0 is made as no onboard wind speed estimator is available, and thus the feedback loop (36) will be used to compensate for this. Then, the control allocation problem can be formulated as follows: where F W (·) = F W (·) from (11), but with the aeroelastic noise component set to zero. This is done by first constructing an intermediate frame A cmd (see Figure 11) as follows: Now the commanded force m a B cmd /I + g must lie on the 2-dimension i A cmd , k A cmd plane. Let The control allocation problem can thus be simplified to the following 2-dimensional problem: where ·, · is the inner product operator. A solution to the above always exists. This can be verified visually in Figure 12, where the range of (41) is plotted over all α cmd and T cmd . It can be observed that the range space of (41) is simply-connected (and also convex), and thus for sufficiently small commanded forces f x = m a I B cmd /I + g I , i I A cmd , f z = m a I B cmd /I + g I , k I A cmd , there will always be a solution T cmd , α cmd that solves the allocation problem (42). Thus Newton's method [67] is used to solve (42). Although C W L (·) and C WD (·) are continuous but not differentiable at the knot points, simple heuristics or the Secant method [4] can be used to circumvent this at the knots. Alternatively, the piece-wise linear curves can be made smooth by using a suitable sigmoid function at each knot.
With T cmd , α cmd determined, the rotation R I B cmd in (37)  Note that this choice is arbitrary due to symmetry of the vehicle, and would not be possible for hybrid VTOL UAVs based on a planar airfoil, where coordinated turn maneuvers are required [47].

Feedforward Body Rates
For accurate tracking, the reference body rates need to be sent to the attitude controller, as shown in Figure 10. This is done by repeating the control allocation steps of Section 3. 3

Attitude Controller
The attitude error is the difference between the estimated body frameB and the commanded body frame B cmd , represented as a rotation matrix as follows: where exp(·) maps a rotation vector to the corresponding rotation matrix [60], R tilt B represents the tilt error, and R B cmd tilt represents a pure twist error. This decomposition can be done efficiently using unit-quaternions, see [78] or Appendix B. This is done because the twist error of the vehicle has no consequence on the tracking of the commanded acceleration, due to symmetry of the vehicle. Thus, the decomposition will be used to prioritize the tilt errors over twist errors, as follows: Let q be the unit-quaternion representation of the rotation matrix R, and q its vector part. Then, the control law used for the attitude subsystem is given by (48) where k p,tilt , k p,ψ are positive scalars representing the proportional gains for the tilt and twist errors, respectively, and K D is a positive definite matrix representing the derivative gains. The stability proof can be found in Appendix C. Another approach would be to replace the full attitude error term q B cmd B in (48) with the pure twist term q B cmd tilt , although the corresponding stability proof has not been found and is a suitable direction for future work. This approach differs from the tilt prioritization controller of [78,79], which includes (additional) feedback terms to linearize the angular dynamics.

Inner Control Allocation
The commanded body torques τ cmd and thrust T cmd are then fed to a control allocation block which determines the corresponding propeller rates Ω cmd . Again, assuming that v A/I = 0 as discussed in Section 3.3, this can be written formally as follows: find Ω cmd s.t.
where the maps F PT (·), M B P (·), F B P (·) are given in (26) (12), but with the aeroelastic noise component set to zero. To solve (49), the approximation of F PH ≈ 0 is made, as discussed in Section 3.3, and additionally are made, although taking this into account is a suitable direction for future work. From a modeling perspective, the approximation M PP ≈ 0 is reasonable since as mentioned in Section 2.2, the pitching moment of the propeller is dominated by the wing. Next, the approximation M PR ≈ 0 is reasonable if the propeller rotation rate of a clockwise rotating propeller is roughly similar to that of a counter clockwise rotating propeller, such that the rolling moments approximately cancel each other out. From an experimental perspective, these approximations provide adequate tracking even at the high speed regime (see Section 5).
With these simplifications, (49) can be written equivalently as: find where (Ω cmd ) 2 is performed element-wise, where l arm is the distance between the vehicle center and a propeller's center, are defined for brevity, and where Ω H is the standard quadrotor hover solution: Note that the C and D matrices are defined if the body frame B is as shown in Figure 3, modulo the ordering of the rotors. If they are instead aligned with the rotor arms, as is typically done in the quadrotor literature, then: One way to solve (52) is to use Newton's method, as done in the outer loop control allocation block. However, this results in excessive computation, requiring the inverse of a 4 × 4 Hessian matrix. As an alternative, let where the square-root is applied element-wise. The goal is to find the fixed point Ω cmd := Ω * such that Ω * = f (Ω * ). Note that where the norm applied to the Jacobian is the induced matrix 2-norm, and λ * c i . f (·) will be a local contraction mapping if (68) is strictly less than 1 [80]. This will always be the case for a sufficiently small operating region in the climb ratio λ c . For example, using the parameter values from Table 3, this means that for the annular wing vehicle, the operating λ c ratio must be less than 0.68. This is a very extreme operating condition and will not be encountered in practice, as a climb ratio of 0.2 corresponds roughly to when UAV propellers produce zero thrust [72]. Since f (·) is a local contraction about the solution, then by Banach's Fixed Point Theorem [80], Algorithm 1 can be used to find the fixed point Ω * .
Algorithm 1: Iterative algorithm to solve the inner control allocation problem (52) input: desired body torques τ cmd and thrust T cmd , reference velocity expressed in the inertial frame v I B ref /I , the desired body frame expressed as a rotation matrix R I B cmd (from which α cmd can be computed), and initial guess Ω (0) . output: the propeller rates Ω cmd .
Ω cmd ← Ω (0) ; for iteration = 1,2,... do Ω cmd ← f (Ω cmd ); end In practice, Algorithm 1 works for λ c larger than 0.68, as Banach's Fixed Point Theorem is only a sufficiency condition, although such a scenario is not so interesting as this would not be encountered in practice. In practice, Algorithm 1 is initialized with the previous time-step's solution.
See Figure 13 for a simulation of Algorithm 1, where the commanded operating condition is the target operating speed of v I B ref /I = 10m/s. The hover solution Ω H is used to initialize Algorithm 1 in this case. This figure demonstrates two things: at this high speed, the hover solution is off by a large amount, as the propellers need to rotate faster to produce the same thrust (under predominantly axial flow), as evident during the modeling done in Section 2.2. The second is that Algorithm 1 converges very quickly, often only requiring a couple of iterations in practice. Figure 14 shows the closed loop tracking results when using the proposed scheme presented herein, and if the hover solution Ω H were applied instead (both schema use the same feedback control gains as used in the experiments in Section 5). Here it can be seen that the proposed scheme tracks the desired high-speed maneuver accurately, whereas the hover solution fails to even hit the target reference speed of 10m/s. Of course, the standard hover approach can be improved with further tuning (i.e., increasing) the feedback gains, however doing so will have undesirable consequences, such as causing the system to be more sensitive to disturbances and noise.  Figure 14. Simulated closed loop response using two different inner control allocation schemes: the hover approach, where Ω cmd = Ω H as defined in (61), and the proposed approach using Algorithm 1.
The feedback gains used are the same as those used in the experiments in Section 5. The difference in the responses can be exaggerated by using a less aggressive outer loop controller.

Motor Controller
The commanded propeller rates Ω cmd are then fed to the motor controller, which uses a combination of a PI controller for feedback compensation and a feedforward term as follows: where k p,Ω , k i,Ω are the PI gains. The feedforward term is again necessary for accurate tracking and compensation of the propeller torque that is encountered in both hover and forward flight.

Trajectories
It can be shown that the current draw of the motors is a function of the crackle (5th derivative) of the vehicle's position trajectory. Thus, to minimize current spikes, trajectories that have bounded crackle (see Figure 15) are used to generate the reference position trajectory and its derivatives. This is done by using the harmonic jerk model of [81]: where j max , T J are design parameters, and t i is the initial time instant of the corresponding segment. These parameters are selected so that a desired position set-point is achieved in each axis.

Software Implementation
An important aspect of the control system is the software architecture used to implement it. The approach here uses the NuttX real-time operating system (RTOS), which is a light-weight, POSIX compliant RTOS [82]. The POSIX compliance is the critical ingredient used to achieve high fidelity software-in-the-loop simulations, as the application software written for the embedded platform can also be compiled to run on a host's POSIX compliant personal computer (PC) (such as a distribution based on GNU/Linux). This allows for certain bugs such as stack overflows, race-conditions arising from context switches and interprocess communication (IPC) to be simulated [83], and thus can be caught in advance instead of having to debug software directly on the embedded platform. This is unlike the approach of [22,84], where only the control algorithms are simulated as common library routines, or a Simulink (or equivalent) approach [7,20,37,[85][86][87][88][89][90], where there is a clear separation or distinction between simulated code and code implemented on the embedded platform.
The software architecture employed is shown in Figure 16. The onboard task can run either on the embedded platform, or on the host PC for simulation. The onboard task is oblivious to whether it is being run as a simulation or not. It reads in sensor values using dedicated polling threads through the virtual file interface via the device drivers [83], and then sends the motor voltages u cmd via the same interface. The device drivers provide the needed abstraction, where they are implemented to communicate with the real sensors and actuators on NuttX. For simulation, they are implemented on a GNU/Linux machine via kernel modules [91], which relay information to and from the simulator process.
All of the control laws discussed in Section 3 are implemented in the fast thread running within the onboard main task. It runs at a fixed frequency of 500Hz that is controlled by a dedicated timer driver /dev/timer0 via signalling [83]. The slow thread communicates with the radio device for sending telemetry and receiving high level commands from the user (such as takeoff, land, or fly a certain trajectory). It is also responsible for monitoring for any sensor faults and bringing the vehicle to an emergency land state. The offboard task runs only on the host PC and runs regardless of whether a simulation or real experiment is being run. The offboard task communicates with the host radio device for sending high level commands to the vehicle and receiving telemetry. It also provides the user with a real time visualization of the received telemetry data. The simulator task runs on the host PC, using the system model described in Section 2 to drive the ordinary differential equation (ODE) solver [92]. Drivers (Linux) Figure 16. Software architecture used to implement the control system. Multiple threads are run within each process/task, as shown, and communicate amongst each other using standard IPC mechanisms [83]. The onboard main task runs either on the GNU/Linux machine for simulation, or on the embedded NuttX platform.

Experimental Results
The experiments were performed outdoors; a video can be found at https://youtu.be/ULc63HrgPro. The experiments were performed in conditions with wind speeds of 2m/s and occasional wind gusts of up to 5m/s, according to local weather forecasts provided by [93].
Two experiments are performed. The first is a straight line trajectory, shown in Figures 17 and 18 i.e., the arccos of the last entry of the rotation matrix, is shown in Figure 19. The tilt angle varies between 0 and 80 degrees during the straight line maneuver. The forward and reverse pass differ due to external wind disturbances and possibly by asymmetries due to manufacturing imperfections. The second trajectory is shown in Figures 20-22, where the vehicle is commanded to travel in a circle of radius 10m at a speed of 10m/s. The results again show good tracking of the reference trajectory throughout the high-speed maneuver. Figure 22 shows that the vehicle's k-axis is not tangent with the circle, as the annular wing is not only providing lift to counteract gravity, but also generating lift in the plane of the circle to provide the required centripetal acceleration, mixed of course with the rotor thrust. This mixing is done automatically by the outer control allocation block as discussed in Section 3.

Conclusions
A quadrotor with an annular wing has been presented that enhances human safety and can provide efficient forward flight. The annular wing and propellers have been characterized through wind tunnel experiments. The designed model-based control system accounts for the dominant aerodynamics of the propellers under general oblique flow conditions using an iterative algorithm, and utilizes a tilt-prioritized attitude controller that takes advantage of the symmetry of the vehicle.
The control system has been demonstrated through outdoor experiments and shown to provide accurate transition maneuvers and high-speed trajectory tracking.
The iterative control allocation algorithm presented in Section 3.6 can be applied to any quadrotor vehicle, and is a suitable direction for future work for improving tracking performance at high speeds. A more optimized vehicle design, using lighter and more rigid materials, is a suitable direction for future work. Other directions include taking into account the propeller H-force and incorporating an onboard wind speed estimate in the control allocation blocks, and seeking the attitude control stability proof when a pure twist term is used in the control law (48). Other directions include performing a rigorous comparison on cross-wind resistance and noise emission of the proposed annular wing vehicle with conventional planar wing tail-sitters.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Theoretical Lift Model
Consider an infinitesimal wing segment and its associate frame dw, as shown in Figure A1. The forward frame rotation from the body frame B to the frame of the infinitesimal wing is given by Also, for the purpose of the analysis, we can assume that n k B Then the velocity of the air expressed in the infinitesimal wing's frame is The angle of attack of the infinitesimal wing element is then where C dwL (·) is the lift coefficient polar of the E169-il airfoil provided by [57] using the following parameterization: C dwL (α) = c dwL,α α(1 − σ(α)) + 2σ(α) sin(α) cos(α) (A6) with c dwL,α = 4. The form for the roll off in the lift coefficient at high angles of attack is provided by [94], where σ(·) is a sigmoid function and is given by where M, α 0 are parameters representing the sharpness and location of the airfoil stall, respectively. Figure A2 shows the resulting lift curve C W L (·), with a comparison to the planar lift curve C dwL (·).
The annular wing has a larger coefficient of lift at all angles of attack than the planar airfoil (note both use the same reference planform area of S = Dc), which is expected. Additionally, the stall behavior is slightly washed out by the annular wing. A similar exercise can be performed for the drag coefficient (see e.g., Equation (4), [59]).

B
Let R(q) denote the rotation matrix mapping from a unit-quaternion q = (q 0 , q 1 , q 2 , q 3 ) (with q 0 the scalar part of the unit-quaternion, and the vector part is denoted q = (q 1 , q 2 , q 3 )): 2q 1 q 2 − 2q 0 q 3 2q 1 q 3 + 2q 0 q 2 2q 1 q 2 + q 0 q 3 q 2 0 − q 2 1 + q 2 2 − q 2 3 2q 2 q 3 − 2q 0 q 1 2q 1 q 3 − 2q 0 q 2 2q 2 q 3 + 2q 0 q 1 which is defined so that for two unit-quaternions q and p, R(q p) = R(q)R(p) where is the quaternion multiplication operator defined using Hamilton's convention. Note the definition (A8) is different to that used in [60], which sparked the NASA JPL convention for the quaternion multiplication. A related discussion as to why the approach used here is preferable can be found in the recent work [95].
Solving the above yields [78]: The kinematics of q is given by: Assume that without loss of generality, the inertia tensor J is a diagonal matrix. Then, the following identity is verifiable [97]: which implies that ω B B/B ref is bounded. It can be shown that the second time derivative,V, is also bounded, which by Barbalat's Lemma [98] implies thatV goes to 0 as t goes to ∞. Equivalently, by (A24) this means that ω B B/B ref goes to 0. This also implies that the attitude error must go to zero, evident from substituting the control law (A18) into the error dynamics (A16). Of course this is true for all initial conditions except when q tilt B,0 = 0, which results in a singularity in the control law as the vector part q tilt B becomes ill-defined (see (A12)).