Fault-Tolerant Control of a Variable-Pitch Quadrotor under Actuator Loss of Effectiveness and Wind Perturbations

The actuator fault-tolerant control problem for a variable-pitch quadrotor is addressed under uncertain conditions. Following a model-based approach, the plant nonlinear dynamics are faced with a disturbance observer-based control and a sequential quadratic programming control allocation, where only kinematic data of the onboard inertial measurement unit are required for the fault-tolerant control, i.e., it does not require the measurement of the motor speed nor the current drawn by the actuators. In the case of almost horizontal wind, a single observer handles both faults and the external disturbance. The estimation of the wind is fed forward by the controller, while the actuator fault estimation is exploited in the control allocation layer, which copes with the variable-pitch nonlinear dynamics, thrust saturation, and rate limits. Numerical simulations in the presence of measurement noise show the capability of the scheme to handle multiple actuator faults in a windy environment.


Introduction
Multirotor helicopters are a class of Unmanned Aerial Vehicle (UAV) that are employed in a wide range of applications, resulting in one of the most appealing classes of vehicles for both research and industry. The standard quadrotor represents the most common configuration, which is equipped with four fixed-pitch blades connected to independent motors. Variable-Pitch Quadrotors (VPQs) can instead vary both the rotation speed of each motor and the blade pitch of each propeller, theoretically doubling the degrees of freedom of each actuator. VPQs are also known as collective pitch quadcopters or heliquads.
This results in many advantages, from which it is possible to list higher thrust rate of change, reverse thrust, reverse flight capabilities, and scaling well with size [1][2][3]. Moreover, the popular ArduPilot control suite supports VPQs using the Heliquad configuration. Despite this, VPQs still represent a market niche, where few commercial devices (e.g., the Assault Reaper 500, the Stingray 500, and the WLtoys Assassin V383, for recreational use) are available. In most of the scientific literature, however, custom prototypes (e.g., [4][5][6]) are developed to test algorithms and control strategies. Fixed-pitch quadrotors and VPQs have similar dynamics, while relevant differences arise in the control effectiveness matrix. This allows for identical control laws for both configurations, but VPQs require a more complex control allocation algorithm. The additional degrees of freedom, provided by the pitch angles, can be exploited to satisfy additional constraints, such as minimizing energy consumption, handling the presence of faults and failures (as defined in [7]), etc.
Specifically, the control allocation algorithm, which manages the redundancy while respecting constraints, strictly depends on two key factors: the thrust model and the choice of the control inputs. on the choice of the control inputs (i.e., centralized or independent propeller speed). On the other hand, the control allocation algorithm must exploit redundancy to achieve fault tolerance, so the choice of the control inputs is relevant; however, the proposed control allocation algorithm is designed such that it suffices to enable an equality constraint to deal with a centrally powered VPQ (instead of a conventional one with four separate motors).
The paper is structured as follows. In Section 2, the mathematical model of the vehicle is detailed, including the characterization of variable-pitch propellers and the model of actuator faults to be faced. Section 3 is devoted to the design of a Non-linear Disturbance Observer (NDO) to estimate faults and disturbances. In Section 4, a double-loop controller to track position and orientation, despite the presence of external disturbances, is introduced. In Section 5, a control allocation strategy that redistributes the control effort by taking into account the presence of actuator faults is detailed. Simulation results to show the effectiveness of the overall fault-tolerant control strategy are proposed in Section 6. Finally, conclusions and future works end the paper.

Mathematical Model
The overall kinematics and dynamics of a generic multirotor, thus including variablepitch quadrotors, can be approximated as that of a unique rigid body with six degrees of freedom [21]. Let us consider an earth-fixed frame R E = (O E , x E , y E , z E ), which is assumed to be inertial, and a body-fixed frame R B = (O B , x B , y B , z B ), which is placed in the center of mass of the vehicle. For simplicity, we assume the center of mass coincides with the geometric center of the quadrotor, as depicted in Figure 1. From now on, the standard basis vectors of R 3 are denoted with e 1 , e 2 , e 3 . Let us denote with p F = col(x F , y F , z F ) the position of the center of mass with respect to R E , and with ω = col(p, q, r) the angular velocity of the vehicle. Neglecting the second-order effects (i.e., rotor dynamics, blade flapping, gyroscopic and inertial effects due to the rotors), the quadrotor kinematics and dynamics can be expressed as [21] where m is the total system mass, J = diag(J x , J y , J z ) is the tensor of inertia along the x B , y B , and z B axes, k t and k r are the linear and angular friction coefficients, η = col(ϕ, θ, ψ) is the vector of the attitude angles (roll, pitch, and yaw, respectively) that characterize the rotation from R B to R E .
where cos(·) = c (·) and sin(·) = s (·) are introduced for the sake of brevity, F B m and τ B m are the overall force and torque due to the motors (decomposed in R B ), F E w is the force due to the wind (decomposed in R E ), and T(η) is the kinematic coordinate transformation related to the adopted roll-pitch-yaw rotation, that is [22]: Let us denote with f i ∈ R the actual lift thrust produced by the i-th actuator. In VPQs, as in conventional quadrotors, all the actuators are oriented upwards and they produce collinear lift forces only, as depicted in Figure 2. Denoting with u z the overall lift force, given by In other words, the overall lift thrust is along z B only, because none of the actuators can generate a horizontal force in R B . Then, substituting (4) in the first row of (1), we obtain the following model: The four-dimensional vector col(u z , τ B m ) = col(u z , τ p , τ q , τ r ) is the so called virtual input vector.

Input Mapping
Consider the i-th motor. We denote with f i ∈ R the lift thrust generated by the propeller attached to the motor, with d i ∈ R the generated drag, and with Ω i ∈ R the angular rate of the propeller, taken to be positive when the rotation of the motor matches that of the propeller (clockwise or counter-clockwise). Let us consider the case of rectangular, untwisted blades (which are suitable for variable-pitch drones) with uniform inflow. Let us assume the blade pitch is limited in order to avoid propeller stalling, which is an undesired effect to be avoided in practice. Then, both the lift force and drag of each rotor can be expressed as a function of the rotational speed [3]: where c L i and c D i are the lift and drag coefficients that depend on the blade geometry. In the case of fixed-pitch propellers, c L i and c D i are constant. In a variable-pitch propeller, however, as the pitch angle α i ∈ R changes, the lift and drag coefficients change as well.
Two main approaches to characterize the dependency of c L i and c D i on α i can be found in [9], and they can be categorized into grey-box (e.g., [3]) and black-box (e.g., [5]) approaches. In this work, we follow the grey-box approach of [3], which makes use of the blade element theory and the momentum theory to obtain where R a is the propeller radius, ρ is the air density, while c T i (·) and c Q i (·) are dimensionless lift and drag coefficients, which are functions of α i . Following [9,23], which assume steadystate linearized blade aerodynamics, we consider the implicit relations where k 1 = 6 , and k 4 = 1 8 σC d0 i . In detail, σ = N b c πR a is the blade solidity, N b = 2 is the number of blades, c is the chord length, C l α i is the 2D liftcurve-slope of the airfoil section, comprising the rotor, and finally, C d0 i is the zero-lift drag coefficient (see [9] for further details). Without loss of generality, we consider a '+' quadrotor configuration (see Figure 1), so the overall input mapping is where l is the arm length. The mapping (12) is nonlinear in α i and Ω i , which are the actual control inputs for each motor i, for a total of eight inputs.

Explicit Lift and Drag Functions
Equations (10) and (11) are implicit functions of α i , so they are not practical for online computation. In the literature, such functions are solved iteratively by means of numerical approaches [9]. In this section, instead, we introduce two equivalent equations that determine the coefficients c T i and c Q i in terms of explicit functions of α i .
The lift function c T i (α i ) in (10) is expected to be an odd function of α i [9]. In fact, from (8), we obtain which is an explicit odd function of α i . The drag function c Q i (α i ) in (11) is expected to be an even function of α i [9]. Then, noting that c T i (|α i |) = sign(α i )c T (α i ), (11) can be rewritten as which is an explicit even function of α i . We also calculate in advance the derivatives of (13) and (14) with respect to α i , as they are going to be employed in the remainder of the paper. We have where the second row is obtained with the use of

Actuator Faults and Control Input Definition
A total of eight control inputs are considered, which are the pitch angles α i and the squares of the desired rotor speeds, denoted as u i , for i = 1, . . . , 4. All the control inputs are collapsed in the vectors In the absence of faults, the relation Ω 2 i = u i holds. For each actuator i, a possible LOE is considered, so that the resulting angular speed may be lower than the desired one, i.e., Each w i , for i = 1, . . . , 4, is an unknown and time-varying quantity which models the fault magnitude on the i-th actuator. It is constrained in w i ∈ [0, 1], where w i = 1 corresponds to the faultless scenario. For the sake of brevity, the fault magnitudes are collected in the vector

Wind and Fault Estimation
We assume the wind acts on the x E − y E plane only, that is, F E w = f wx e 1 + f wy e 2 . In other words, we neglect the updraft, and we claim that the horizontal wind does not generate any significant torque due to the symmetry of the vehicle. Thus, six target independent variables need to be estimated; these are f wx , f wy and w 1 , w 2 , w 3 , w 4 , and they are summarized in the extended uncertainty vectorw = col( f wx , f wy , w 1 , w 2 , w 3 , w 4 ). By handling the joint vectorw as an external disturbance, a disturbance observer for the estimationŵ ofw is then designed, which in our purpose can be defined as a dynamical system in the form˙ŝ for which ŵ −w → 0 as t → ∞, and whereŝ is the observer state,x = col(p F ,ṗ F , ω, η) andū = col(u, α) compactly refers to the quadrotor state and the control inputs, and where F(·) and H(·) are smooth function to design. Several disturbance observers have been proposed in the literature. Here, we follow the NDO proposed in [24], which can be summarized as follows.

1.
Define an auxiliary variable s which directly depends on the estimation targetw. Usually, the auxiliary variable s =w − λ(x,ū) is chosen, for some smooth function λ(·).

2.
Find the dynamicsṡ of the auxiliary variable s. In this step, the disturbance observer design problem is turned into a state observer design problem.

3.
Based on the auxiliary dynamicsṡ, design a state observer for the estimationŝ of s.

4.
Defineŵ by using the inverse function of the chosen auxiliary variable definition.
For design purposes, the model is conveniently rewritten in Section 3.1. Then, the design of the NDO is proposed in Section 3.2.

Remark 1.
The estimation of the actuator faults and the external force generated by the wind is obtained by the disturbance observer, where only kinematic data are required, i.e., without the need of measuring motor speeds, motor currents, etc.

Remark 2.
The design of λ(·) involves the solution of a set of partial differential equations.

Remark 3.
If the dynamics ofw is uncertain, the ultimate boundedness of ŵ −w is demanded instead of convergence.

Model Manipulation
Note that the input mapping (12) and the action of the fault (18) can be rewritten as and  resulting in the overall relation col(u z , τ B m ) = B c (α)diag(u)w. The quadrotor kinematics and dynamics in (5) can be equivalently expressed as d dt Then, defining z 1 = col(p F , η) and z 2 = col(ṗ F , ω), we can rewrite the overall quadrotor model in order to show that the non-linear model is affine with respect to the uncertaintȳ w. In fact, we haveż where and P(z 1 , u, α) = P 1 P 2 (z 1 , u, α) ∈ R 6×6 , with Finally, P is almost everywhere nonsingular, since where det(B c (α)) = 0 almost everywhere (assuming l > 0). Noting that each u i is positive in real flight conditions, cos(ϕ) cos(θ) is different from zero in non-aerobatic flight, and det(B c (α)) = 0, then the inverse of P is well-defined in practical flight conditions.

Disturbance Observer Design
The NDO is designed on model (26). Indeed, the NDO design in DOBC is well-known for control affine systems [25], so the idea is actually based on the affine uncertainty, which isw in our problem. Consider the auxiliary variable s =w − λ(z 1 , z 2 , u, α), where λ(·) is a function to design. For the sake of brevity, we denote λ = λ(z 1 , z 2 , u, α). By differentiation, we haveṡ where R s is a known term given by The disturbance observer is then designed aṡŝ whereŝ andŵ are estimations for s andw, respectively. The function λ(·) is investigated through the estimation error dynamics. On this purpose, let us consider the error variables s = s −ŝ andw =w −w. The error dynamics iṡs With the choice λ(z 1 , z 2 , u, α) = HP(z 1 , u, α) −1 z 2 , where −H is Hurwitz, the overall estimation error boils down to˙w = −Hw +ẇ, (36) which is an asymptotically stable linear time-invariant system, perturbed by the unknown inputẇ.

Control Law Design
A classical inner/outer loop approach is employed to calculate the desired virtual inputs [26,27], where the small angle approximation is adopted (ϕ ≈ 0 and θ ≈ 0). Moreover, in order to avoid a conflict between the inner and the outer loop, a spectral separation is required; to achieve this, in practice, a higher rate is employed for the inner loop with respect to the outer loop [28].

Outer Loop Design
Given the time-varying position references x F r and y F r (for x F and y F , respectively), the outer loop is designed to calculate the references ϕ r and θ r (for the angles ϕ and θ, respectively) in order to track the position reference. Exploiting the small angle approximation [26,27], that is, ϕ ≈ 0 and θ ≈ 0, we have Substituting (37) in (5), it follows that As the inner loop runs at a higher rate with respect to the outer loop, we can assume the inner loop is fast enough to track ϕ r and θ r [26]. Hence, by solving (38) for the variables ϕ and θ, we obtain the ideal references ϕ r and θ r for u z = 0, that is, where v x and v y are injection inputs, defined as The parameters α x1 , α x0 , α y1 , α y0 are chosen according to the desired pole placement problem. In fact, if the angles ϕ and θ follow references (39)-(41), we can show by substitution in (5) that the outer loop tracking error dynamics becomes: where e x = x F − x Fr , e y = y F − y Fr ,f wx = f wx −f wx andf wy = f wy −f wy .

Inner Loop Design
Let z Fr and ψ r be the reference signals for z F and ψ that we want to track, and let ϕ r and θ r be the reference signals for ϕ and θ (provided by the outer loop in (39)). From (5), it follows thatz so we can set the reference values u z,r and τ B m,r (for u z and τ B m , respectively) as The terms v z ∈ R and v η ∈ R 3 are auxiliary inputs, and they are set as where the scalars α z1 , α z0 and the matrices α η1 , α η0 are set according to the desired pole placement. If the variables u z and τ B m follow references (46)-(48), we can show by substitution in (5) that the inner loop tracking error dynamics becomes where e z = z F − z Fr and e η = η − η r .

Control Allocation
In this section, we show how the virtual input (i.e., u z and τ B m ) can be distributed among the motors (i.e., u) and the pitch servos (i.e., α). Let us consider a non-linear optimization problem in the optimization variables α and u: where J(α, u) in (51) is a cost function to be determined, while (52) and (53) represent generic non-linear and nonconvex constraints. In the literature, the control allocation problem in multirotors is usually tackled with lightweight methods, such as pseudo-inverse and redistributed pseudo-inverse [29,30]. Such methods are motivated by the need for high rates, but they are not effective in the presence of constraints and rate limits. Moreover, they need a significant simplification of the problem formulation (i.e., linear or quadratic cost functions, linear equality constraints, potentially also removing inequality constraints).
In the following, we set up the control allocation problem and we detail a method to solve it efficiently by means of a locally convex Quadratic Programming (QP) reformulation (as carried out, for example, in [31] for a marine vessel). The method that is detailed in this section to solve the non-linear allocation problem (51)-(53) represents a simplified version of the generic Sequential Quadratic Programming (SQP) strategy (see, for example, [32]) to solve non-linear optimization problems. The main differences are • SQP employs the Hessian of the Lagrangian in the cost function, while we employ the Hessian of the original cost function. In other words, we always set the Lagrange multipliers' initial guess to zero for simplicity, so that the Hessian of the Lagrangian coincides with the Hessian of the original cost function. • Once the search direction is found, a full step is performed in the proposed method, and no merit function is computed. • The stopping criterion in the proposed method boils down to performing just a single iteration (i.e., a single QP is solved as in Real-Time Iteration [33], instead of a sequence of QP problems as done in SQP), and so no convergence criteria are needed.
These modifications introduce some loss of accuracy in the solution, but they reduce the number of evaluations, in order to perform the allocation faster in view of online implementation.

Equality Constraints
The equality constraints (52) originate from (22)-(23), namely where Please note that c L (α i ) and c D (α i ) have been defined in (8), while c L (α i ) and c D (α i ) have been defined in (15). The time-varying numeric values for u z , τ B m , and w in (54) are given by (46) (u z,r , τ B m,r ) and (33) (ŵ), respectively; in other words, the actual lift force and torques must satisfy those required by the inner loop, also taking into account the actual fault estimation.
Please note that the equality constraints are non-linear and nonconvex with respect to the optimization variables α and u. Non-linear optimization problems are hard to solve online due to computational limits, as well as because the existence of the solution is not guaranteed in general. To solve efficiently the problem online, the equality constraints (52) are linearized, as in SQP, so that a QP problem can be obtained. Using the first-order Taylor approximation around (α 0 , u 0 ), we have where α 0 and u 0 represent the current value of α and u, respectively. Then, (54) can be approximated by where ∆α = α − α 0 , ∆u = u − u 0 , and As w is not available, to build the control allocation problem online we replace it in (60) with its estimation contained inŵ.

Remark 4.
In the case of a centralized motor, it suffices to add the constraint ∆u 1 = ∆u 2 = ∆u 3 = ∆u 4 ; provided that the same initial value is employed for u, then u 1 = u 2 = u 3 = u 4 holds. In other words, (59) is replaced by

Inequality Constraints
The inequality constraints of interest to be included in (53), in practice, are saturation constraints and rate limits, i.e., Such inequality constraints are affine with respect to the optimization variables α and u, so they are tractable in online optimization.
Replacing (52) with its approximation (59) (or (61), for a single centralized motor) and expressing the constraints (62) in the variable x = col(∆α, ∆u), we obtain an optimization problem in the form min where are defined to express the inequality constraints (62) in a more compact way. Please note that (63) is a QP problem by construction, provided that the cost function J(x) is quadratic.

Cost Function
The main term to be considered in a cost function for UAVs, especially in the case of multirotors, is the energy consumption E(x). A quadratic form for E(x) is often assumed in the UAV literature [29,30]; such a choice is mainly motivated by simplicity, as it usually leads to a QP problem. Motivated by the need for autonomy in UAVs, more detailed representations, coming from physically accurate consumption models, have been proposed in the literature as well. The most employed models are second-order polynomials [9] and irrational functions [34], as well as more generic non-linear functions [6]. All of the previous choices fit well with the proposed framework, as it is always possible to locally approximate the cost function with a quadratic cost function, as done for example in [31]. For the purpose of this work, the steady-state thruster consumption from [34] is employed. Given a generic matrix M = [m ij ] ∈ R m×n and k ∈ R, let us introduce, with a slight abuse of notation, the short |M| k = [|m ij | k ] (i.e., the k-th power of the absolute values of the entries of M, in an element-wise fashion). Then, the steady-state thruster consumption for a quadrotor can be rewritten as where |u| 3/4 and |u T | 3/4 should be intended element-wise as already anticipated. Furthermore, q u,i is assumed to be constant for simplicity (see [34] for further details) and with q u,i > 0 for i = 1, . . . , 4. Hence, we can define the cost function where Q ∆α , Q ∆u , Q α , Q u are symmetric positive definite matrices. The first term of (67), which takes into account energy consumption, is detailed in (65). The second and the third quadratic terms of (67) discourage variations in the control inputs to minimize actuator wear and tear. The last term of (67) defines preferred values α p for the pitch angles α. This term becomes necessary in order to compensate for the first term (i.e., the one that discourages energy consumption), because the actuator is more energy-efficient when the pitch angle is large in magnitude. Thus, in the absence of the fourth term, the pitch angle is often stuck to its saturation limit to minimize energy consumption. In such a case, the only way to increase the lift further is to increase the motor speed, as the pitch angle cannot be increased in magnitude, so the maneuverability of the vehicle worsens.
Considering the second-order Taylor approximation for E(x), i.e., and performing some manipulation, (67) becomes where c is a constant term (i.e., a term that does not depend on x), which is irrelevant for the solution of the optimization problem, so it can be neglected.

Scaling and Infeasibility
The problem formulation according to (59) (or (61), for a single centralized motor), (64), (69), that is, shows two drawbacks. First of all, the existence of a feasible solution is not guaranteed, so problem infeasibility must be dealt with. Moreover, the problem is ill-conditioned, because the magnitude of ∆α and ∆u is very different; according to (17)- (18), α is typically in the range of decimals (in radians), while u is typically in the range of tens of thousands (in radians per second squared), and the same occurs for ∆α and ∆u.
To improve the problem conditioning, diagonal scaling is performed [32]. We define a new optimization variable where P = diag(max(|x|, |x|)) is a diagonal matrix such that the components of z share the same magnitude, as they become constrained in [−1, 1] according to the rate limits in (62) (please note that the box constraint estimation [−1, 1] is conservative; tighter constraints may apply).
In order to deal with infeasible problems, we soften the equality constraint by adding a vector of unconstrained artificial variables ζ into the equality constraint, while minimizing the norm of the artificial variable. The overall problem (70), after the change of variables and the softening of equality constraints, becomes min z,ζ which can be finally rewritten in the classical form that is usually required by QP solvers, i.e., y * = arg min y 1 2 y TȞ y +f T y s.t.Ǎ eq y =b eq Ay ≤b y ≤ y ≤ y, where ζ and ζ are large enough to guarantee feasibility (e.g., ζ = −ζ = M, M → ∞), whileǍ anď b are empty matrices, because no inequality constraints are left, with an exception made for the upper and lower bounds that are usually introduced separately. Once the solution y * of (72)-(74) is calculated, the solution to the control allocation problem is Remark 5. The choice of a quadratic penalty for ζ is made for simplicity, as it introduces a small number of additional variables and it makes the cost function smooth and hence easier to solve online. This strategy is similar to the quadratic penalty method in [32] and, more precisely, to the weighted least squares in [35]. As a drawback, such a penalty term is not exact, i.e., the solution of the penalized problem is not equivalent to the solution of the original problem [32,36]. In particular, the higher the penalty through Q ζ , the more the solutions of the two problems coincide, but also the worse is the conditioning. Thus, Q ζ is a trade-off between ideal solution quality and acceptable conditioning. We claim that this issue is not crucial, because several approximations are introduced to simplify the problem, and so the solution is not exact anyway. Alternatively, it is possible to introduce additional positively constrained slack variables and to employ exact penalty functions, such as 1 and ∞ penalties, at the price of dealing with the nonsmoothness of the cost function (see [32] for further details).

Numerical Simulation
A numerical simulation is performed using MATLAB, where the inner loop and the outer loop run at 400 Hz (i.e., ArduPilot's loop rate) and 10 Hz, respectively. A zeroorder hold is performed between consecutive samples. The quadrotor parameters are set according to [18] and they are reported in Table 1. Different from [18], the friction is neglected and the thrust-to-weight ratio is slightly increased (from 2.5 to 3.131) in order to face severe faults. The control parameters were set empirically through a coarse grid search (as such, better performances could be achieved with a finer tuning). The closed-loop poles for the inner loop are placed in −2 and −2.1 for z F , in −10 and −11 for ϕ and θ, and finally in −5 and −5.1 for ψ. The outer loop poles are instead placed in −1 and −2 for both x F and y F . A properly scaled white Gaussian sensor noise is simulated according to an MPU-9250 IMU TDK InvenSense [37], which is a common device in quadrotor flight controllers (i.e., in Pixhawk 4).

Independent Motors' Speed
Let us first consider a VPQ equipped with four independent motors and four independent pitch servos. The simulation is performed for a total of 30 s under the external wind force reported in Figure 3, which also reports the correspondent estimations provided by the observer, where the gain matrix has been set as H = 5I, with I being the identity matrix. Due to the stochastic nature of the wind, achieving the boundedness of the estimation error replaces the ideal goal of a perfect convergence, which is here achieved in the case of bounded derivatives of the wind force. The estimation of the faults and the wind is practically started after a finite time from the beginning of the flight (i.e., after 1 s of simulation in this example), in order to make the initial conditions of the observer practically irrelevant.
Faults components and their estimations are reported in Figure 4. Two trapezoidal faults are injected. The first one affects motor 1 after 10 s (the LOE magnitude is 60%, i.e., w 1 = 0.4) and a second one affects motor 3 after 15 s (the LOE magnitude is 30%, i.e., w 3 = 0.7). In the range t ∈ [0, 1], eachŵ i is set equal to 1, matching the reasonable assumption of no faults at the initial time (or, without loss of generality, we assume that we know the fault status at the beginning of the flight). Note that, since each estimation is provided independently from the NDO, the same results can be achieved in the case of simultaneous faults.
The tracking reference components and the related state variables are reported in Figure 5. The references on x F and y F are sinusoidal signals, the reference on z F forces the quadrotor to face both ascending and descending profiles, while the yaw angle is kept at zero. The system is forced to change the velocity constantly in both magnitude and direction, with a linear speed magnitude up to 7 m/s. All the references are smoothly tracked with an acceptable error (the plotted variables are the actual state variables instead of the noisy ones), with the exception of the yaw component (ψ), which exhibits a temporary increasing error when the fault worsens.
The references ϕ r and θ r generated by the outer loop are reported in Figure 6. Please note that, when the fault is rapidly worsening, the tracking capabilities in the inner loop are affected. The inner/outer loop approach, which employs a higher frequency of the inner loop and a zero-order hold, allows us to track the references, despite the presence of noise, faults, and disturbances.
The control inputs u and α are reported in Figure 7. All the actuators are subject to the same saturation constraints. This allows us to plot the bounds with horizontal dashed lines in Figure 7 (upper bounds in this case). In particular, the saturation limits have been set as α = col(α min , α min , α min , α min ) α = col(α max , α max , α max , α max ) u = col(u min , u min , u min , u min ) u = col(u max , u max , u max , u max ) ∆α = −∆α = col(∆α max , ∆α max , ∆α max , ∆α max ) where u min , u max , α min , and α max are detailed in Table 1. Considering α = α, u corresponds to a thrust-to-weight ratio equal to 3.131.
At the beginning, when no fault is present, the opposite actuators (1 and 3, 2 and 4) show similar motor speeds and pitch angles, due to symmetry. As soon as the fault on motor 1 occurs, the opposite actuators 1 and 3 behave differently; in response to the severe fault acting on motor 1, the control allocation algorithm demands a high speed and a large pitch angle to compensate for the missing lift force.

Centralized Motor
In this section, we consider a VPQ equipped with a single centralized motor and four independent pitch servos. The simulation conditions and the wind components are the same as the previous simulation, including the same random seed. Figure 8 reports the wind effect and its estimations provided by the observer; no noticeable differences can be noted, as expected, because the fault estimation does not directly depend on the actual inputs. The same consideration holds for the faults and their estimations ( Figure 9). The tracking reference components and the related state variables are reported in Figures 10 and 11. The performances are similar; minor differences in the outer loop tracking capabilities can be detected.
Finally, the control inputs u and α are reported in Figure 12.
First of all, please note that the motor speeds coincide, as expected. Furthermore, in response to the fault acting on motor 1, the control allocation algorithm demands a high motor speed to compensate for the missing lift force. As there is only one centralized motor, the speed of every propeller must increase. Analogously, please note that, the larger the fault, the higher the average value of the pitch angle.   As the motor speed is higher on average in the case of the centralized motor (Figures 7 and 12), the energy consumption is higher. This is unavoidable, because in the case of the centralized motor, the input redundancy is lower than in the case of independent motors. In Table 2, we summarize the overall normalized (i.e., Q u = I) consumption during the flight, as proposed in (65).

Conclusions
Unlike standard quadrotors that are equipped with four fixed-pitch blades and four independent motors, VPQs can vary both the rotation speed and the blade pitch of each propeller, thus increasing the degrees of freedom for control purposes. In addition to higher thrust rate of change, reverse thrust and reverse flight capabilities, and scaling well with size, VPQs enable for fault-tolerant control strategies thanks to their inherent input redundancy, whether they are equipped with four independent motors or a centralized one. We have shown that a single observer can estimate both actuator loss of effectiveness and the external disturbance given by horizontal wind. According to the disturbance observerbased control framework, the estimation of the wind is fed forward in the outer (position) control loop. The estimation of the actuator loss of effectiveness is instead employed in the control allocation algorithm. It manages each propeller's pitch and motor speed to generate the force and the torques commanded by the inner (attitude and altitude) control loop. Simulation results show that the estimation of both the disturbance and the actuator loss of effectiveness is practically feasible in the presence of conventional measurement noise in commercial devices. We remark that the fault estimation is performed from kinematic data of the onboard inertial measurement unit only, without the need to measure the motor speed nor the current drawn by the actuators.
Comparing centralized and independent motors, the former is less efficient in terms of energy consumption, especially in the case of faults, because all the propellers must spin faster. Clearly, many other factors should be considered to evaluate the energy consumption properly, such as possible weight and cost reduction. In any case, both centralized and independent motors can tolerate actuator faults using the proposed scheme. The proposed control scheme could also deal with pitch lock-in-place, provided that reliable information on the stuck fault is available (e.g., measuring the current pitch angle or estimating it through a fault detection algorithm). Practically, it is sufficient to constrain the faulty pitch servo to the current value, and the control allocation algorithm can accommodate for the stuck fault using the remaining control inputs. However, in the case of a single centralized motor, the reduced amount of redundancy makes it less robust to lock-in-place faults and severe loss of effectiveness.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript: