Real-Time Nonlinear Control Allocation Framework for Vehicles with Highly Nonlinear Effectors Subject to Saturation

Hybrid Unmanned Aerial Vehicles UAV are vehicles capable of take-off and landing vertically like helicopters while maintaining the long-range efﬁciency of ﬁxed-wing aircraft. Unfortunately, due to their wing area, these vehicles are sensitive to wind gusts when hovering. One way to increase the hovering wind-rejection capabilities of hybrid UAV is through the addition of extra actuators capable of directing the thrust of the rotors. Nevertheless, the ability to control UAVs with many actuators is strictly related to how well the Control Allocation problem is solved. Generally, to reduce the problem complexity, conventional (CA) methods make use of linearized control effectiveness in order to optimize the inputs that achieve a certain control objective. We show that this simpliﬁcation can lead to oscillations if it is applied to thrust vectoring vehicles, with pronounced non-linear actuator effectiveness. When large control objectives are requested or actuators saturate, the linearized effectiveness based CA methods tend to compute a solution far away from the initial actuator state, invalidating the linearization. A potential solution could be to impose limits on the solution domain of the linearized CA algorithm. However, this solution only reduces the oscillations at the expense of a lag in the vehicle acceleration response. To overcome this limitation, we present a fully nonlinear CA method, which uses an Sequential Quadratic Programming (SQP) algorithm to solve the CA problem. The method is tested and implemented on a single board computer that computes the actuator solution in real time onboard a dual axis tilting rotor quad-plane. Flight test experiments conﬁrm the problem of severe oscillations in the linearized effectiveness CA algorithms and show how the only algorithm able to optimally solve the CA problem is the presented Nonlinear method.


Introduction
In the last twenty years, a boom in the (UAV) industry has been experienced all over the world. UAV prove to be flexible vehicles for several different tasks ranging from Mapping & this was proposed in [1], where the gust rejection capabilities of a conventional quad-plane were compared to the one of a dual-axis tilting rotor quad-plane, showing how effective the thrust vectoring capability is to minimize the influence of the wind on the vehicle in the hovering configuration.
Unfortunately, the large number of actuators increases the complexity of the Control Allocation (CA) problem for these vehicles. CA is the problem of distributing control effort over a set of actuators to achieve certain desired control forces and moments. These vehicles must solve the CA problem optimally to fully exploit their capabilities.
Several CA problem statements and resolution methods have been proposed over time. A good overview of the different CA methods is available in [2][3][4]. The unconstrained CA problem is most commonly solved with the pseudo inverse method, which inverts the mapping from control commands to the desired motion effects [5]. As for the constrained CA problem, redistributed pseudo-inverse [6] and direct allocation [7] are a few CA strategies based on the pseudoinverse method, able to optimize the computed solution in case of actuator saturation. More elaborate iterative constrained CA methods, formulate the problem as a constrained linear or quadratic programming optimization problem [8][9][10]. Note that all these problem formulations presented so far still share the same linearized control effectiveness assumption.
Regrettably, when vehicles with strong nonlinear behavior are employed, linearized effectiveness approaches tend to fail. This happens especially when a large control objective needs to be realized with ineffective non-linear actuators.
A much more effective method, potentially capable of dealing with such a problem, is the nonlinear programming method. The very first attempt to approach the CA problem with a nonlinear programming method was presented in [11], with the aim of solving the allocation problem for reentry vehicles. During reentry, these vehicles are exposed to very high angle of attack, where the nonlinear dynamics of the effectors are too evident that the CA problem cannot be properly solved with a linearized approach. Simulation results confirmed the weaknesses of the linearized control effectiveness approaches in solving the CA problem for these kinds of vehicles, highlighting the potential of the nonlinear programming method.
Following the same idea, the nonlinear resolution of the CA problem was also applied to overactuated ground vehicles [12][13][14], gaining more and more popularity over time. Recently, learning-based approaches for the resolution of the nonlinear CA problem were proposed with the aim of reducing the computational load [15,16]. Overall, although computationally intensive, the nonlinear programming approach for the resolution of the CA problem proved to be an extremely flexible and powerful method to determine control effectors solution for highly nonlinear and overactuated vehicles. However, so far, limited computational power of embedded computers, and the necessity of solving the problem in real time, prevented this CA method to be implemented and tested on a flying vehicle.
In this paper, we address the CA problem of vehicles with highly nonlinear control effectors, such as thrust vectoring UAV. Within the paper we initially show the ineffectiveness of state-of-the-art linearized effectiveness approaches for these type of vehicles. Then, we propose a method to solve in real time the nonlinear CA problem using a Sequential Quadratic Programming (SQP) algorithm. The proposed method uses the nonlinear Equations Of Motion (EOM) inversion to compute the actuator solution, thus overcoming the control effectiveness linearization limitation. The control framework used to implement such a CA method is a modified Incremental Nonlinear Dynamic Inversion (INDI), adapted to include the incremental law through nonlinear EOM evaluation.
To prove the applicability of such a control strategy on a real time and computationally constrained platform, the control framework was implemented and tested on the dual-axis tilting rotor quad-plane depicted in Fig. 1. The vehicle has a 2.4 kg take-off mass and features 8 servos and 4 motors for a total of 12 actuators controlling the 6 degrees of freedom. A commercial Raspberry Pi 4B Single Board Computer (SBC) was used to solve the Nonlinear CA problem in real time onboard the drone. Thanks to hardware and software optimization, it was possible to run the algorithm computing the actuator solution at an average refresh rate of 350 hz, while using just one of the four available cores of the SBC. Another of the available cores was dedicated to the serial communication with the main Flight Control Board (FCB), leaving a total of 2 free cores available for running algorithms supporting additional tasks.
To summarize, the research highlights of the work are: • We demonstrated the inapplicability of state-of-the-art linearized effectiveness Control Allocation algorithms on overactuated UAVs with thrust vectoring capability. • We proposed a Control Allocation algorithm based on Nonlinear EOM inversion able to adequately control overactuated UAVs with thrust vectoring capability. • We demonstrated the applicability of the proposed Nonlinear Control Allocation method on a flying vehicle by implementing the algorithm on a commercial SBC onboard of a dual-axis tilting rotor quad-plane.
The paper is structured as follows: In Section 2, the nonlinear system is presented and the INDI control framework is introduced. Within this section, the linear and quadratic programming formulation of the CA problem are also derived. In Section 3, the Nonlinear CA algorithm and the nonlinear INDI control framework are presented. In Section 4, the Nonlinear CA module is firstly implemented in simulation and then tested on the flying double axis tilting rotor quad-plane. In Section 5, some limitations of the Nonlinear CA methods are discussed. Finally, in Section 6 conclusions are drawn.

System Description and INDI Derivation
Considering the generic nonlinear system: Where x(t) ∈ R n is the state vector, u(t) ∈ R m is the control input vector and y(t) ∈ R l is the output vector. f nx1 : D x,u → R n is a nonlinear function representing the vehicle dynamics while h lx1 : D x → R l is a nonlinear function representing the output dynamics. For the derivation of the linearized CA methods, let's now consider the first order Taylor expansion of the first term in Eq. 1, centered around the current state x 0 and control input u 0 : where F x is the state dependent term and F u is the control input dependent term. Under the assumption of small controller sampling time and very fast actuators, it is possible to apply the time scale separation principle, letting us neglect the state dependent term x [17]. Based on the aforementioned assumptions, Eq. 2 can than be be rewritten as follows: From Eq. 3, it is now possible to derive the associated INDI control law for the vehicle. Assuming a direct expression between the state vector and the output vector (i.e. y = h(x) = x) and replacingẋ with the desired state derivativeẋ d , the control law becomes: where B † (x 0 , u 0 ) represents the generalized inverse of the effectiveness matrix, linearized around the current state x 0 and current control input u 0 .ẏ d −ẏ 0 is the output derivative errorẏ e , representing in this case the vehicle acceleration. The desired output derivativeẏ d is also known in literature as pseudo-control ν. The term u s represents the control input solution of the CA problem. Concerning the current output derivative valueẏ 0 , it represents the vehicle linear and angular accelerations and a sensor based estimation can be employed [18]. A scheme of the INDI control framework is depicted in Fig. 2.
For completeness, it is possible to manipulate the terms in Eq. 3 to directly compute the control input solution, rather than the control input increment: Then, by bringing the current actuator state on the right side and including it into the inversion law it becomes: where the desired acceleration vectorẏ d was replaced by the pseudo-control vector ν. Instead of minimizing the actuator solution norm, it is also possible to minimize the solution norm with respect to a desired control input state u d [3]: with The problem formulation in Eq. 7 is not yet a constrained problem and the solution can't be optimized in case of actuator saturation. Within the paper, we will refer to the CA strategy in Eq. 7 as the Pseudo Inverse Unconstrained (PIU) CA algorithm.

Quadratic Programming Formulation of the CA Problem
Equation 3 can be reformulated as quadratic programming problem [10]. This brings two main benefits. The first benefit is the optimization of the control input solution in case of actuator saturation. The second benefit is the possibility to include multiple objectives in the cost function to be minimized during the resolution process. Following the formulation presented by [18], the CA problem becomes: with where u s is the control input solution. The two diagonal matrices W u and W ν represent respectively the weight for the control input and the pseudo-control. The role of these matrices will be more clear later on, when we will apply the CA methods on the flying vehicle. The term γ ν is the scale factor assigning the preference on the two objectives included in the cost function. In order to prioritize the minimization of the control objective ν over the energy efficiency of the solution, this term should be chosen to be very high (γ ν >> 1). The solution of the least squares problem presented in Eq. 8 can be efficiently determined with the active set method presented in [18]. Within this paper, we will refer to this CA strategy as the Weighted Least Squares (WLS) CA algorithm.

Nonlinear Formulation of the CA Problem
So far, the methods presented to solve the CA problem associated with the INDI control framework make use of a linearized approximation of the actuator effectiveness. However, in some occasions, the vehicle effectors behavior cannot be fully represented by a linear system. For those types of platforms, it is possible to formulate and solve the CA problem with an iterative nonlinear optimization process [11]. If we still assume very fast actuators, the Nonlinear CA problem can be reformulated as follows: where γ u is the control input optimality scale factor. In contrast to Eq. 8, in this cost function, the scale factor is assigned to the control input rather than to the control objective. It is therefore suggested to assign a very small value (γ u << 1) to keep the problem priority on the primary control objective. The main noticeable difference between the WLS CA problem formulated in Eq. 8 and the Nonlinear CA problem formulated in Eq. 9 is the replacement of the linearized control effectiveness matrix with the nonlinear vehicle dynamics expression.
Concerning the two diagonal matrices W u and W ν , as in the WLS CA, they penalize respectively the associated actuator command and control objective. Within the paper, we will refer to this CA strategy as the Nonlinear CA algorithm.

Resolution Method of the Nonlinear CA Problem
The problem presented in Eq. 9 is a constrained nonlinear optimization problem whose solution can be determined using the nonlinear programming approach. Among the different nonlinear programming algorithms proposed in literature [19], SQP has been shown to be the most efficient and accurate algorithm in almost all circumstances [20]. We therefore decided to employ this algorithm to solve the nonlinear problem presented in Eq. 9.
Among all the different SQP implementations available, we decided to employ the line search based SQP algorithm implemented in the MATLAB fmincon function. A detailed derivation of the SQP method for the resolution of the Nonlinear CA problem is beyond the scope of this paper. For the interested reader, [19] explains in detail how the optimization process works and how it can be implemented.
In brief, the strategy of the SQP optimization process is to reduce the nonlinear problem complexity by the formulation of a series of unconstrained sequential quadratic programming sub problems. Each sub problem solution determines the searching direction while the step size is determined by evaluating a merit function containing both the cost function and the problem constraints. Specifically, two main characteristics of the chosen SQP optimization strategy are particularly beneficial for the implementation of our CA problem: • Strict feasibility of bounds: at every iteration, the implemented SQP algorithm never presents a solution that is out of bounds. • Robustness to NaN or Inf results: In case a step is taken which results in an invalid evaluation of the cost function, the algorithm tries to take a smaller step with the aim of computing a valid solution.

Implementation of the CA Methods on the Dual-Axis Tilting Rotor Quad-Plane
In this section we will derive the INDI control framework for the dual-axis tilting rotor quad-plane in Fig. 1. The vehicle was presented in our previous work [1] and has a total of 12 actuators including 8 servo-rotors. The disposition and notation of the actuators can be seen in Figs. 3 and 4. The number actuators gives the vehicle a full 6 Degrees Of Freedom (DOF) and makes it overactuated. Furthermore, the presence of tiltable rotors gives the aircraft highly non-linear actuator dynamics. All these characteristics give the aircraft a high potential maneuverability but at the same time make the CA problem particularly complex to solve. We will use the vehicle platform as a test bed to compare the different previously presented linear and non-linear CA methods (PIU, WLS and Nonlinear). The first step for the implementation of the presented CA methods on the flying vehicle is the derivation of the vehicle EOM. Once derived, the EOM can be replaced by the generic nonlinear system considered in Eq. 2.

Equation of Motion Derivation
For the sake of brevity, we defer to the appendix for the complete derivation of the equations of motion, the definition of the reference systems and the assumptions made. The equations of motion for the dual axis tilt rotor quad-plane are as follows: whereP e are the linear acceleration in the earth reference frame andω represent the body rates derivative. The forces vectors F a and F p represent respectively the aerodynamics forces and the thrust produced by the propellers, expressed in the earth reference frame. As for the moments equilibrium, the moment term M t represents the torque generated by the rotors due to the propeller thrust, M d represents the torque generated by the rotors due to the propeller drag, M i represents the torque induced by the i-th propeller inertia due to the motor rotational rate change˙ i . The term M p represents the torque generated by the rotor precession term due to the tilting rotation, M r models the inertial term of the rotors due to the vehicle rates and M tilt represents the torque generated by the rotor inertia due to the tilting rotation. A more detailed explanation of each element of Eq. 10 is covered in the Appendix.

INDI Framework for the Dual-Axis Tilting Rotor Quad-Plane
For a vehicle such as the dual-axis tilting rotor quad-plane, the control module scheme is slightly different from the one of a conventional under-actuated vehicle. In the hovering configuration, the dual-axis tilting rotor quad-plane has full 6 DOF control. Therefore, in this flight configuration, an outer loop is not required and a single control loop can be employed to directly control attitude and position. The resulting INDI control scheme is presented in Fig. 5, where a linear error controller is employed for the generation of the desired linear and angular acceleration. A representation of the error controller scheme is depicted in Fig. 6.
As for the pseudo-control vector, it is composed of the following desired linear and angular accelerations, expressed in the earth reference frame e : For the control input vector, it is constituted of a total of 12 elements containing both the motor rotational speed and the rotor tilting angles:

Derivation of the Effectiveness Matrix for the Linearized CA Methods
For the PIU and the WLS control allocation methods, an estimate of the linear actuator effectiveness of the vehicle is required. As described in the previous section, the estimation of the linear control effectiveness matrix will be evaluated in the current state x 0 and current control input u 0 and then inverted to determine the control input set. The effectiveness matrix is composed of the partial derivative of the EOM in Eq. 10 with respect to the control input array defined in Eq. 12, for a total matrix size of 6x12: For the determination of the solution, both the WLS and the PIU have to invert the constant effectiveness matrix in Eq. 13.
Unfortunately, for some current actuator values and states, the effectiveness matrix may contain small elements. This means that the associated actuator solution will be far away from the current actuator state, invalidating the linearization. The PIU and WLS CA algorithms have no information about the system's nonlinearities and therefore have to wait for the actuator movement to realize that the computed solution does not give the desired result. This behavior is clearly visible by running the WLS CA algorithm with a desired vertical acceleration and all motors slightly tilted outward, as depicted in Fig. 7: Due to the linearization, the effectiveness matrix calculated with this actuator state contains some non-zero elements for the terms ∂z ∂ g i . This means that the CA algorithm may use the lateral tilting to control the vertical acceleration without increasing the RPM of the motors, which would lead to an increase in the cost according to Eq. 8. Unfortunately, the effectiveness matrix only contains a coefficient and does not have any information of the nonlinearities of the platform actuators. This means that the WLS CA algorithm will Fig. 7 Initial actuator state defined in Eq. 14 eventually compute a big lateral tilting angular change as an alternative to increasing the motor power. The actuator solution computed by the WLS CA algorithm, for the problem conditions in Eq. 14 is reported in Eq. 15 and represented in Fig. 8.
By comparing the computed solution with the saturation points of the actuators in Table 2, it is visible that the proposed WLS actuator solution saturates both the motor rotational speed and the lateral tilting angle in an attempt to generate a vertical acceleration. This condition is clearly non-optimal, as a zero tilting angle would yield more vertical acceleration, and the controller will detect the ineffectiveness of the solution only once the tilting actuator will move. This will eventually trigger permanent oscillations in the angular solution computed by the WLS CA algorithm.
For the same scenario, the PIU CA algorithm computes an unfeasible actuator solution, with a motor rotational speed 1.
In this case, the CA algorithm only saturates the motors, with both azimuth and elevation tilting angles equal to zero.

Simplification of the EOM for the Nonlinear CA Method Formulation
To decrease the computational load of the Nonlinear CA method, only the main terms of the EOM of Eq. 10 are considered. The vehicle dynamics considered for the nonlinear optimization process in Eq. 9 is the following: The secondary terms of the EOM not considered for the inversion will be treated as unmodeled dynamics and handled by the error controller.
Note that this formulation of the EOM also considers the aerodynamic effects not considered in the linearized actuator effectiveness matrix formulated in Eq. 13. It is now possible to insert the simplified EOM expression of Eq. 18 into the cost function derived in Eq. 9 to complete the Nonlinear CA problem formulation.

CA Parameters Choice
For the implementation of the previously introduced CA methods, a few parameters have to be set and described. The parameter in common between all the CA algorithms is the choice of the desired actuator state. For the tilting system, during the hovering test, the desired state of the tilting angles is set to the vertical condition (b i = g i = 0). Regarding the motor desired state, it was set to 100 rad/s, equivalent to the minimum applicable value.
For the WLS and Nonlinear CA algorithms, the weighting matrices W u , W ν and the penalty parameters γ ν and γ u have to be chosen. Starting with the control input weighting matrix W u , it is a square diagonal matrix of the same size as the number of actuators (12x12 in that case). Each diagonal element of the matrix assigns a penalty on the usage of the respective actuator and it is located on the secondary objective of the cost function. Ideally, we prefer to prioritize the usage of the tilting system over the motor for a more energy efficient solution. Therefore, we decided to apply a unitary penalty coefficient to the tilting system and a penalty coefficient equal to 3 to the motors. The weighting matrix W u used for the vehicle then becomes: As for the pseudo-control weighting matrix W ν , as already explained in [18], it should be set in such a way for the problem to firstly prioritize the pitch and roll changes, then the linear acceleration changes and the yaw angle. This leads to a pseudo-control weighting matrix W ν of this form: Concerning the penalty parameter γ , we set the value of 1e5 to γ ν and 1e-5 to γ u while for both the Nonlinear and WLS problem formulations, a maximum number of 60 iterations was set and the current actuator state u 0 was used as starting point for the optimization process. The last parameter that needs to be discussed for the constrained CA problems is the constraint definition. Ideally, the problem constraints are constant at every iteration and equal to the physical actuator limits. However, one may claim that with a limit on the maximum change in control input, the error that would be made by linearizing the system would become smaller. To investigate this, when using the WLS CA algorithm for the resolution of the CA problem, a different choice of the problem constraints will also be tested. The modified constraint, denoted as "local tilting constraint" later in the paper, is of the following form: Where the superscripts "0" indicates the current value, the superscripts max and min indicate respectively the maximum and minimum physical value of the actuator and δ c represents the constant angular region of constraint. The maximum and minimum actuator values for the dual-axis tilting rotor quadplane are reported in Table 2.

Actuator Model Identification
The previously proposed CA methods all rely on an incremental control law and therefore require an estimate of the current actuator state to work. Accurate modeling of the actuator dynamics is also needed for the vehicle simulation. For the motor propeller actuation system, we determined its dynamics through bench tests using an external RPM sensor. The result of a step test campaign on the motor propeller system is shown in Fig. 9A. In the same plot, the motor rotational speed evolution was also fitted with a first and a Step test actuator response for the identification of the dynamics. The left plot (A) is the motor step test, the central plot (B) is the azimuth tilting angle step test and the right plot (C) is the elevation tilting angle step test. The blue curve represents the desired step change and the red curve reports the variable evolution. The purple and yellow curves are a fit of the variable evolution using respectively a first and a second order discrete transfer function second order discrete transfer function. As visible, the motorpropeller system can be modeled using a first order discrete transfer function. The analytical expression of the first order discrete transfer function associated with the motor dynamics is the following: The sampling frequency used to determine the transfer function in Eq. 22 is 500 Hz. This is the update frequency of the main flight computer running onboard the vehicle.
Concerning the identification of the tilting dynamics, the approach of [1] was used. An IMU was mounted solidly to the tilting structure and angle step changes were induced on the tilting servos. The dynamics is then recorded and analyzed to determine the tilting rotor dynamics at different motor speeds.
The motor spinning power chosen to identify the motor dynamics discrete transfer function was 75%. The tilting evolution for the azimuth and elevation rotor angle step test is displayed in Fig. 9B and C. The tilting angle evolution was fitted with a first and second order discrete transfer function featuring also a rate saturation limit. In this case, the second order transfer function gives a slightly better representa-tion of the tilting dynamics. The analytical expressions of second order discrete transfer function used for the azimuth and elevation tilting angle dynamics are the following: Where R l represents the rate saturation limit value of the tilting system. The same sampling frequency of 500 Hz was used for the determination of the discrete transfer functions in Eqs. 23 and 24. For completeness, in Table 1, the characteristics of the equivalent continuous time system associated to the actuators are displayed.

Actuator Constraint and Problem Normalization
For a complete CA problem formulation, two more steps are required. The first step is the determination of the physical  limits of the actuators on the vehicle. The second step is to use the total travel of the actuators to normalize the CA problem, in such a way to equally compare the actuator effects during the CA resolution process. The physical limits of the actuators are displayed in Table 2.
With the maximum and minimum actuator limits determined, it is possible to scale the control input vector presented in Eq. 12 for the normalization of the CA problem: Where the scaling factors G m , G el and G az are displayed in Table 3. Using the normalized control input array u n presented in Eq. 25, the overall actuator travel will be the same for all the actuators and will be equal to 2. For a complete normalization of the CA problem, the maximum, minimum and desired actuator value, the effectiveness matrix B and the EOM were normalized as well.

Experimental Setup
The implementation of the previously presented CA methods on a real-time setup requires a reasonable amount of computational power that a microcontroller alone cannot provide. In order to solve the computational deficit problem, we decided to distribute the computational load among two boards. The main FCB used is a Pixhawk 4 flight controller running the open-source paparazzi software 1 [21]. The task of the Pixhawk 4 is to interact with the vehicle sensors and instruments, control the actuators and to generate the pseudocontrol vector. Once the pseudo-control vector is generated, it is communicated over a UART interface to the Raspberry Pi 4B, 2 whose task it is to solve the CA problem and communicate the actuator solution back to the Pixhawk 4. The development of the CA modules operating on the SBC was characterized by a few steps. Initially, the PIU, WLS and Nonlinear CA algorithms were implemented and tested in Simulink. Secondly, using the Simulink Coder support package, 3 the C functions implementing the CA methods were generated. Finally, the generated C code of the CA routines were manually optimized, compiled and then deployed on the Raspberry Pi 4B SBC. For the communication between the Raspberry Pi 4 and the Pixhawk, a serial communication was employed. In order to maximize the computational speed of the CA solver, the serial communication on the raspberry Pi was implemented on a different thread. As for the rest of the components employed for the vehicle hardware such as the type of tilting servos and motors used, the reader can refer to [1].

Simulation Results
With the knowledge of the vehicle dynamics, the actuator characteristics, and a proper identification of the control framework, it was possible to set up a simulation environment in Simulink. The simulation outcome could then be compared to the flight test result for the model validation.
A mixed maneuver engaging attitude, altitude and position changes was chosen. Initially, a pitch and roll angle of 20 degrees are commanded. Then, the vehicle is asked to track a desired altitude and position step change.
As a reference for the model validation, only the Nonlinear CA algorithm was used to compare the simulation outcome with the flight test results. In Fig. 10 the attitude and position evolution between flight test and simulation is compared. As for the actuator evolution, for the sake of simplicity, only the tilting angles and motor power of the rotor number 2 are compared in Fig. 11.

Comments
The flight test results and the simulation outcome are matching, thus validating the modeling of the vehicle and actuator dynamics. The only visible mismatch happens around time 20 seconds, during the climbing phase of the maneuver. As also confirmed by the logs, this mismatch may have been caused by a sudden voltage drop related to a spike in power demand of the motors. This condition also influences the motor rotational speed response of the flight test which presents a damped oscillation not present in the simulation.

Flight Test Results
To assess the behavior of the different CA methods on the flying vehicle, the same maneuver proposed for the simulation was programmed on board the drone. The flight test was carried out in the TUDelft Cyberzoo flight arena and the flying vehicle was always secured with a rope from the arena ceiling. The rope was manually handled in such a way to minimize the interference with the flight.
For an accurate estimation of the vehicle position and orientation, an Inertial Navigation System (INS) composed of an OptiTrack optical motion system, a 3-axis accelerometer, a 3-axis magnetometer and a 3-axis gyroscope was employed. Since the controller is designed to track accelerations in the earth reference frame e , an estimation of the inertial accelerations of the vehicle in the earth reference frame is needed. This estimation is obtained by projecting the inertial body accelerations recorded by the accelerometer to the earth reference frame through Eq. 26.
An overview of the different CA algorithms tested on the flying vehicle are summarized in Table 4 while in Fig. 12, the attitude and position tracking for the successful and partially successful flights are displayed. The partially successful flight curves are truncated at the time the flight test was stopped. In Fig. 13, the actuators evolution of the rotor number 2 are shown. As done for the analysis of the simulation results, the rotor number 2 was chosen as reference for the analysis of the solution computed by the different CA algorithms during the flight tests. Finally, in Fig. 14, the plot containing the run-time of the three main CA algorithms is displayed. The run-time was defined as the time needed by the CA algorithm running on the Raspberry pi to compute the actuator solution. This value defines the actual update period of the desired actuator value and it is generally not constant.

Comments
Several interesting observations can be made by looking at the flight test results for the tested CA methods. Both the Nonlinear and WLS CA methods provide a valid actuator solution at every time step. These methods are constrained optimization processes and therefore always provide an actuator solution within the actuator limits. On the other hand,  During that rapid maneuver, the PIU CA algorithm computes a motor speed solution way above the saturation point. This means that the controller tries to achieve the desired acceleration set by clipping a value to an already saturated actuator state, leading to a visible loss of attitude tracking.
Another major consideration has to be done when applying equal weight in the control input matrix W u , as done for the PIU CA algorithm and the WLS CA algorithm with unitary W u . With those two CA methods, there is no penalty on the usage of the motor over the tilting system. This means that it is mathematically more convenient to use motor power differential instead of rotor tilt to achieve yaw changes. This condition is observed around time 15 seconds, where the vehicle corrects a yaw error by decreasing the motor power of rotor 2 and 3 while increasing the motor power of rotors 1 and 4. Unfortunately, this practice leads to the saturation of rotors 1 and 4, therefore reducing the maneuverability of the vehicle.
Being able to prioritize the usage of the tilting system over motor power is a very important feature for the determination of an energy efficient solution. However, while it was possible to fly the vehicle using the Nonlinear CA method with the control input weighting matrix W u in Eq. 19, it was impossible to do so with the WLS CA algorithm. Assigning penalty on the motor power highlights the limitations of the linearized approximation of the actuator effectiveness. Under this condition, the WLS CA method tends to allocate the control objective to the tilting system as much as possible, leading to rapidly oscillating tilting angles solution coupled with an inaccurate motor power solution. To better understand the reasoning behind the tilting oscillations, the reader It is possible to partially bound the oscillatory actuator solution by applying a local constraint to the tilting solution in the WLS CA problem formulation, as presented in Eq. 21. However, this condition induces a delay in the acceleration tracking proportional to the speed of the actuators. The delay in the acceleration tracking leads to an overshoot in the position and attitude tracking of the vehicle. When a δ c = 5 deg is applied, the overshoot in the position tracking becomes so pronounced that the flight test had to be stopped prematurely to avoid a crash with the lateral net.
Concerning the proposed Nonlinear CA method, it always computes a smooth, valid and effective actuator solution, correctly prioritizing the desired control input and actuator use according to the provided weighting matrices W ν and W u . Among the three analyzed algorithms, the proposed Nonlinear CA is the Control Allocation strategy providing the best tracking for both position and attitude, as visible from Fig. 12.
As for the computational load of the different algorithms, the run-time needed by the algorithms to compute the actuator solution is displayed in Fig. 14. The WLS CA algorithm is the most computational-intensive algorithm with an average of 4.7 milliseconds needed for the determination of the actuator solution. The large run-time needed by the WLS algorithm for the determination of the solution is mainly associated with the frequent actuator saturation condition experienced during the flight. In case of actuator saturation, the iterative activeset optimization process run by the WLS algorithm needs multiple iterations to determine the solution, leading to an increase in execution time. The second most computation-ally intensive algorithm is the Nonlinear CA algorithm, with an average run-time of 2.9 milliseconds. Interestingly, the Nonlinear CA algorithm, as opposed to the WLS CA, experiences a reduced run-time when a set of actuators is saturated. Finally, the computational time of the PIU CA algorithm is the smallest one, requiring an average of 11 microseconds to compute, and it is more or less constant during the whole flight.

Limitations of the Proposed Nonlinear CA Method
There are two main limitations concerning the applicability of the proposed Nonlinear CA on flying vehicles. The first limitation has to do with the computational power required by the Nonlinear CA algorithm to run. Depending on the number of actuators and the complexity of the EOM, the cost function formulation and its gradient might increase in complexity, leading to a longer run-time needed by the SQP algorithm for the computation of the actuator solution. A longer run-time will reduce the phase margin of the controller if it becomes excessive. A similar issue can be encountered if a less powerful SBC is used to run the Nonlinear CA algorithm. Note that the nonlinear CA algorithm determines the actuator solution iteratively. Therefore, there is a direct relation between the number of iterations and the accuracy of the computed actuator solution. At the same time, the number of iterations also affects the overall run-time and computational load of the algorithm. The choice of the number of iterations is therefore a key parameter that must be chosen as a compromise between the computational load and the accuracy of the solution. The second limitation affecting the Nonlinear optimization process is to ensure to find an actuator solution which is not in a local minimum of the cost function. The easiest way to prove the absence of local minima in the cost function, is to ensure global convexity of the function. Unfortunately, the complexity of the proposed cost function makes it difficult to analytically prove the global convexity. Even if the convexity of the cost function cannot be analytically determined, it is still possible to investigate if starting the optimization process using a different initial condition would change the cost function's final value. This gives a rough idea about the shape of the cost function and the presence of local minima.
In order to check if the cost function contains local minima, a monte-carlo simulation was performed varying the initial actuator set of the optimizer u init , for different acceleration set, actuator condition and vehicle states. A total of 7000 tests were performed each time using 300 different initial actuator sets, for a total of more than 2 million optimization runs. A uniform statistical distribution was used to determine states and control objective values at every itera-tion. The constant and randomized variables used during the statistical analysis are displayed in Table 5. The choice of the variable ranges was set in such a way as to reproduce extreme conditions for the vehicle state, actuator state and for the control objective.

Statistical Analysis Results
Out of the 7000 tests performed during the simulation, 98.6% of the time, the cost function value determined using u init = u 0 is within 10% of the minimum cost function value determined over the 300 runs with randomized u init . Only 93 tests out of 7000 provided a final cost function mismatch of more than 10%, highlighting a potential local minimum condition.
In order to evaluate the primary objective term of the cost function, the residual norm of the solution determined using a different initial actuator set was compared with the solution providing the minimum residual norm for that test. The residual vector is defined as the difference between the desired accelerations and the accelerations achieved by the determined actuator solution. The norm of this vector will therefore be composed by a sum of angular and linear accelerations with units rad/s 2 and m/s 2 .
In Fig. 15, the residual norm of the "minimum norm actuator solution" is compared with the solution computed using different u init . We define as "minimum norm actuator solution", the actuator solution providing minimum residual norm out of the 300 randomized u init runs. The blue histogram in Fig. 15 represents the residual norm error between the "minimum norm actuator solution" and the "maximum norm actuator solution" of the 300 optimization runs. In the same figure, the red histogram is the residual error between the "minimum norm actuator solution" and the solution coming from an arbitrary initial actuator set. The arbitrary actuator set was chosen as the value coming from the first of the 300 random optimization runs. Finally, the red histogram shows the residual norm error between the "minimum norm actuator solution" and the solution computed using the current actuator state as initial point for the optimizer.
Out of the 7000 tests performed, the maximum residual norm between the "minimum norm actuator solution" and the solution computed using u init = u 0 is 3.5. This value increases to 7.2 when the arbitrary u init is considered and reaches 10.4 when it is compared with the "maximum norm actuator solution". Interestingly, it was observed that in a total of 20 tests, the solution computed using u init = u 0 leads to the lowest residual solution.
To evaluate the secondary objective term of the cost function, the actuator solution displacement from the current actuator state was also analyzed. A bigger actuator displacement in the solution involves a rapid and sudden change of the actuator state. In certain circumstances we would like to limit the actuator displacement in favor of a slightly bigger residual to avoid sudden jumps in the actuator solution and to achieve a smooth actuator evolution over time. This is particularly important in our case, where the actuators are mainly constituted of servo motors. To evaluate this parameter, the scaled actuator displacement of the solution computed using different u init conditions is shown in Fig. 16. In this figure, the magenta histogram shows the scaled solution displacement from u 0 of the "minimum norm solution", the red histogram represents the scaled solution displacement using the arbitrary initial actuator set while the green histogram is the scaled solution displacement when using the current actuator state as initial value for the optimizer.
Starting the optimization from different initial actuator states also provides different actuator solutions and therefore different residuals. However, the initial actuator state leading to a solution with minimal residual norm is not known in  Fig. 15 clearly shows that starting the optimizer from the current actuator state provides a lower residual norm than starting it from an arbitrary actuator point.
At the same time, it is clear from Fig. 16 that initializing the Nonlinear CA algorithm from an initial actuator set different from the current one, leads to a solution with higher displacement from u 0 . Therefore, even if by starting the optimization process at the current actuator state we might not be able to find the solution with the lowest norm, we expect to find the sub-optimal solution with minimum actuator dis- Fig. 16 Scaled displacement of the computed actuator solution from the current actuator value for different u init conditions. Notice that the y-axis is logarithmic. The displacement is defined as the norm of the difference in normalized control input u n − u n0 , with u n as defined in Eq. 25 placement. This choice should mitigate sudden jumps in the actuator solution, guaranteeing a smooth actuator evolution over time.

Conclusions
In this paper we addressed the Control Allocation (CA) problem of hybrid Unmanned Aerial Vehicle (UAV) with nonlinear effectiveness actuators, such as tilting rotor vehicles.
Through flight tests and simulations we proved that stateof-the-art linearized effectiveness CA algorithms like the Pseudo Inverse Unconstrained (PIU) and Weighted Least Squares (WLS) are not capable of determining a smooth and effective actuator solution for tilting rotor vehicles. We proposed a Nonlinear CA algorithm capable of computing an optimal and smooth actuator solution for such vehicles. Furthermore, the vehicle performance in terms of attitude and position tracking obtained using the Nonlinear CA algorithm outperformed the one obtained using any other tested CA method. It was observed that under high control objective commands or when an uneven control input weighting matrix W u is used, the solution computed by linearized effectiveness algorithms is characterized by uncontrolled and wide oscillatory tilting angles, coupled with an inaccurate motor power solution. In an attempt to mitigate the oscillatory behavior of the WLS CA algorithm, a local tilting constraint was included in the WLS problem structure. Although potentially beneficial, this condition only bounds the oscillations of the tilting solution at the expense of an added acceleration response delay, inversely proportional to the applied angular region of constraint δ c .
Concerning possible future development, even if no criticalities were found during the flight tests, a statistical analysis on the cost function used in the Nonlinear CA algorithm revealed the presence of some local minimum points. Future work should focus on this aspect, trying to modify the cost function in such a way to eliminate the presence of local minima in it. Furthermore, the Nonlinear CA algorithm could be extended to the forward flight, with the aim of developing a unified control strategy for hovering, transitioning and forward flight of tilt rotor hybrid vehicles.

Acknowledgements
The work was carried out within the Unmanned Valley Project. The authors would like to thank the "Europees Fonds voor Regionale Ontwikkeling(EFRO)" who is founding the Unmanned Valley project under grant code KvW-00168 for the South-Holland region.
Author Contributions all authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Alessandro Mancinelli. The first draft of the manuscript was written by Alessandro Mancinelli and Ewoud J.J. Smeur and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Data Availability the datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations
Ethics Approval the work does not require ethics approval since it does not involve human or animal subjects. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix: Mathematical model of the dualaxis tilting rotor quad-plane
In this section, we will report the EOM derivation for the dual-axis tilting rotor quad-plane as described in our previous paper [1].

A.1 Reference Frames and Notation
Firstly, the rotor disposition and spinning direction have to be characterized. In Fig. 3, a scheme containing the rotor disposition and the motors spinning direction is shown.
Secondly, it is important to define the different reference frames used for the characterization of the vehicle dynamics: • Earth Frame e : Origin on the Earth surface, x e aligned with Earth north, y e axis aligned with Earth east and z e axis pointing towards the center of Earth. • Wind Frame w : Origin in the airplane CG, x w axis in the vehicle plane of symmetry pointing in the wind speed direction, z w in the vehicle plane of symmetry and perpendicular to x w , y w perpendicular to x w and z w , pointing to the right wing.
An overview of the Earth, Body and Propeller frames and the identification of the rotor tilting angles is shown in Fig. 4. The transformation matrices for the projection of a vector from one reference frame to another can be identified as follows: For the coordinate transformation between the body reference frame b to earth reference frame e the following matrix is used: where c and s represent the abbreviation respectively of the cosine and sine function, while φ, θ and ψ are the Euler angles in the traditional ZYX order. For the coordinate transformation between the propeller frame p to body reference frame b the following matrix is used: where the angles b i and g i are the i-th rotor tilting angles. Conventionally, within the paper we will refer to b i as ele-vation tilting angle and to g i as azimuth tilting angle. For a visual representation of the tilting angles, the reader can refer to Fig. 4. Concerning the coordinate transformation between the wind frame w and the body frame b , the following matrix is used: where α is the angle of attack and β is the sideslip angle. Finally, we define the matrix T used to obtain the rate of change of the Euler angles from the body rates ω: where represents the Euler angle vector composed by φ, θ and ψ.

A.2 Assumptions
In order to facilitate the EOM derivation, a few assumptions are made: • Inflow into the propeller is assumed not to influence it's performance. • The thrust generated by the rotor is always perpendicular to the propeller disk and it is applied in the center of the propeller disk. • The change in the body inertia due to the rotor tilting is negligible and x b , y b and z b are vehicle principal axes. • x p , y p and z p are principal axes for the propeller, and the inertia terms I p x x and I p yy are negligible. • The inertia tensor of the tilting mechanism in the propeller reference frame p is a diagonal matrix.

A.3 Equations Of Motion Derivation
With the reference frames and assumptions defined, it is possible to analyze all the forces and moments contributing to the system dynamics for the development of the EOM: whereP e are the linear acceleration in the earth reference frame andω represent the body rates derivative.