Feedback linearization control of quadrotor with tiltable rotors under wind gusts

Article history: Received 6 June 2017 Received in revised form 15 August 2017 Accepted 10 September 2017 Quadrotors are popular unmanned aerial vehicles that have a plethora of applications for civilian and military purposes. This is due to their superior agility and maneuverability which widens the span of applications. In this paper, a feedback linearization controller is designed to control all the states of the over actuated quadrotor with tilting rotors that was developed by the author and his colleagues. The controller is introduced in a novel approach that overcomes the problem of nonlinear inputs and that decouples the system into completely two independent subsystems while rejecting wind gusts. In addition, an optimization algorithm is introduced to choose among the possible sets of inputs based on energy consumption minimization. The results demonstrate that the quadrotor with tilted rotor can effectively attain the desired trajectory in the presence of wind disturbance.


Introduction
*Quadrotors have been extensively used for various civilian and military purposes.Examples of such applications are surveillances, traffic monitoring, rescue missions, patrolling forests in case of fire outbreak, warfare, and other risky missions.Modeling for the conventional quadrotor exists in numerous literatures for instance in Bouabdallah et al. (2007), Pounds et al. (2010), and Voos (2009).These model structures possess limitations regarding the orientation coupling, because of insufficient inputs, hence making it incongruous for some particular tasks.In light of this, different modeling schemes have been employed Ryll et al. (2012) and Oner et al. (2008) in quadrotors with tilt wing to address some deficiencies of the conventional quadrotor to improve the actuation capacity.Additional inputs have only been successfully aimed at decoupling the quadrotor orientation either to achieve forward flight mode or moving sideways while hovering as the case may be.This is particularly useful in applications requiring more maneuverability.
However, a thrust tilting approached Elfeky et al. (2013Elfeky et al. ( , 2014) ) and Hua et al. (2013) were applied to a conventional quadrotor to decouple the orientation, with the aim of reducing additional input.Each rotor is designed in a way such that it can tilt around two axes with respect to fixed body frame.In this design, the number of inputs becomes twelve instead of four in the conventional quadrotors.Therefore, it is enough to have six inputs to fully actuate this modified quadrotor.Also twelve inputs may be used to impose arbitrary trajectories and other requirements such as disturbance rejection.With this design, each of the twelve states (outputs) (6 positions/orientations -6 translational/rotational speeds) can be controlled independently and freely.In addition, the novelty of this system is that, it can perform all the desired control objectives with half of its actuators faulty.Various control techniques have been used successfully to control conventional quadrotor within the range of their capabilities.In Bouabdallah et al. (2007), Pounds et al. (2010), and Ryll et al. (2012), linear control laws have been extensively employed.However in Voos (2009), Altug et al. (2005), Senkul and Altug (2014), Morel and Leonessa (2006), Liu et al. (2013), Li and Wang (2013) and Saif et al. (2012) non-linear control methods were also successfully employed.For the quadrotor with tilt-wing different control methods have been applied for instance, an output feedback control in Ryll et al. (2012), LQR control in Oner et al. (2008) and PID control in Salih et al. (2010) and Elfeky et al. (2013) to achieve different flight modes.
In this work, we shall exploit the capability of quadrotor with tilted rotors to tackle the issues of orientation decoupling.Each rotor possesses two degrees of freedom, thus increasing the number of input to twelve.Although, six inputs are sufficient for full actuation, additional inputs can be utilized for different flight modes and wind gust rejection.Linearized feedback controllers are used to control the quadrotor.This paper is organized as follows.Section 2 presents the dynamic model of quadrotor with two degree of freedoms tilting rotors.The wind gust modelling is stated in Section 3. The design of feedback controllers based on feedback linearization is given in Section 4. Section 5, discusses the optimization techniques that are used to select 3 inputs out of 12 such that an objective function related to the power consumption is minimized.Simulation results that show the effectiveness of the designed controller are shown in Section 6.Finally, the paper conclusion is given in section 7.

Frames and rotation matrices
The quadrotor consist of five rigid bodies that are connected together and are in relative motion around them.These five bodies are the quadrotor body itself Band four propellers P i attached to the body.Denote the world inertial frame by be the quadrotor body frame attached to its center of gravity.In addition, the frames of the rotors are taken to be parallel to each other and parallel to the quadrotor body frame and are given by F P i ∶ {O P i , X P i , Y P i , Z P i }, i = 1, … ,4.
The orientation of each of the rotors is controlled by two rotations with respect to the rotor-fixed frame; α i , a rotation about Y P i , and β i , about Z P i .This rotation creates another rotating frame for the rotors denoted by F P ̅ i ∶ {O P ̅ i , X P ̅ i , Y P ̅ i , Z P ̅ i }, i = 1, … ,4 as shown in Fig. 1.
Unlike the conventional quadrotor, in this proposed design (Elfeky et al., 2013(Elfeky et al., , 2014)), when the rotors are aligned alongZ P i , say, rotor 1 and rotor 4 are assumed to rotate counter-clock-wise CCW, while rotor 2 and rotor 3 rotate clock-wise CW as shown in Fig. 2. Let the forward direction be alongX B .This design can serve for many critical applications.The fact that the motions are completely decoupled and that the quadrotor doesn't need to pitch to go forward nor to roll for lateral motions; this fact makes our design very suitable for sensitive payload as it provides a very smooth ride.Surveillance and monitoring could be improved as the quadrotor can fly at precise attitudes with precise speeds and orientations which make it very suitable for military application.
Let R P ̅ i P i be the rotational matrix from the rotorsrotating frame O P ̅ i to the rotors-fixed frameO P i .Since the rotors' fixed framesO P i are parallel to the bodyfixed frame O B at the center of gravity, then where c(. ) and s(.)denote cos(. ) and sin(. ) respectively.A full rotation of the quadrotor body with respect to the inertial frame can be described by three consecutive rotations about the three body axes, roll rotation Φ about the body x-axis, pitch rotation Θ about the body y-axis and yaw rotation Ψabout the body z-axis.Then R B E is the body transformation matrix with respect to the earth inertial frame, and is given by The relationship between the body-fixed angular velocity vector [ p q r ] T and Euler-Angles rates [ Φ ̇Θ̇Ψ̇] T is given by

Quadrotor dynamics
Let ωi be the rotational speed of the rotor i.Then the lifting thrust is given as bω 2 i and the drag moment is given as dω 2 i, where b and d are the thrust and drag moment constants respectively.Therefore, the thrust components (force) of the i th rotor at the body C.G. are then given by The moments of a titled rotor consist of the drag moment and the moments generated by the thrust components.These two components are expressed as where δ(i) = [1,1,−1,−1] and ri is the vector from center of gravity to the reference point of the rotors, i.e. r 1 = [l, 0, −h], r 2 = [0, l, −h], r 3 = [−l, 0, −h], r 4 = [0, −l, −h] where h and l are the vertical and horizontal displacements from the center of gravity to the rotors respectively.Denote by η the quadrotor position vector and by Ω the body angular velocities vector which are given byη = [x y z] T and Ω = [p q r] T .Then the summation of forces acting on the quadrotor body are given by the dynamic equation: where m is the mass of the quadrotor, gz = [0 0 g] T and K is the matrix of drag constants, and is given by T is the disturbances components of wind gust that will be discussed later.The rotation dynamic equation is then given by: where I is the inertia matrix of the quadrotor, and is given by I = diag{ I x I y I z }.
MG is the gyroscopic forces, and is given by IR is the rotor moment of inertia and M d represents a random disturbance moment, M f is the drag/friction moments with K 4 , K 5 and K 6 representing the drag coefficients (Habib et al., 2011) and is given by: Finally, the equations of motion are where Fi and Mi are given by Eqs. 4 and 5 respectively.

Fig. 2:
Quadrotor with each rotor tilting about two axes

Wind gust modeling
The modeling approach used in this work is based on Solovyev et al. (2015).This approach takes the following into consideration: The wind force expression, depending on the effective influence area on the quadrotor, is also derived.However, this model is suitable, based on the finding in Habib et al. (2011), where the effect of wind gust in small quadrotors is significantly correlated to the rate of increase or duration of a gust rather than the magnitude of the gust.
At any point in time the effect of wind felt at the different elements of the body is assumed to have equal magnitude and direction.The wind model velocity V takes the form: where, t m ; represents the maximum flight time.n; represents a discrete random variable to determine the number of wind steps for t m .
V oi ; represents the wind velocity before each step.
t oi ; represents a discrete random variable to determine each wind step start.d ni ; represents a discrete random variable to determine each duration of gust.
V m ; represents each gust magnitude.
Simulation example of the model given in Li and Wang (2013) for t 0 = [0; 9; 16; 19]s, V m = [1; 4.5; 0; 1]m/s, d n = [7; 5; 2; 5]s, t m = 25s and V 0 = 0.5m/s is shown in Fig. 3. However the following limitations apply when generating random values: < a (a: restriction of the rate of step rise) v i ; represents a discrete random variable to determine each gust magnitude.Also, the point at which wind blows as a wind direction be the azimuth (Ψ w ), measured from the north through east.Wind direction changes at each wind velocity step is given by: where ∆Ψ w is the random value of wind direction change.Since the wind velocity changes with altitude, the average wind velocity is determined by: where, V cz ; wind velocity at the altitude of z.V 0z ; specified wind velocity at the altitude of z o .p; energetic wind profile index.The wind force is given by: where, S e is the effective area influenced by the wind and A is a conversion factor to Nm 2 .
With reference to the influence force, we decompose into the following components for more appropriate or easier application: 2 cos(Ψ w ); F w y = S e AV cz 2 sin(Ψ w ) For simplicity, the quadrotor surface area is represented as a cylinder.So the surface area: The right hand representing the sum of lateral area and bases and μ, σ representing the fill factors here.Therefore, if wind affects only half of the quadrotor the effective area will be given by: S e x = μπrhcos(θ) + σπr 2 sin(θ) S e y = μπrhcos(φ) + σπr 2 sin(φ) with θ and φ representing the pitch and yaw angles.

Problem formulation
A quick preview of the modeling of this tiltable UAV reveals that the system presented has more inputs than degrees of freedom; that is to say, it is overactuated.Taking that into account, it might be desirable to find six control inputs -forces and moments-that reflect a direct effect on the six degrees of freedom.The choice that's most related to the degrees of freedom is the set of forces and moments along the three axes.Let Fv = [Fx Fy Fz] T and Mv = [Mx My Mz] T be the control inputs which are related to the actual 12 inputs through the following static equations: Note that the 6 control inputs [Fx Fy Fz Mx My Mz] are composed of nonlinear combination of 12 actual inputs [ωi αi βi],i = 1,...,4.This means that for each set of control inputs, the system has many combinations of actual inputs.The way to find the best combination is discussed in the next section.
For convenience, the equations of forces and moments ( 6) and (7) will be treated as two subsystems in terms of the new control inputs for the rest of our analysis where Eq. 13 represents subsystem Σ1 and Eq. 14 represents subsystem Σ2 where To construct the formulation of feedback linearization for tracking problem for Σ1, let eη be the error in the position e η = η d − η where η is the position vector, ηd is the desired position vector.Dividing by m and subtracting ηd from both sides of Eq. 13 The control input for this subsystem is F. To introduce feedback linearization for this part of the system through the control input, denoted by Fv, let A similar procedure is carried out for Σ2.Let where the constants KvF, KpF, KvM and KpM are to be determined.The closed loop configuration is shown in Fig. 4.

Feedback linearization and stability analysis
In this subsection, the stability of the system is studied and stated by the following two theorems.
Theorem 1: Subsystems Σ1 and Σ2 are asymptotically stable under control inputs ( 16) and ( 17 A necessary and sufficient condition for the error eη to be locally stable is that V ̇(e η ) ≤ 0. To check this condition, the linearization controller is introduced through the input Fv as follows Substituting Fv into F in Eq. 19 yields Comparing the expression of Fv in Eq. 20 with Fv in Eq. 16 yields KpF = I while KvF is required to be a positive definite matrix.
Similarly, to check the stability condition, the feedback linearization controller is introduced through the input Substituting Mv into M in Eq. 23 Similarly again, comparing the expression of Mv in Eq. 24 with Mv in Eq. 17 yields KpM = I while KvM is required to be a positive definite matrix.
Theorem 2: Assume the drag constant K is not known exactly and is estimated with K ̂, subsystem Σ1 is asymptotically stable if the following conditions are satisfied Introducing the feedback linearization controller through the input Substituting in V ̇(e η , K ̃) and after some manipulations ˙V̇( e η , K ̃) = −K υF (ėη T ėη) − ėη T (K ̃/m)η̇+ s T K ̃TK ̂̇s (30) Note that K ̃ and K ̂̇˙ are 3x3 diagonal matrices while e˙ and η˙ are 3x1 matrices.
To analyze the second and the third terms of Eq. 30 further, the matrices are broken down to their basic elements

Optimization
We mean by Optimization here is the process of selecting the best (optimum) element -or set of elements-from a set of available alternatives according to some preset criteria.In the situation presented here, it is a question of finding the bestaccording to some cost function-combination of inputs that satisfy the control equations.
As mentioned in the previous section, the forces F and moments M are used as inputs to the system to avoid dealing with the nonlinearity present in the actual inputs ωi, αi and βi Feedback linearization control is also achieved through F and M which means that the controller chooses values of F and M to be fed as inputs to the system.It has been also shown in Elfeky et al. (2013) that with only two opposite rotors running, there's always a combination of ωi, αi and βi that satisfy any arbitrary values of F and M.This means that with full actuation, there are many combinations of ωi, αi and β i that satisfy any arbitrary values of F and M. The question addressed in this section is how to choose one combination of ωi, αi and βi from the set of many combinations.
Several optimization techniques are available in the literature.However, few of them deal with nonlinear problem with constraints.Two optimization techniques are tested and compared as a proof of concept and to complete the controller design discussed in the previous section.Those techniques are Genetic algorithm and Interior-point algorithm.This part specifically has high potential for improvement in both theory and implementation.It must be clear that this is not a typical optimal control problem where the objective is to find optimal gains for the controller, instead, the controller gains are already determined and the purpose of optimization is to map the forces and moments used in the controller to the actual 12 inputs of the system through the static equations relating them.
Interior point (IP) is an optimization class of algorithms that can solve linear and nonlinear convex optimization problems.The algorithm reaches an optimal solution by traversing the interior of the feasible region rather than the surface of the region.Convex optimization and particularly interior point are covered in Nesterov and Nemirovskii (1994) and Renegar (2001).
Genetic Algorithm (GA) belongs to a class of algorithms called Evolutionary Algorithms (EA).EA generates solutions using nature-inspired techniques such as inheritance, mutation, selection, and crossover.GA is covered in details in Goldberg (1989), Homaifar et al. (1994), Fonseca and Fleming (1995), and Michalewicz et al. (1992).

The optimization problem formulation
It's first desired to find a cost function for the optimization problem, which means finding some criteria to choose among the many sets of inputs.The most typical and critical objective in UAV application is the minimization of energy consumption during the flight.The cost function can be written as follow: where the first term is the power consumption by the rotors.The cubic exponent comes from the assumption that the torque is proportional to the square of the angular velocity, while the power is the torque time the angular velocity.The second and the third terms penalize tilting movements (cost of energy) and n is the number of samples over which the optimization is performed.
The optimization problem can then be written as follows Two algorithms are tested to perform the optimization problem, Genetic Algorithm and interior-point algorithm.These are implemented using MATLAB functions ga and fmincon.Two algorithms were developed for each method, an offline and recursive algorithms.The objective of the recursive algorithm is to minimize the time for optimization to converge by optimizing for one sample with 12 inputs only at a time instead of optimizing for the 12 inputs in the whole flight samples.In Genetic Algorithm, the recursive technique is developed by setting the final population in one time sample to be the initial population in the next time sample.While in interior point optimization, the recursive algorithm is developed by setting the final solution at one time sample to be the initial set point in the next sample.Those techniques work only if the sampling time is small enough so that the there's no big difference between two consecutive samples.

Simulation results
In this section, the study results will be reported.First the simulation, while considering no wind gust and disturbance is presented.After that the results with wind gust is stated.
First, an initial comparison is carried out to choose which optimization algorithm would be more suitable.The criteria for comparison are convergence time and optimality in terms of cost function.However, the first priority here is time for practical implementation purposes.The four algorithms are tested for 1, 2 and 3 samples of a flight simulation.The results are listed in Table 1 which includes the convergence time and the optimal value of the objective function.It's evident from the results that Recursive Interior-Point (RIP) technique is more practical for our purpose.Although conventional IP technique is superior in optimality, which is expected, RIP takes Significantly, less time to converge which makes it the best candidate for online implementation?Further study is carried out on the developed RIP to investigate its behavior in terms of convergence time.Fig. 5 shows the convergence time against number of samples optimized.The curve is almost linear except near the first few samples as shown in Fig. 6.The slope of the linear part is equal to 0.270 m/sample.This means that the algorithm takes around 270 ms to optimize one sample of 12 inputs.This result is not suitable for online implementation, therefore all flight simulations are implemented offline.However, this result is very promising for future research.
To test the control algorithm, the quadrotor is commanded to elevate up to 10m, tilt forward with 30 degrees and perform a circle while pointing to its center.To perform a circle, the quadrotor is commanded to follow sinusoidal paths along x and y axes.The flight is an example of a surveillance mission where the quadrotor might be taking panoramic photos of a certain target.This complex flight demonstrates the strength of both the design and control technique.The optimization is done offline, i.e., and prior to the flight.The parameters of the system are assumed to be known.The control parameters are chosen, according to Theorem 1. KpF = KpM = I while KvF = KvM are taken to be equal to 3I, where I here is the 3x3 identity matrix.Fig. 7 shows both the desired and actual positions of the quadrotor in x, y and z axes while Fig. 8 shows the first derivative of the position vector and the 3-D position is showed in Fig. 9. It's evident from the figures that the path was followed with high accuracy.Finally, the simulation is done with wind gust and disturbances.Fig. 15 shows the elevation position of the command to take off vertically and then thrust further in x direction to a height of [10, 0, 50] m.In Fig. 16, the 3D plots are shown.Here, the controller is able to track the coordinate under wind disturbance with a negligible error due to wind.Fig. 17 which depicts the coordinates of the quadrotor moving in the x direction with a slight error in the y position as shown due to wind gust disturbance.Fig. 18 shows some fluctuation in the velocities due to effect of wind disturbance.Fig. 19 and Fig. 20 shows the angles and the angular velocities are unaffected under wind disturbance.This shows the capability of the quadrotor itself to move towards a certain direction without compromising the attitude under wind disturbance.

Conclusion
Feedback linearization proved to stabilize and control the attitude of the quadrotor with tilting rotors.All the degrees of freedom can be controlled freely and independently.This allows for more complex and precise maneuvers that weren't possible with conventional quadrotor.In summary, this control algorithm highlights the advantages of the novel design by taking advantage of the high mechanical flexibility available.In fact, the controller decouples the model of the quadrotor into two completely separate systems, systems of forces and positions, and a system of moments and orientations.The two systems can be studied and controlled absolutely independently.The only connection between the two systems appears in the final stage of optimization where the six equations are solved together to find the 12 inputs.The optimization part is developed to choose the best inputs according to energy consumption which is very critical in all UAV applications.A nested PD control loop is adopted to control the quadrotor position and cancel disturbances with additional inputs due to the quadrotor structure.The wind model is derived from an approach which takes into account the wind velocity change, the wind gust step change, the variation of wind velocity with increase in height and the changes in wind direction.
Finally, we can say that the control and optimization part show very promising results.The results demonstrate that the quadrotor with tilted rotor can effectively attain the desired trajectory in the presence of wind disturbance.However, there is still a substantial space for improvements.Future research may include a wide variety of angles and considering different scenarios of uncertainty in the parameters, states and disturbance.The optimization part as well can be further improved to consider different techniques and online implementation of optimization.
Other control techniques can also be developed for this system.Particularly, model predictive control MPC is a candidate due to the fact that it combines control and optimization under constraints.


The effect of wind velocity change (increasing or decreasing). Gust duration. Wind velocity change with respect to altitude. Wind direction change.

Fig. 3 :
Fig. 3: Simulation showing wind velocity before and after wind gust

Fig. 7 :
Fig. 7: x, y and z positions of the quadrotor The orientation vector of the quadrotor is shown in Fig. 10 in radians.It can be seen that the angle reached π/6 and the angle went from 0 to 2π to keep pointing towards the center of the circle while the quadrotor is completing a full circle.This means that Yaw angle is also performing a complete circle simultaneously.In Figs.11 and 12, the input angles α i and β i generated by optimization are shown.One objective of optimization was to limit the change in

Table 1 :
Convergence time and objective function comparison