A Comparative Study of Nonlinear MPC and Differential-Flatness-Based Control for Quadrotor Agile Flight

Accurate trajectory-tracking control for quadrotors is essential for safe navigation in cluttered environments. However, this is challenging in agile flights due to nonlinear dynamics, complex aerodynamic effects, and actuation constraints. In this article, we empirically compare two state-of-the-art control frameworks: the nonlinear-model-predictive controller (NMPC) and the differential-flatness-based controller (DFBC), by tracking a wide variety of agile trajectories at speeds up to 20 m/s (i.e., 72 km/h). The comparisons are performed in both simulation and real-world environments to systematically evaluate both methods from the aspect of tracking accuracy, robustness, and computational efficiency. We show the superiority of the NMPC in tracking dynamically infeasible trajectories, at the cost of higher computation time and risk of numerical convergence issues. For both methods, we also quantitatively study the effect of adding an inner loop controller using the incremental nonlinear dynamic inversion method, and the effect of adding an aerodynamic drag model. Our real-world experiments, performed in one of the world’s largest motion capture systems, demonstrate more than 78% tracking error reduction of both NMPC and DFBC, indicating the necessity of using an inner loop controller and aerodynamic drag model for agile trajectory tracking.

However, NMPC is computationally extremely demanding compared to the state-of-the-art non-predictive method: the differential-flatness-based controller (DFBC) [12,13].This method also shows great potential in accurately tracking agile trajectories autonomously.A recent study has used DFBC to track trajectories up to 12.9 m/s with 2.1 g accelerations, with only 6.6 centimeters position tracking error [13].Although the state-of-the-art DFBC has not achieved agile flights as fast as NMPC showed in [4], its high computational efficiency and tracking accuracy render the necessity of applying NMPC for agile trajectory tracking questionable.Therefore, a comparative study of NMPC and DFBC is needed to provide insights for future research on fully autonomous agile flights, in order to further improve their efficiency and reliability.

B. Contribution
This article presents the first comparative study of the two state-of-the-art tracking controllers: NMPC and DFBC method, in fast trajectory tracking with speed up to 20 m/s (i.e., 72 km/h) and 5 g accelerations.Specifically, we select the NMPC method that uses the full nonlinear dynamics with proper actuator bounds and regards single rotor thrusts as control inputs.This NMPC has been applied in previous works, such as [3,5].For a fair comparison, we improve the DFBC method used in [12,13] with a control allocation approach using constrained quadratic programming [14] also to consider control input limits.
All experiments are conducted in both simulations and the real world.The simulations compare in both ideal and practical conditions with model-mismatch, estimation latency, and external disturbances.The real-world experiments are conducted in one of the world's largest motion capture systems, with 30×30×8 m 3 flight volume.Multiple agile trajectories are selected as reference, including not only dynamically feasible trajectories, but also dynamically infeasible trajectories that require thrusts exceeding the maximum capacity of the quadrotor motors, which is likely to happen in high-speed flights due to model mismatch.These tests investigate the performance of both methods in the presence of significant high-speedinduced aerodynamic effects, inevitable system latency, and on an onboard embedded computer.
The experimental results reveal that NMPC is considerably more computationally demanding, and more prone to suffer from numerical convergence issues in the presence of large external force disturbances.However, NMPC also excels in tracking dynamically infeasible trajectories, making it more suitable for tracking time-optimal trajectories that violate the rotor thrust constraints.
In addition, this study also highlights the importance of implementing an inner-loop controller for robustification.Specifically, we select the incremental nonlinear dynamic inversion (INDI) method as the inner-loop angular controller, thanks to its simplicity in implementation and demonstrated robustness in various real-world experiments [15,16] including a combination with DFBC [13].As for NMPC, differently from the state-of-the-art using a PID as the low-level control [4], we propose a method to hybridize NMPC with INDI that considers the real input limits of the quadrotor instead of constraints on virtual inputs.Real-world flight results demonstrate more than 78% position tracking error reduction of NMPC and DFBC with an INDI inner-loop.We also reveal that a well-selected inner-loop controller is more crucial than simply considering the aerodynamic effects, as is compared in Fig. 1.
Apart from the technical contribution, this paper can also be regarded as a tutorial for non-expert readers on agile quadrotor flight.We encourage the practitioners to use the presented results as a baseline for further development of both DFBC and NMPC approaches.

II. RELATED WORK
In this paper, we classify the trajectory tracking controller into two categories: non-predictive and predictive methods.While the predictive methods encode multiple future timesteps into the control command, the non-predictive methods only track a single reference step.In the following, we review works towards improving quadrotor trajectory tracking accuracy from these two different categories.A more comprehensive survey of quadrotor position and attitude control methods can be found in the literature (e.g., see [17,18]).

A. Non-predictive Quadrotor Trajectory Tracking Control
Unlike most fixed-wing aircraft, quadrotors are inherently unstable.Therefore, the initial work of quadrotor control aimed at achieving stable hovering and near-hover flights.Thanks to the small-angle assumptions in these conditions, linear control methods such as PID and LQR demonstrate sufficiently good performance (see, e.g., [19,20]).
However, as the requirements for agile flight emerges, these assumptions are no longer valid.Nonlinearities from the attitude dynamics are the first problem to tackle.For this reason, nonlinear flight controllers are proposed, such as feedback linearization [21] and backstepping [22].In order to cope with the singularities of Euler angles as the nonlinear attitude representation, quaternions are widely adopted to parametrize the attitude [23].In addition, the authors of [24] propose the geometric tracking controller to directly control the quadrotor on the manifold of the special Euclidean group SE(3), showing almost globally asymptotic tracking of the position, velocity, and attitude of the quadrotor.
A seminal work [25] reveals that quadrotors are differentially flat systems.By virtue of this property, given the time-parameterized 3D path, one can derive the reference attitude, angular rate, and accelerations.These references can be sent to a closed-loop flight controller as feedforward terms, while additional feedback control is required to address model mismatch and external disturbances.As such, differentialflatness based controller (DFBC) has significantly improved the tracking performance at relatively high speeds [12,13].
As the flying speed increases, a quadrotor starts experiencing non-negligible aerodynamic effects, including drag [26], aerodynamic torque [27], and variation of thrusts [28].Authors of [12] show that the aerodynamic drag does not affect the differential flatness of a quadrotor.Thus they adopt a first order aerodynamic model with feedforward terms derived from the reference trajectory to improve the tracking performance.Since an accelerometer can directly read external forces, [13] leverages accelerometer measurements to improve the tracking accuracy instead of resorting to an aerodynamic drag model.This method also demonstrates remarkable performance in disturbance rejection and platform adaptability.
Effectively handling control input limits is a remaining challenge for non-predictive methods, including DFBC.So far, existing methods have prioritized the position tracking over heading using various approaches, such as redistributed pseudo inversion [29], weighted-least-square allocation [30], controlprioritization method [31,13], and constrained-quadratic-programming allocation [32].While these methods can mitigate the actuator saturation effect when the trajectories are dynamically feasible, its performance in tracking dynamically infeasible trajectories is still questionable.In this work, we will push the limit of DFBC to conduct agile trajectory tracking tasks to a much higher flight speed and acceleration than the state-of-the-art, and study its performance in tracking dynamically infeasible trajectories where violating rotor thrust limits becomes inevitable.

B. Model Predictive Control for Quadrotor Trajectory Tracking
Model predictive control (MPC) is a prevalent method in robots thanks to its predictive nature and ability to handle input constraints [10,33].It generates control commands in receding horizon fashion, which minimizes the tracking error in the predicted time horizon by solving constrained optimization problems.
However, MPC is very computationally demanding compared with non-predictive methods.Especially, using nonlinear-MPC (NMPC) with a full-state nonlinear model of a quadrotor was computationally intractable onboard earlyage unmanned-aerial vehicles (UAVs).For this reason, linear-MPC (LMPC) is adopted in many studies for only position control [34], or control linearized model based on smallangle assumptions [35].Therefore, these LMPC methods cannot fully capture nonlinearities in rotational dynamics, and underperforms NMPC methods [10,33] Separating flight controllers into high-level position control and low-level attitude control is another common approach to simplify the model and alleviate the computational load of NMPC [6,7,9].However, in such a cascaded control architecture, the high-level controller cannot precisely capture the real quadrotor capability because they often ignore the effects of system limitation (such as [6]), or introduce a virtual constraint on states [9].Consequently, the command fed into the lower-level controller is either too conservative to achieve agile flights, or over-aggressive that causes instability.
Thanks to the recent development in hardware and nonlinear optimization solvers [36][37][38], running NMPC with a nonlinear full dynamic model becomes computationally tractable on an embedded computer.Hence, recently, some studies start using the full-nonlinear dynamic model, and regarding single rotor thrusts as inputs for NMPC, thus are able to fully exploit quadrotors capability [3][4][5].These methods either directly use the optimized single rotor thrust commands [3], or send intermediate states from the solution (such as the angular rates) to a low-level controller [4,5].A recent study [4] demonstrates the ability of the full-model NMPC with a PID low-level controller in tracking a pre-planned race trajectory at speed up 20 m/s which surpasses the top speed of 12.9 m/s reported in [13] using DFBC, in spite of a much larger tracking error (0.7 m v.s.0.066 m).
Although running NMPC is realizable on modern embedded computers, it still requires significantly more computational resources than non-predictive methods represented by DFBC.For this reason, NMPC may suffer from numerical convergence issues, especially when the platform lacks a sufficient computational budget.Moreover, the advantage of NMPC becomes questionable as DFBC can also address input limits using the control allocation technique and generate feedforward control leveraging differential-flatness property.Therefore, it is necessary to compare these two methodologies and understand at what conditions each approach is preferable to provide insights and recommendations for future applications.

A. Notations
Throughout the paper, we use subscription r to indicate the reference variables derived from the reference trajectory.Subscription d indicates the desired value, which is calculated from a higher-level controller.We use bold lowercase letters to represent vectors and bold uppercase letters for matrices; otherwise, they are scalars.Two right-handed coordinate frames are used in this paper: they are the inertial-frame F I : {x I , y I , z I } with z I pointing upward opposite to the gravity, and the body-frame F B : {x B , y B , z B } with x B pointing forward and z B aligned with the collective thrust direction (see Fig. 2).Vectors with superscription B are expressed in the body-frame; others without any superscription are expressed in the inertial-frame.The rotation from F I to F B is represented by rotational matrix R(q) = [x B , y B , z B ] ∈ SO(3) parameterized by quaternion q = [q w , q x , q y , q z ] T ∈ S 3 .We use subscript {x, y, z} to represent the imaginary components of a quaternion, namely q x,y,z = [q x , q y , q z ] T .Operator diag(A 1 , A 2 , ..., A n ) denotes a diagonal matrix with scalars or matrices (A 1 , A 2 , ..., A n ) as diagonal entries.
B. Quadrotor Model 1) Quadrotor Rigid-Body Model: The quadrotor model is established using 6-DoF rigid body kinematic and dynamic equations.For translational dynamics, we have where ξ denotes the position of the quadrotor center of gravity (CoG); T and m are the collective thrust and total mass respectively; g ∈ is the gravitational vector; f a indicates the exogenous aerodynamic drag force during high-speed flights.
The rotational kinematic and dynamic equations are expressed as where ⊗ denotes the quaternion multiplication operator.Ω is the angular velocity of F B with respect to F I .Its derivative, namely angular acceleration, is denoted by α .Throughout the paper, we use its projection on F B , namely Ω B = [Ω x , Ω y , Ω z ] T , since its directly measurable from the inertial measurement unit (IMU).I v indicates the inertia matrix of the entire quadrotor.τ is the resultant torque generated by rotors.d τ is the model uncertainties on the body torque, which can come from high-order aerodynamic effects, center of gravity bias, or distinction among rotors.
The collective thrust and rotor generated torques are functions of rotor speeds: where represents the thrust generated by each rotor and • indicates the Hadamard power.c t is the thrust coefficient.ω is the vector of angular rates of each propeller.G 1 to G 3 are matrices defined as where c q is the torque coefficient; β and l are geometric parameters defined as per Fig. 2. I p is the inertia of the rotor around z B .Terms G 2 ω and G 3 (Ω)ω are torques respectively due to angular acceleration of rotors and gyroscopic effects, which are usually neglected for controller design.However, in Sec IV-C, we will revisit the INDI method [13] for angular acceleration control that takes into account effects of the inertial torque term G 2 ω.
2) Aerodynamic Drag Model: Quadrotors during highspeed flight experience significant aerodynamic drag forces, which need to be precisely modeled to improve tracking accuracy while minimizing the computational overhead.In this work, we use an aerodynamic drag model which captures the major effects, and is proved effective in works such as [12].
where [v x , v y , v z ] = R(q) T ξ (i.e., the projection of velocity in the body frame; here we assume zero wind-speed).k d,x,y,z and k h are positive parameters can be identified from flight data.
IV. METHODOLOGIES This section elaborates the two control methods compared.A nonlinear NMPC method is selected that considers the thrust limits of each rotor, the full nonlinear dynamics, and the aerodynamic effects.As for the DFBC method, we improve the state-of-the-art such that these factors are also addressed, which ensures a fair comparison with NMPC.Finally, both methods are augmented with an INDI controller [15] to convert the single rotor thrust commands to rotor-speed commands, while improving the robustness against model uncertainties and disturbances on the rotational dynamics.The control diagrams of both methods are illustrated in Fig. 4 and Fig. 3.

A. Nonlinear Model Predictive Controller
NMPC generates control commands by solving a finite-time optimal control problem (OCP) in a receding horizon fashion.Given a reference trajectory, the cost function is the error between the predicted states and the reference states in the time horizon, meaning that multiple reference points in the time horizon are used.In order to conduct numerical optimizations, we discretize the states and inputs into N equal intervals over the time horizon τ ∈ [t, t + h] of size dt = h/N with h denoting the horizon length, yielding a constrained nonlinear optimization problem: Note that, we use the thrust of rotors u defined in ( 4) and ( 5) as the control input of NMPC.The state vector is defined as x = ξ ξ q Ω B , and The reference state vector x r and input u r can be obtained from a trajectory planner which generates full states such as the one introduced in [4].Function f (x k , u k ) is the discretized version of nonlinear quadrotor model ( 1)-( 5).The same as many other works (see, e.g., [3,39]), we omit G 2 and G 3 related terms in (4) as their effects are negligible.
x init is the current state estimation when solving the OCP.u min and u max are minimum and maximum values of motor thrusts.Constraints on angular velocities are also added, which is found beneficial for improving the stability of NMPC according to [5].Note that in the above optimization problem, the following abuse of notation is used when calculating quaternion error: The above NMPC solves the full nonlinear model of a quadrotor, instead of resorting to a cascaded structure, or linear assumptions.In this paper, this quadratic nonlinear optimization problem is solved by a sequential quadratic programming (SQP) algorithm executed in a real-time iteration scheme [40].We implement this algorithm using ACADO [41] toolkit with qpOASES [42] as the solver.More implementation details are given in Sec.V.

B. Differential-Flatness Based Controller
Quadrotors are differentially flat systems [25], namely, all the states and inputs can be written as algebraic functions of the flat outputs and their derivatives.This allows a direct mapping from the flat outputs (positions ξ and heading ψ) to the angular rates and angular accelerations, which are leveraged by DFBC as feed-forward terms to improve the tracking accuracy.
In the following, we introduced the DFBC method improved from a previous work [25], where the original geometric attitude controller is replaced by the tilt-prioritized method [31].We also use the quadratic-programming based control allocation [14] to address input constraints.These modifications are beneficial in tracking dynamically infeasible trajectories, as will be discussed in Sec.VI-D.Fig. 4 presents an overview of this method.
1) Desired Attitude and Collective Thrust: First of all, we calculate the desired acceleration from a PD controller: where K ξ and K v are positive-definite diagonal gain matrices.Then from ( 1) and ( 9), the desired thrust T d and thrust direction z B,d are obtained as Given reference heading angle ψ r , we get an intermediate axis x C,d , by which the desired attitude can be obtained by the following equations: where q d is the desired attitude expressed by the quaternion.
2) Reference Angular Velocity and Acceleration: We leverage the differential flatness property of a quadrotor to derive the reference angular velocity and angular accelerations.Using them into the attitude control can help in tracking jerk and snap (3rd and 4th order derivatives of position ξ), which is found beneficial to the tracking performance [13].
Taking the derivative of both sides of (1) and assuming a constant external aerodynamic force f a , we have Then given reference jerk ... ξ r and substitute it for ... ξ in ( 18), we get by which the reference angular velocity can be obtained as where ψ r is the reference heading angle.
We further take the derivative on both sides of ( 18) and uses snap reference .... ξ r to replace .... ξ , yielding Then the desired angular acceleration can be obtained as Note that in (20) to (22) we use the current attitude {x B , y B , z B }, angular velocity Ω, and collective thrust T instead of their references.Hence, the drone still follows the reference jerk ... ξ r and snap .... ξ r even though its attitude, angular rates, and collective thrust have been deviated from the reference.
Above derivations for Ω B r and α B r need the value of collective thrust T and its derivatives.While T can be calculated from (4)(5) with measured rotor speed ω, its derivatives ( Ṫ and T ) are unable to be directly measured.For this reason, we approximate them by using reference jerk ... ξ r and snap .... ξ r .Multiplying (dot-product) both sides of (18) by z B , we have and its derivative Equation ( 23) and ( 24) are then substituted into ( 18)-( 22) to calculate Ω B r and α B r .3) Tilt-Prioritized Attitude Control: Quadrotors use rotor drag torques to achieve heading control.The control effectiveness of heading is around one order of magnitude lower than pitch and roll.As a consequence, heading control is prone to cause motor saturations.Fortunately, the thrust orientation of a quadrotor (namely tilt) is independent of its heading angle.Thus tilt-prioritized control has been proposed in [31] that regulates the reduced-attitude (pitch and roll) error qe,red and yaw error qe,yaw separately as follows: [q e,w , q e,x , q e,y , q e,z ] qe,red = 1 q 2 e,w + q 2 e,z   q e,w q e,x − q e,y q e,z q e,w q e,y + q e,x q e,z 0 qe,yaw = 1 Subsequently, the desired angular acceleration can be obtained by the following attitude control law: qe,red + k q,yaw sgn(q e,w )q e,yaw where k q,red and k e,yaw are positive gains for reduced-attitude and yaw control respectively.Setting a relatively high k q,red over k e,yaw is helpful in improving position tracking accuracy while preventing input saturations.α B r in (28) performs as a feed-forward term.It is worth noting that, the inclusion of α B r seems to be reasonable from a theoretical perspective [31], although removing it has been found to have almost no effect in our real-world experiments.
4) Quadratic-Programming-Based Control Allocation: The control allocation module generates thrust commands of each individual rotor from desired collective thrust T d and angular acceleration α B d .The same as NMPC, we also neglect the G 2 and G 3 terms in (4).Then from (3) and ( 4), we obtain the direct-inversion control allocation: This, however, does not consider input limits and may cause loss of control.For instance, an excessive collective thrust command can saturate all motors and consequently disable the attitude control.
An effective alternative to address saturations is the quadratic-programming based allocation that solves the fol-lowing optimization problem: where W ∈ R 4×4 is a positive-definite diagonal weight matrix.Each diagonal entry respectively indicates the weight on the thrust, pitch, roll and yaw channels.Generally, setting a relatively high weight on pitch and roll (see Table .I) is advantageous to prevent quadrotor loss-of-control when motor saturations are inevitable (e.g., tracking dynamically infeasible trajectories).If the solution is originally within control bounds, the result is the same as the direct-inversion allocation from (30).As for the implementation details, we solve this quadratic programming problem using an Active-Set Method from the qpOASES solver [42].

C. Incremental Nonlinear Dynamic Inversion
After obtaining thrust commands from (10) or (31), one can use (5) to directly obtain the rotor speed command.However, the above-mentioned controllers neglects the unmodeled term d τ in the rotational dynamics (3), which are found detrimental to the overall control performance.
Modeling d τ is very challenging for real-world systems.Therefore, we resort to incremental nonlinear dynamic inversion (INDI), a sensor-based controller that uses instantaneous sensor measurement, instead of an explicit model, to represent system dynamics, thus being robust against model uncertainties and external disturbances.The performance and robustness of INDI have been demonstrated in previous research (see, e.g., [15,13,16]) with proven stability [43].
We use INDI as the inner-loop controller of both NMPC and DFBC for fair comparisons.The hybridization of INDI and DFBC is similar to a related work [13], except for the attitude controller and control allocation introduced in the previous section.Here, a method is proposed to effectively combine INDI with NMPC to improve its robustness against rotational model uncertainties.
After knowing the single rotor thrust command u DFBC or u NMPC from ( 31) and ( 10) respectively, we can retrieve the desired collective thrust, and desired angular acceleration using ( 3) and ( 4), yielding Td Note that, for the DFBC method, Td and αB are different from those derived from ( 14) and ( 28) if the optimal cost of ( 31) is non-zero.Then from INDI, we can get the desired body torque (see (30) and ( 31) of [13] for detailed derivations) Hence, the effect of unmodeled d τ is captured by filtered angular acceleration measurement ΩB f and filtered body torque τ f , where is calculated from rotor speed measurements.Ω B f and ω f are low-pass filtered body-rate and rotor speeds with the same cut-off frequencies.Thus they have a similar amount of delay and can be synchronized.In this paper, we use a second-order Butterworth filter with a 12 Hz cut-off frequency.Note that we use subscript k − 1 to indicate the last sampled variable, and ∆t is the sampling interval.Ḡ1 and Ḡ2 respectively indicate matrices formed by the last three rows of G 1 and G 2 .
Then from (14), the following equation is obtained to solve rotor speed command ω c : with the only unknown ω c which can be solved numerically.
Differently from [13], the motor time constant is not needed in (35).
The INDI inner-loop controller converts the original thrust inputs from the high-level controller to rotor speed commands.Because this conversion is through algebraic equations, system delay typically seen in the cascaded control structures can be effectively circumvented.The advantage of INDI over classical PID inner-loop controller will be demonstrated in Sec.VII.

V. IMPLEMENTATION DETAILS
Before elaborating on the simulation and real-world experiments, we introduce the implementation details and experimental setups.Both NMPC and DFBC flight controllers are implemented in an open-sourced flight stack Agilicious 1 programmed in C++, with ACADO [41] toolkit and 1 https://agilicious.dev/qpOASES [42] as external libraries.A wrapper in ROS environment is also written which allows data logging, interfacing, and network communication.The controller gains are listed in Table I and the inertial and geometric parameters of the tested quadrotor are listed in Table II.These parameters are used for both simulation and real-world experiments.
In the simulation experiments, the flight software runs on a laptop at 300 Hz while the frequency of NMPC is limited to 100 Hz for consistency with the real-world experiments.The rotor speed command is sent to a simulator included in Agilicious that uses a 4th-order Runge-Kutta integrator to propagate quadrotor dynamic equations ( 1)-( 4) at 500 Hz.To simulate the motor dynamics, rotor speed commands generated by the controllers pass through a low-pass filter with a 30 ms time constant.The drag model ( 9) is also included to study the effect of aerodynamics, of which the parameters are included in Table II.The quadrotor states from the simulator are directly fed into the controllers unless a dedicated state estimation error is simulated.
The real-world experiment is performed in an indoor arena equipped with a motion capture system (VICON), with a 30 × 30 × 8 m 3 tracking volume, as is shown in Fig. 1.The flight software runs at 300 Hz on an onboard NVIDIA Jetson TX2 embedded computer.It includes the control algorithms (DFBC/NMPC and INDI inner-loop), and an extended Kalman filter fusing VICON (400 Hz) and IMU measurements (500 Hz) to obtain the full-state estimates.While DFBC runs at 300 Hz, the frequency of NMPC is limited to 100 Hz due to its computational complexity.The INDI inner-loop also runs at 300 Hz no matter which controller is being used.The onboard computer sends rotor speed commands of individual rotors to a low-level RadixFC board.The latter runs a customized firmware based on the open-source NuttX real-time operating system at 500 Hz to perform closed-loop rotor speed control, and sends rotor speed and IMU measurements to the onboard computer.We refer the users to Agilicious for more details about the hardware.VI.SIMULATION EXPERIMENTS In a set of controlled experiments in simulation, we examine DFBC and NMPC to address the following research questions: (i) Given the computational budget and a well-identified model, which control approach achieves better tracking performance?(ii) How do these two approaches perform in tracking dynamically infeasible trajectories, which inevitably lead to actuator saturation?(iii) How do these two approaches perform in simulated practical situations with a model mismatch, external disturbances, state estimation latency?Note that in this section, both methods are tested with the augmentation of the INDI inner-loop controller.In this section, we use NMPC and DFBC to represent those with INDI and an aerodynamic drag model for readability.

A. Reference Trajectories
Multiple elliptical reference trajectories on both horizontal and vertical planes are generated for testing the tracking performance.The results of tracking these simple trajectory primitives can further reflect the performance of tracking more complex 3D trajectories.These trajectories are parameterized by the maximum velocity V max , maximum acceleration a max , and ellipticity n.Specifically, the horizontal reference trajectories are expressed as: and the vertical reference trajectories are ξ r (t) = [r max sin(kt), 0, 5 + r min cos(kt)] T where The heading angle of each horizontal reference trajectory is selected to always point at the center of the ellipse; for vertical trajectories, the reference heading angles are kept constant.Specifically, 144 reference trajectories are generated by combining parameters a max ∈ {10, 20, 30, 40, 50, 60} m/s 2 , V max ∈ {5, 10, 15, 20} m/s, n ∈ {1, 2, 5}.
Tracking some of these trajectories may require the quadrotor to exceed its thrust limitations.These trajectories are referred to as dynamically infeasible.In many time-critical applications, such as autonomous drone racing, rotors need to reach their full thrust limits to fully exploit the agility of the platform and achieve faster lap times.Designing such an agile time-optimal trajectory requires a perfect model that captures real thrust limits and aerodynamic effects in high-speed flights.Without this model, the generated trajectory may exceed the quadrotor capabilities and becomes dynamically infeasible to accurately follow.Since such an accurate model is usually unattainable, studying the performance when tracking these infeasible trajectories is necessary.
For the simulated quadrotor with the configurations given in Table II, there are 68 infeasible trajectories and 76 feasible trajectories.They are determined by being tracked by a modified NMPC without thrust limits imposed.If the applied thrust exceeds the maximum thrust of the real drone, the trajectory is marked as infeasible.The feasibility of each reference trajectory is given in Table III

B. Evaluation Criteria
In the following comparisons, we use the root mean square error (RMSE) of position and heading as the precision metric of a method.The crash rate is another criterion to show the robustness of a method.A certain flight is defined as crashed if its position ξ violates the following spatial constraint at an arbitrary instant: We select b = [5, 5, 5] T meters in this study.

C. Tracking Dynamically Feasible Trajectories
First of all, we compare the performance of NMPC and DFBC in tracking 76 feasible trajectories.We perform the test in an ideal condition with perfect model knowledge and state estimates.Since only feasible trajectories are tracked, there is no saturation of single rotor thrusts and both methods succeed in tracking all trajectories without a single crash.
In this set of tests, we fine-tune the parameters of both methods (listed in Table.I) such that they achieve a similar position tracking error.Fig. 5 compares the boxplots of NMPC and DFBC when tracking trajectories with different reference maximum accelerations.Both methods have similar position RMSE in these flights.As for the heading tracking, NMPC has an average heading RMSE of 2.0 deg, which is better than 5.8 deg of DFBC.

D. Tracking Dynamically Infeasible Trajectories
1) Tracking Accuracy: Fig. 6 shows the box plot of position and heading RMSE of NMPC and DFBC in tracking the previously generated 68 infeasible trajectories.Crashed flights are excluded from this plot (see Fig. 7 for the crash rate).We find that the DFBC method shows smaller position RMSE in tracking trajectories with lower acceleration.However, as the reference acceleration grows, NMPC significantly outperforms 2) Crash Rates: Fig. 7 compares the crash rates of DFBC and NMPC in tracking all the infeasible trajectories.NMPC (solid-blue) outperforms DFBC (red-dash-dot), especially in tracking trajectories with high reference accelerations.
Here, we also demonstrate contributions of the tiltprioritized attitude controller and the quadratic-programming based allocation presented in Sec.IV.We replace them respectively with the conventional geometric attitude controller [24] and the direct-inversion allocation (30).The resultant crash rates are shown in Fig. 7 as well.It is obvious that the DFBC method tested in this study shows a significantly lower crash rate in tracking infeasible trajectories.

E. Robustness Study
While previous simulations are carried out in an ideal condition, this section studies the robustness of NMPC and DFBC methods in the following non-ideal conditions: • Model mismatch, including the center of gravity bias, mass mismatch, thrust and torque coefficients model mismatch, and error of aerodynamic drag model.• External disturbances.In details, we simulate constant external forces along x I , or external torques along both x B and y B .These disturbances last for 5 seconds.• Position, velocity, and attitude estimation latency.This aims at studying the effect of latency from such as signal transmission, state estimation algorithms, and sensors.As angular rates are measured directly by IMU with negligible latency, we only study the latency on pose and velocity estimates.These robustness studies are conducted using only the feasible trajectories since NMPC already shows its advantage over DFBC in tracking infeasible trajectories.We set the position RMSE as 5.0 m and heading RMSE as 90 deg for those crashed flights.
Table IV and Table V respectively present the results without and with INDI inner-loop.NMPC shows substantially higher robustness than DFBC against model uncertainties and disturbances on rotational dynamics due to CoG bias or external moments.Adding an INDI inner-loop controller can improve the robustness against these uncertainties.Even so, NMPC still slightly outperforms DFBC when both are augmented by an INDI inner-loop controller.
Unlike the rotational dynamics, NMPC shows a higher crash rate while experiencing model uncertainties and disturbances acting on the translational dynamics.For example, when the real mass is 30% higher, both controllers generate fewer thrusts  than required to track the trajectory, resulting in a constant position tracking error.With this tracking error, NMPC fails to converge in over 10% of flights and crashes the drone.The same reason also explains the significantly higher crash rate of NMPC in the presence of large external disturbances.Fig. 8 shows an example where NMPC failed in converging and crashed the drone after experiencing a 10 N external force disturbance.Note that adding INDI inner-loop does not improve the performance of either controller to reject disturbances in the translational dynamics.A previous research [13] implemented INDI in the position control (INDI outer-loop) for DFBC by virtue of its cascaded structure and improved its translational robustness.On the contrary, hybridizing INDI outer-loop with NMPC seems challenging owing to its noncascaded nature.
The two methods also perform differently in the presence of system latency.While NMPC slightly outperforms DFBC when the system latency is lower than 30 ms, as the latency grows, NMPC shows a much higher crash rate (68.0%) than DFBC (6.7%).As Fig. 9 shows, we repeat the tests under 30 ms and 50 ms latency with a high-gain setup and a lowgain setup.In both setups, these gains are tuned such that both controllers have identical average position RMSE in the ideal condition.As is expected, reducing gains will alleviate the effect of estimation latency.Interestingly, we observe that NMPC is more sensitive to the gain selections when system latency is larger than 30 ms.

F. Effect of Controller Parameters
Controller performance is dominated by the parameters or gains.However, finding optimal controller gains can be tedious and highly dependent on the drone parameters.Thus this section performs a sensitivity analysis to provide insights into the effect of controller parameters on the tracking performance.We set the control parameters given in Table I as the baseline, and compare the tracking performance with altered parameters with respect to the baseline.
1) Nonlinear MPC: Fig. 10 shows the effect of changing elements of the weight matrix Q.The position tracking performance is highly correlated with Q ξ and Q v .Among all the parameters, the position weight Q ξ plays the most important role.A larger penalty on the position error results in a smaller tracking error.However, as Q ξ continues growing, the NMPC starts to destabilize the drone.We believe that this is due to the presence of actuator dynamics in the simulation, which introduces system delay and causes instability with high gain controllers.Note that in real-world experiments, system latency can be larger due to state estimation algorithms and signal transmissions, which further reduces the permitted maximum Q ξ .Another observation is that increasing Q Ω Q Fig. 10: Position (top) and heading (bottom) tracking RMSE of the NMPC method in different control parameters.The vertical axis represents the ratio between RMSE and that from a baseline parameter setup.The horizontal axis represents the ratio between the control parameter and the baseline parameter.The RMSE of those crashed flights are clamped at 300% of RMSE base also adversarially affects the position tracking performance as it dilutes the power of position weights.Regarding heading control, as is expected, increasing Q q,z can improve the heading tracking performance while causing little effect on position tracking.
2) Differential Flatness Based Control: The controller gains of DFBC consists of position, attitude gains and weight matrix for control allocation.The gains are selected to make the position error ξ e = ξ − ξ r a second-order system: where ζ ξ and ω n,ξ are damping ratio and natural frequency, respectively.Then substituting ( 13) into (40) and replacing ξd with ξ (assuming a perfect inner-loop attitude control), we can calculate the position control gains as follows: Similarly, the attitude gains can be designed as: K Ω = 2ζ q diag( k q,x , k q,y , k q,z ).
where ζ q is the damping ratio, ω n,xy and ω n,z are desired natural frequencies of reduced attitude and yaw errors respectively.Fig. 11 presents the relationship between tracking error and closed-loop natural frequency.Both time constants and tracking errors are normalized by the baseline gains given in Table I.A higher frequency indicates a higher control gain, leading to less tracking error.However, high gains are also prone to suffer from system delay and subsequently cause instability.
Weighting matrix W is selected based on tuning in the simulation without sophisticated calculation.It aims at leaving An example of inappropriate selection of W is an identical matrix, which makes the QP-based allocation identical to a direct inversion.The performance in terms of crash rate has been compared in Fig. 7 necessity of using a well-tuned W in tracking dynamically infeasible trajectories.

VII. REAL-WORLD EXPERIMENTS
We conduct the experiments on a quadrotor in an instrumented tracking arena.A set of aggressive trajectories are executed to compare the closed-loop tracking performance of NMPC and DFBC in the presence of joint effects such as model mismatch, state estimation error, and system latency.The real-world experiments also highlight the contributions of the INDI low-level controller, the aerodynamic force model, both NMPC and DFBC methods.Experiments are conducted in tracking different reference trajectories, ranging from simple loop trajectories to FPV drone racing tracks (see Table .VI). Trajectories with the same 3D shape but different velocities and accelerations are also tested.The 3D paths of these trajectories are illustrated in Fig. 12.
We evaluated NMPC and DFBC, with and without the INDI inner-loop controller.NMPC with a PID low-level controller is also tested for comparison, implemented to track time-optimal trajectory in [4].In the following, these methods consider aerodynamic drag models by default, unless further mentioned.

A. Contribution of the INDI and Drag Model
Table VIII compares the average RMSE of both methods and shows the effect of the INDI inner-loop controller, and the effect of the aerodynamic drag model.For NMPC and DFBC, neglecting aerodynamic drag would increase position tracking error by 144% percent and 122%, respectively.
In contrast to the aerodynamic drag model, removing the INDI low-level controller influences more on the tracking accuracy.For NMPC, we see more than 364% and 115% increase in position and heading RMSE if INDI is not used.For DFBC, removing INDI leads to a more severe consequence: more than 705% and 492% of position and heading RMSE increase.In fact, DFBC without INDI cannot successfully track some of these trajectories without crashing the drone.These results indicate that, an adaptive/robust inner-loop controller plays a more important role than the drag model in accurately tracking aggressive trajectories.
We also make comparisons against the control method used in [4], which uses PID as the inner-loop controller of NMPC.This comparison is made in tracking a time-optimal trajectory generated off-line using the algorithm proposed in [4].In this task, the quadrotor needs to fly through multiple gates in a predefined order.Fig. 13 demonstrates the tracking results, including the 3D path and the position error.Clearly, both NMPC and DFBC augmented by INDI significantly outperform this NMPC with PID inner-loop controller.Comparative results on other reference trajectories can be found in Table VI.On average, the proposed NMPC with the INDI approach shows a position RMSE of 0.102 m, which is 70% lower than the position RMSE of NMPC with PID low-level controller (0.343 m).

B. Tracking Dynamically Infeasible Trajectories
In these real-world experiments, we also track three dynamically infeasible trajectories that exceed the maximum thrust of the tested quadrotor.Comparative results of NMPC and DFBC are given in Table VII.Clearly, NMPC shows significantly higher tracking accuracy than DFBC when tracking these dynamically infeasible trajectories, which is consistent with the simulation results in Sec.VI.Fig. 14 presents the 3D path tracking the Loop C reference trajectory.Tracking such a reference loop trajectory requires a higher collective thrust than the maximum thrust of the tested quadrotor.Hence the quadrotor cannot follow the reference trajectory and both control methods make the quadrotor take the shortcut to fly inside the reference trajectory.While DFBC results in a chaotic trajectory, NMPC makes the drone fly a much more regular loop trajectory with a smaller radius than the reference, thus requiring less collective thrust.The difference on the side-view is more distinctive: NMPC has much higher altitude tracking accuracy than the DFBC method.As is expected, NMPC requires significantly longer time to generate a single control command (2.7 ms) compared with DFBC (0.020 ms).Be that as it may, both methods can run onboard at sufficiently high frequency (≥100 Hz) to achieve accurate tracking of agile trajectories.

VIII. DISCUSSION
According to the simulation and real-world flight results, we conclude that both NMPC and DFBC, with an INDI inner-loop controller, can track highly aggressive trajectories with similar accuracy as long as the reference trajectories are dynamically feasible.
The advantage of NMPC appears when tracking dynamically infeasible trajectories that violate the rotor thrust constraints.Even though DFBC also uses the constrainedquadratic programming for control allocation, it only considers a single reference point; thus, the actions taken are too shortsighted to avoid future violations of these constraints.By    contrast, NMPC tracks the infeasible trajectory using future predictions, including multiple reference points that minimize the tracking error and respect single rotor constraints.Such difference helps NMPC outperform DFBC by 48% and 62% on position and heading accuracy, respectively, if the reference trajectories are dynamically infeasible.NMPC also demonstrates higher robustness against the rotational model mismatch.In real-world experiments, we also performed tests without INDI low-level controller for comparisons.As such, both methods suffered from model uncertainties on the rotational dynamics.In this condition, NMPC was still able to track all the trajectories, whereas DFBC experienced several crashes and had much higher tracking RMSE against NMPC.
However, NMPC also has limitations.For example, the biggest disadvantage of NMPC is that it requires significantly higher computational resources.On our tested hardware, the average solving time of the nonlinear NMPC is around 2.7 ms, while DFBC needs only 0.020 ms, which is around 100 times faster.This renders it impractical to run NMPC on some miniature aerial vehicles with a limited computational budget, such as Crazyflie [44].By contrast, we can deploy DFBC on light-weighted drones with low-end processors and conduct agile flights at speeds over 20 m/s, as long as the platform has enough thrust-to-weight ratio.
Another disadvantage of NMPC is that it potentially suffers from numerical convergence issues.Unlike DFBC, which has proof of stability or convergence of each sub-module, the nonlinear NMPC used in this comparison relies on the numerical convergence of the nonlinear optimization algorithm.Unfortunately, rigorous proof of its convergence is still an open question.In fact, nonlinear NMPC tends to fail in converging when the current position is too far from the reference, either caused by large external force disturbances or an error in the thrust-to-weight ratio model.For example, our robustness study shows a 10% higher crash rate of NMPC than DFBC when the real mass is 30% higher than the model.In addition, we also observe that NMPC is more prone to fail in converging than DFBC in the presence of large system latency.
In summary, NMPC is not superior to DFBC in all scenarios.If the reference trajectories are dynamically feasible, DFBC can achieve similar tracking performance while consuming only 1 percent of the control resource that NMPC requires.However, it is difficult to justify the feasibility of a trajectory near the physical limitations of the platform because of the model uncertainties such as aerodynamic effects.
Hence, NMPC excels at exploiting the full capability of an autonomous drone in a safer and more efficient manner, due to its advantage in handling the infeasibility of reference trajectories.Future work may further hybridize NMPC with differential flatness property to reduce the NMPC computational cost while retaining its advantage in handling dynamically infeasible trajectories.The reason why NMPC failed in converging also needs further investigation in order to deal with larger amount of model uncertainties.

IX. CONCLUSION
This work systematically compares NMPC and DFBC, two state-of-the-art quadrotor controllers for tracking agile trajectories.Simulation and real-world experiments are extensively performed to evaluate the performance of both methods in terms of accuracy, robustness, and computational efficiency.We report the advantages and disadvantages of both methods according to the results.This work also evaluates the effect of INDI inner-loop control on both methods.Real-world flight results at up to 20 m/s demonstrate the necessity of applying an INDI inner-loop on both NMPC and DFBC approaches.The effect of the high-speed aerodynamic drag model is also evaluated, which is found less influential compared to the inner-loop controller.

Fig. 1 :
Fig. 1: Top: Our quadrotor tracking a race trajectory.Bottom: Box plot comparing the position tracking root-mean-square-error (RMSE) of NMPC, DFBC, and their variations with INDI innerloop, with or without considering aerodynamic drag effects.For each method, data are collected from real-world flights tracking different reference trajectories with speeds up to 20 m/s (i.e., 72 km/h) and accelerations up to 5g.

Fig. 3 :Fig. 4 :
Fig. 3: The control diagram of the model predictive controller with an INDI inner-loop controller.

Fig. 8 :Fig. 9 :
Fig.8: Position tracking error for a loop trajectory (Vmax = 15m/s, amax = 40m/s 2 , n = 1).Grey area indicates the period with 10 N lateral disturbances acting on the drone.NMPC failed to converge and crashed the drone, while the DFBC method succeeds in recovering the drone after the external disturbance disappeared.

Fig. 11 :
Fig. 11: Position (top) and heading (bottom) tracking RMSE of the DFBC method in different control gain settings.The vertical axis represents the ratio between RMSE and the RMSE from a baseline gain setup.The horizontal axis represents the ratio between the controller gains and the baseline gains.RMSE of the crashed flights are clamped at 300% of RMSE base

Fig. 15
Fig. 15 and 16 respectively show the CPU time of DFBC and NMPC with INDI as the inner-loop controller while tracking the Race Track C trajectory.The computation time of DFBC is generally faster than 0.025 ms, which is significantly faster compared with NMPC that takes around 3 ms.In addition, both QP and INDI modules spend less than 0.01 ms CPU time.

Fig. 17 compares
Fig.17compares the average computational time of gen-

Fig. 13 :
Fig.13: Tracking performance of the Race Track C with a maximum speed up to 20 m/s.NMPC and DFBC with INDI inner-loop show similar tracking performance.Both outperform the state-of-the-art NMPC with a PID inner-loop[4]

Fig. 16 :
Fig. 16: CPU time of NMPC+INDI while tracking the Race Track C trajectory.Note that the vertical axis is presented in the log scale.

Fig. 17 :
Fig. 17: Processing time to generate each control command of NMPC and DFBC methods.

TABLE I :
Controller gains and parameters

TABLE IV :
Tracking performance comparisons in different types of non-ideal conditions commonly seen in practice without INDI inner loop.Each number shows the mean and standard deviation (SD) across tracking 76 dynamically feasible trajectories.

TABLE V :
Tracking performance comparisons in different types of non-ideal conditions commonly seen in practice with INDI inner loop.Each number shows the mean and standard deviation (SD) across tracking 76 dynamically feasible trajectories.
The position tracking RMSE and heading tracking RMSE are listed in Table.VI and Table.VII.

TABLE VI :
Tracking performance of different methods in real-world experiments

TABLE VII :
Performance in tracking three dynamically infeasible trajectories in real-world experiments.Position and heading RMSE comparisons between NMPC and DFBC with INDI inner-loop.

TABLE VIII :
Position and heading tracking RMSE of NMPC and DFBC methods.Adding INDI as the inner-loop improves the performance significantly for both methods.