Effect of Rotor Tilt on the Gust Rejection Properties of Multirotor Aircraft

: In order to operate safely in windy and gusty conditions, multirotor VTOL aircraft require gust resilience. This paper shows that their gust rejection properties can be improved by applying a small amount of ﬁxed outward rotor tilt. Standard aerodynamic models of the rotors are incorporated into two dynamic models to assess the gust rejection properties. The ﬁrst case is a conceptual birotor planar VTOL aircraft. The dependence of the trim and stability on the tilt angle are analyzed. The aircraft is stabilized using a pole-placement approach in order to obtain consistent closed-loop station-keeping performance in still air. The effect of gusts on the resulting response is determined by simulation. The second case study is for a quadrotor with a 10 ◦ outward rotor tilt. The aerodynamic coefﬁcients are analyzed for trimmed station-keeping over a range of steady wind speeds. An LQR controller is used to apply station-keeping that includes integral action, and the gust responses are again obtained using simulation. The results show that the outward rotor tilt causes the aircraft to pitch down into a lateral gust, providing lateral force that opposes the gust and so signiﬁcantly improving the gust rejection properties.


Introduction
Multirotor Vertical Take-Off and Landing (VTOL) aircraft have become very popular for many applications such as surveillance, monitoring, search, aerial surveying and photography, structure inspection and so on. They are also being proposed for personal transport aircraft (air-taxis) and low-cost delivery services. However, many of the most useful and high-value applications occur in environments with gusty wind conditions, such as off-shore structures, rugged and mountainous terrain and urban areas. Because most multirotor VTOL aircraft are very light, they are very sensitive to wind gusts. This limits their operation in gusty conditions, particularly for applications where accurate hovering and/or landing is required.
Most literature on quadrotor control does not explicitly consider the problem of gust rejection. Sometimes the gust rejection properties of the controller are evaluated in simulation, for example [1][2][3][4]. Some explicitly consider the gust rejection properties in the controller design-but by a reactive scheme utilizing feedback and often a gust observer (for example [5][6][7][8]), or even on-board gust measurements [9]. The performance of such systems will always be limited by the actuator magnitude, slew rate and bandwidth constraints.
Some work has been done to analyze the aerodynamics of quadrotor vehicles or to include the aerodynamics in control system models. Analytical modelling of the rotor aerodynamics is usually based on the established methods for rotorcraft [10,11]. The effects of horizontal wind and gusting have been analyzed with this approach [12]. Another noted work is by Bristeau et al. [13], where the effect of the centre of gravity location on the stability and the lateral gust rejection were analyzed. Other quadrotor work that incorporates aerodynamic modelling includes [14][15][16], while Del Cont Bernard et al. [17] have studied ground effect. In particular, Lopez et al. [18] have explored the gust rejection properties of some multirotor configurations using flight test.
In this paper, the effect of using a fixed tilt on the rotors as a means of reducing the gust sensitivity is analyzed. The idea is that by a small fixed-angle outward tilt of the rotors, the aerodynamics can be tailored so that a horizontal gust on the vehicle causes the aircraft to pitch or roll down into the wind and so reducing the horizontal position displacement. The approach was proposed by Whidborne and Cooke [19] where the effect on a birotor Planar VTOL (PVTOL) aircraft was analyzed and shown to have some potential in reducing the gust sensitivity. A PVTOL aircraft is a simplified birotor aircraft model that is constrained to move only in the vertical plane. It is commonly used as a model for testing multirotor control concepts, e.g., [20]. Note that the idea of reducing the gust sensitivity by rotor tilt was also proposed by Otsuka and Nagatani [21] who conducted some experimental studies to measure the pitch/roll moment effect [22].
Standard techniques are used to include aerodynamic effects in the equations of motion so that the effect of rotor tilt on the gust response of multirotor vehicles can be analyzed. Two multirotor configurations are considered, a birotor PVTOL aircraft and a standard + configuration quadrotor. The proposed analysis methods can easily be extended to other multirotor configurations. The effect of the rotor tilt on the stability of the PVTOL aircraft is analyzed using small perturbation methods. Feedback control is introduced to stabilize the aircraft, and the effect of gusts on the lateral position is investigated by simulation. For the quadrotor model, the influence of the rotor tilt on the aerodynamic derivatives is investigated. Feedback control to maintain the position of the vehicle is applied and the influence of gusts is investigated for the quadrotor with and without outward rotor tilt. Some parts of the paper appeared in a similar form in [19], but are included here for completeness and comparison purposes. The paper also corrects an error in Equation (8) of [19] to present more accurate results for the birotor PVTOL aircraft. The analysis is extended beyond that of [19] to a quadrotor aircraft. The aerodynamic model used in the simulations is that of a Draganfly X-PRO quadrotor that has been validated in a wind tunnel [23]. Furthermore, the MATLAB codes for implementing the models and simulations are available for the community in [24], where an animation of the effect of the rotor tilt on the PVTOL is also included.
In the next section the aerodynamic models used to determine the forces and moments generated by a rotor are presented, and results for a particular rotor example are shown. In Section 3, the dynamic model of the PVTOL aircraft is presented and trim conditions using the rotor example are given. In Section 4, the linearized dynamics are analyzed and in Section 5 the PVTOL controller and gust responses are presented. The quadrotor controller, aerodynamic derivatives analysis and gust responses are given in Section 6, and finally some conclusions are given in Section 7.

Rotor Aerodynamic Model
In order to develop a model of a multirotor gusting response, the aerodynamics of a single rotor are first considered. Figure 1 shows a rotor disc in a horizontal airflow, V. The axial direction rotor thrust, T, is given by [25] where ρ is the air density, A is area of the rotor disc, Ω is the angular velocity of the rotor, R is the rotor radius and C T is the rotor thrust coefficient. The in-plane horizontal force, H, (i.e., the force in the plane formed by T and V) is given by where C H is the rotor horizontal force coefficient.
Note that in this paper, we assume that the in-plane moment caused by the air flow is insignificant. This is generally true for hinged blades, however most multirotors use unhinged blades.
Finally, the rotor torque, Q, about the axial direction is given by [25] where C Q is the rotor torque coefficient.

Thrust Coefficient
The thrust coefficient can be determined from momentum theory which was originally developed by Glauert [27] building on the earlier work of Froude on aircraft propellers. Ignoring the flapping term, the thrust coefficient is derived as ( [26], p. 59) (also see [28], p. 98): where σ is the disc solidity factor, a is the rotor blade lift slope, θ is the blade pitch angle, µ is the rotor advance ratio (also the lateral wind component), and λ is the rotor inflow ratio. For a rotor disc with b blades of constant chord, c, and radius R, the solidity factor is defined as σ = bc/(πR) ( [28], p. 48). For tapering blades, the chord can defined as the chord at the 3/4 radius ( [28], pp. 49, 98). For blades with a linear twist, the equation is still valid if the pitch angle at 3/4 radius is used ( [28], p. 98), that is θ = θ 0 + 3/4 θ 1 where θ 0 is the blade pitch at the roots and θ 1 is the twist. The rotor inflow ratio is the ratio of the axial flow velocity of the air to the tip speed of the rotor blade. Here it is considered positive in the direction of induced velocity, which is opposite to the thrust force and is given by ( [26], p. 56) where V is the incident flow velocity, v i is the induced velocity and α is the angle between the incident flow and the tip-path plane as shown in Figure 1.
The rotor advance ratio is the ratio of the relative air velocity in the plane of the rotor and co-linear with the horizontal force to the tip speed of the rotor blade. It is considered positive in the direction of the horizontal force unlike the convention for the in-flow ratio and, ignoring the horizontal component of the induced velocity, is given by ( [26], p. 56) In this study, it is assumed that µ takes a maximum value of 0.5 ( [28], p. 96).
If it is assumed that the induced velocity is constant over the whole disc, the induced velocity can be calculated using the expression ( [26], p. 52) where v h is the induced velocity with the rotor in hover at the same thrust and is given by

In-Plane Force Coefficient
Along with vertical thrust, the rotor disk also has an in-plane horizontal force when the incident flow has a lateral component. This results mainly from the blade profile drag. The horizontal force coefficient is given by (note that this equation is incorrectly given in [19]) ( [29], p. 149) where c d0 is the blade profile drag coefficient and a 1 is the backward flapping angle which can be calculated from ( [26], p. 65, Equation (5.50))

Blade Torque Coefficient
If the blade profile drag is represented as the polar where α a is the angle of attack of the blade, then the torque coefficient can be determined by a detailed expression [30] that if c d2 is insignificant can be simplified to

Force and Moment Mappings
Multirotor vehicles are controlled by varying the rotor speed, Ω of each rotor. To design the controller and analyze the dynamics, the mappings from Ω to the triple (T, H, Q) are required. However, the triple also depend upon the airspeed magnitude, V, and incidence angle, α. Hence the mappings from the triple (Ω, α, V), onto (T, H, Q) must be calculated.
The mappings (Ω, α, V) → (T, H, Q) are described by Equations (1)- (12). However, the presence of T in (8) means the equations are implicit and it appears that there is no tractable explicit solution. Calculation of the mapping thus requires a numerical method. In this paper, the MATLAB [31] routine, fzero, is used. The algorithm is based on the Brent-Dekker algorithm [32] and uses a combination of bisection, secant, and inverse quadratic interpolation methods [33,34] and includes an extension to check for sign changes. Other nonlinear function root-finding methods could be used.
In this paper, we provide an example that uses the rotors for the Draganflyer X-Pro quadrotor shown in Figure 2. The aircraft has undergone extensive wind tunnel and other experimental testing [23] from which we obtain the relevant parameters, a = 5.5, b = 2, c = 0.04 m, R = 0.258 m, δ = 0.05, and θ 0 = 0.3025 rad. The mapping for (Ω, α) → T for constant airspeed values V ∈ {0, 10, 20} ms −1 is shown in Figure 3. Similarly, the mapping for (Ω, α) → H for constant airspeed values V ∈ {0, 10, 20} ms −1 is shown in Figure 4, and the mapping for (Ω, α) → Q for constant airspeed values V ∈ {0, 10, 20} ms −1 is shown in Figure 5. The variation of thrust with airspeed is shown by the mapping (α, V) → T for constant rotor speed Ω = 150 rad s −1 and is shown in Figure 6.

PVTOL Equations of Motion
In order to understand how rotor tilt affects the gust rejection properties of multirotors, the dynamics of a 2-rotor PVTOL vehicle are modelled, the rotor aerodynamics models from Section 2 are included. The aircraft is a conceptual vehicle that is constrained to operate in the vertical (x, z)-plane as shown in Figure 7. The rotors are symmetrically tilted outwards with an angle Υ. The weight of the vehicle is W, x and z are the horizontal and vertical positions of the centre of gravity of the vehicle, θ is the pitch angle, u and w are the horizontal and vertical velocity components in the body axes, q is the pitch rate, V W is the wind which is assumed to be horizontal,  The equations of motion in the body axes are where s (·) denotes sin(·), c (·) denotes cos(·), I is the moment of inertia, m is the vehicle mass, is the arm-length, and W = mg where g is the gravitational acceleration constant. We also have the navigation equations that relate the rates to the aircraft position in the airspace:

Vehicle Trim
The aircraft steady state (trim) conditions are determined next. This is a prerequisite for performing linear analysis of the aircraft dynamics. The aircraft is trimmed at hover in the earth coordinates with a steady horizontal wind. The trimmed variables are the pitch angle and the left and right rotor speeds denoted byθ,Ω andΩ r , respectively. Setting the position and pitch rates and accelerations to zero in (13)- (22) gives Note that T i and H i are dependent upon the rotor speed, Ω i , but the implicit nature of the thrust mapping means a numerical method is required to solve (23)- (25). The approach used here is to minimize the sum-of-squares residual of (23)-(25) using the MATLAB routine, fminsearch, which applies the Nelder-Mead simplex search method [35,36]. Note that other methods could be used. Additionally, note that for the still air condition (V = 0) the solution is trivial and is given byθ = 0 andT = W/(2 cos Υ).
The PVTOL geometry is based on the Draganflyer X-Pro quadrotor, and is defined as I = 0.0625 kg m 2 , m = 2.36/2 kg, = 0.45 m and Υ = 0 [23]. The resulting trim maps for V W ∈ [0, 20] ms −1 are shown in Figure 8. Note that for Υ = 0, from (25)T l =T r and henceΩ l =Ω r =Ω. In particular, it can be seen that for low horizontal wind speeds (below about 15 ms −1 ), the trim rotor speed is below the hover value. This may seem counter-intuitive, but is caused by the lift that the wind provides combined with the fact that the vehicle drag is ignored, which would cause energy loss in a real vehicle.  Figure 9 shows the trim map for Υ = 10 • . The difference between the mean trim rotor rotation rate, (Ω l +Ω r )/2, and the trim rotor rotation rate for Υ = 0 is just 1 rad s −1 for hover in still air, is about 5 rad s −1 at V W = 5 and is much more significant at high wind speeds, with about a 17 rad s −1 difference at V W = 20 ms −1 . This means that performance degradation resulting from the rotor tilt is increased at higher wind incident speeds. Figure 10 shows the dependence of the mean trim rotor speed on the tilt angle.

Model without Rotor Aerodynamics
The effect of the rotor tilt angle, Υ on a linearized model of the dynamics in still air, that isV W = 0, is analyzed first. In addition, the aerodynamics of the rotor are ignored, so the control is effected purely by the thrusts, T i , with the horizontal components, H i , assumed to be zero. From (13)-(15) and (20)-(22), the state vector is defined as x = [z, w, x, u, θ, q] T and the control vector Assuming small perturbations about hover trim, the linear model is given bẏ For the output y = [z, x, θ] T , the system Laplace transfer function matrix is given by The lateral position channel has a pair of zeros, s = ± mg /(I tan Υ). Provided Υ is small, the zeros are high frequency and hence have little impact on the closed loop system performance. If Υ > 0, then there is a nonminimum phase zero. Note that A T and hence the open-loop stability is independent of Υ. However B T is dependent on Υ, and some direct control of the lateral position rate is achieved for Υ = 0, but this is adverse for Υ > 0, resulting in the nonminimum phase zero.

Model with Rotor Aerodynamics
The linearized model is extended to include the aerodynamic effects. The control is taken as and the aerodynamic forces, T i , H i , are dependent upon the state, x, control, u and wind disturbance V W . The dynamic equations, (13)- (22), are linearized at a set of trim points using the second order finite difference method ( [37], p. 200). This gives the small perturbation model where u w = V W −V W is the external disturbance. The resulting state matrix is where Z w = ∂ẇ/∂w, X u = ∂u/∂u and M q = ∂q/∂q are the concise aerodynamic damping derivatives, X q = ∂u/∂q and M u = ∂q/∂u are concise cross-coupling aerodynamic derivatives. The remaining aerodynamic derivatives are all nearly zero. Figure 11 shows the variation of the derivatives {Z w , X u , X q , M u , M q } with Υ. Note that for Υ = 0, the cross-coupling derivatives X q and M u are both zero and the resulting eigenvalues are M q , X u , Z w , 0, 0, 0 , as expected from (30). The dependence on Υ of the real and imaginary parts of the eigenvalues is shown in Figures 12 and 13, respectively. Note the high eigenvalue sensitivity when Υ is near zero. Using the Routh-Hurwitz criterion, stability conditions can be determined. The system characteristic equation, Q(λ) = det(λI − A), is For Υ > 0, the final term in (31), M u g, is negative hence instability. For Υ < −4.04 • the system is also unstable, which can be shown to occur for The control input matrix B u is where Z u1 = ∂ẇ/∂u 1 , X u2 = ∂u/∂u 2 and M u 2 = ∂q/∂u 2 are concise aerodynamic control derivatives. The remaining derivatives are all nearly zero. The variation of the control derivatives, {M V W , X V W , Z V W }, with Υ is shown in Figure 14. Derivatives Z u1 and M u 2 are fairly constant, but X u2 changes sign to negative for Υ > 0. Taking y = [z, x, θ] T as the output, the system Laplace transfer function matrix is given by Hence there is a zero pair in the horizontal position channel, G 2,2 (s), and a zero in the pitch channel, G 3,2 (s). For Υ < 0, the zeros of G 2,2 (s) are left half plane and complex. The real part is fairly constant, approximately −0.4, and the imaginary parts approach infinity as Υ approaches zero. For Υ > 0, the zeros are real, with one in each half plane. As Υ increases, the zeros move from ±∞ towards each other, and at Υ = 15 • the zero pair is at (−18.0, 17.2). The zero in G 3,2 (s) is located at s = (M u2 X u − M u X u2 )/M u2 and it does not vary much and has a value of about −0.21.
Inspection of the disturbance input matrix, B w , shows the sensitivity of the pitch rotation rate to the wind disturbance, V W , where Figure 15 shows the variation of the disturbance derivatives, {M V W , X V W , Z V W } with Υ. Derivative Z V W is zero, X V W is always positive as expected, and increases with increasing tilt angle magnitude. The moment derivative, M V W is proportional to Υ, being zero at Υ = 0. We see that although Υ > 0 causes instability, it also provides a pitch down moment to the vehicle that will cause it to tilt into the wind when subjected to a lateral gust. This is potentially very beneficial for gust rejection. This property is explored next.

PVTOL Linear Control
In this section, the dependence of the gust rejection properties of a PVTOL aircraft on the tilt angle, Υ, is assessed. A set of four stabilizing controllers are designed. State feedback is assumed. The linearized model for hover trim in still air given by (29), (30) and (33) is used as the nominal model. Good control of multirotors for non-aggressive manoeuvres is easily obtained using standard Linear Quadratic Regulator (LQR) control (for example [38,39]). This was the approach used in [19], but it was difficult to compare the controllers for the different values of Υ, so a different approach is used here. Because the vertical position channel and the horizontal position/pitch channel are very decoupled, the controllers for each channel can be designed independently as two SISO problems. We use pole-placement in order to obtain a consistent closed-loop performance for stationkeeping in still air. The resulting gust-rejection properties can then be fairly compared.
The closed-loop pole locations are set to be close to the closed-loop pole locations for the LQR controller for Υ = 10 • from [19].  Figures 16 and 17, respectively. The responses are nearly identical; the poles have been set so the position keeping is not too aggressive, so with a maximum pitch angle of about 16 • , the system dynamics remain fairly linear. The responses remain close until the step demand exceeds around 25 m. Note that there is a small amount of initial undershoot in the step response of x for Υ > 0 resulting from the non-minimum phase, but it is barely perceptible in Figure 16. The vertical interaction is small, as shown in Figure 18, where small differences can be seen in the responses for the different values of tilt angle, Υ. The rotor rotation rates are shown in Figure 19. It can be seen that the rotors have to work the hardest for Υ = 10 • and the least for Υ = 0 • , which is expected.   Responses of the horizontal position to a lateral wind velocity step of 5 ms −1 are presented in Figures 20 and 21. Video visualizations of the responses are available [24]. The gust rejection properties of the vehicle with Υ > 0 are remarkable. For Υ = 10 • , the vehicle actually pitches into the wind, so much so that the final lateral position is negative, meaning that the vehicle actually moves up-wind. For Υ = 5 • , the vehicle almost exactly maintains horizontal position by pitching into the wind. Note that the controller does not contain integral action, so it does not compensate for the steady state position error. For Υ = −5 • , the vehicle pitches away from the wind, so that the effect of the gust on the lateral position is severe, resulting in a 12 m displacement. Figure 22 shows the resulting vertical position response; for all angles of Υ, the vehicle finally gains an altitude of nearly 1 m. All the responses initially gain altitude because of the additional lift provided by the wind. The Υ = 10 • response then actually loses some altitude because the vehicle has pitched over 20 • and so loses enough vertical component of the rotor thrust to drop below the initial altitude; however this loss is recovered later. Figure 23 shows the rotor speed responses. As expected, the Υ = 10 • case rotors have to work the most to stabilize the vehicle because of the greatest pitch angle. The Υ = −5 • rotors have to work harder than the Υ = 5 • case because the negative pitch has to be recovered.

Equations of Motion
In this section, the analysis is extended to assess the gust resilience of a quadrotor with fixed rotor tilt. The schematic of the quadrotor configuration is shown in Figure 24. The rotors are symmetrically tilted outwards with an angle Υ. The earth reference frame is denoted by E, the body-attached frame by B and the four rotor frames by R 1 to R 4 . To model the vehicle dynamics, it is assumed that the structure is rigid and symmetric, the rotor Coriolis forces are insignificant, the motor dynamics bandwidth is high and can be ignored, and the centre of mass coincides with the origin of the body axes, B. The equations of motion are then the standard Newton-Euler equations where F is the total force vector acting on the center of mass in the body axes, m is the mass of the body, I 3 is the 3 × 3 identity matrix, v b is the velocity vector of the center of mass, M is the total moment vector acting about the center of mass, I b is the inertia matrix about B and ω is the angular velocity vector of the body. The body frame translational rates are related to the earth frame rates bẏ p = R(φ, θ, ψ) v b where p = [x, y, z] T is the position of the centre of gravity in the earth frame, (φ, θ, ψ) is the standard aerospace Euler angle (Tait-Bryan) triple and is the Direction Cosine Matrix (DCM). Finally the Euler angle rates are related to the body angular velocities by  φ θψ The rotor aerodynamics are modelled as explained in Section 2. The thrust, horizontal disc force and torque triple, (T i , H i , Q i ), for the ith rotor depend on air flow speed and incident angle on the rotor tip plane. Denoting the ith rotor axis unit vector by z i as shown in Figure 24 and the corresponding air flow vector by v a i , the incident angle on the ith rotor tip plane, α i , is given by where v w is the wind vector and r i is the ith vehicle arm vector, that is the vector from B to R i as shown in Figure 24. The resulting force and moment vectors are respectfully given by The resulting total force, F is given by where z e is the earth down axis unit vector shown in Figure 24, and the total moment vector, M, is

Quadrotor Steady State Analysis
In the same manner as the birotor in Section 4.2, the model is linearized and obtained in the form of (29). However, to aid the analysis, here we consider a reduced order form considering just the dynamic Equation (36), with a state variable vector x r = [u, v, w, p, q, r] T . This removes the kinematics from the analysis whilst leaving the aerodynamics.
The resulting state matrix for the reduced order model is where Z w = ∂ẇ/∂w, etc. The remaining derivatives are zero or very small and can be ignored. The model with a tilt angle Υ = 10 • is trimmed with an x-direction for a headwind of V w ∈ [0, 15] ms −1 . For incident airspeeds above approximately 16.6 ms −1 , the rotor advance ratio (see (6)) exceeds 0.5, and so the reverse flow region becomes significant and the aerodynamic model is less reliable ( [28], p. 96). The pitch trim variation shows similar trends to Figure 9, so is not shown here. The variations of the aerodynamic derivatives of (42) are of more interest. Figure 25 shows the variation of the aerodynamic derivatives, X u , X w , Y v , Z u , Z w , with V w . The lateral y-direction force term, Y v , is almost constant and negative, indicating stability and independence from the incident airspeed. The x-direction force term, X u , is equal to the Y v for V w = 0 with an increasing magnitude with V w , resulting from the increasing incident angle required to maintain trim. The variation in Z w results mostly from the variation in rotor speed, with a lower magnitude with lower rotor speed. There is no cross-coupling for V w = 0. The z-direction force initially increases providing lift in moderate winds, but levels-off as the pitch increases. The z-direction force derivative increases linearly with wind speed, and pitch angle.  Figure 26 shows the variation of the aerodynamic force derivatives with respect to the rotation rates, X q , Y p , Y r , Z q , with V w . The terms are small compared to {X u , X w , Y v , Z u , Z w } shown in Figure 25. The terms all arise from the rotor tilt causing forces that are not parallel or orthogonal to the body axes.  Figure 27 shows the variation of the aerodynamic moment derivatives, L p , M q , L r , N p , N r with V w . The pitch and roll rate damping derivatives, L p and M q , are equal and fairly constant with additional damping at high speed after an initial small reduction. The yaw rate damping derivative, N r , is, as expected, smaller than the roll and pitch derivatives and fairly constant. The cross-coupling terms are zero in still air. The roll moment derivative with respect to yaw rate shows the same trend as for Z u . The magnitude of the N p term increases approximately linearly with V w , these terms arise from the rotor tilt.  Figure 28 shows the variation of the aerodynamic moment derivatives with respect to the body velocities, L v , M u , M w , N V , with V w . These terms arise from the rotor tilt. The pitch moment derivative with respect to z-direction velocity increases in magnitude with the wind speed. The yaw moment derivative is small.
The plots show that rotor tilt introduces some dynamic cross-coupling. The effect on the closed-loop system is small as is seen in the next subsection.

Quadrotor Control
In order to investigate the gust rejection properties of the tilted rotors, LQR controllers are synthesized. The LQR control problem is well known. Givenẋ(t) = Ax(t) + Bu(t), a control input u(t) = −K c x(t) is determined such that closed loop systemẋ = [A − BK c ]x(t) is stable with a gain K c that minimizes

)Qx(t) + u(t)Ru(t))dt
where Q and R are constant weighting matrices.
The objective is to control the positions, x, y, z, and yaw, ψ, to maintain hover. Thus, the system output vector is set to y = [x, y, z, ψ] T and to ensure zero steady state errors, integral action in the controller is implemented by augmenting the linear model with the integral of the error, e i = e dt, whereė i = e = r − y where r is the reference input vector. Hence, given y = Cx, then the augmented system is The linearized control schematic is shown in Figure 29. Note that for a practical LQR controller, an additional nonlinear element should be included in the LQR architecture to rotate the x and y errors by ψ to align them with the body axes in the horizontal plane. That is not required in this study because ψ is small for the simulations. Figure 29. Draganflyer X-Pro quadrotor LQR control schematic with integral action.
Linearized models for hover trim in still air are used to obtain the model, and the weighting matrices are set to Q = I 16 and R = I 4 /2. Nonlinear system simulations are performed to evaluate the responses to a lateral position step demand for two values of tilt, Υ = 0 • and Υ = 10 • . Figure 30 shows the position responses to a 10 m step demand in x for the nonlinear model with 0 • tilt and 10 • outwards tilt. The responses are very similar, thus providing a fair means for evaluating the effect of the tilt. The resulting attitude responses are shown in Figure 31. Here the 0 • tilt pitch response is noticeably greater than the 10 • tilt response. This is because when the aircraft pitches to move forward, the front rotor (R 1 ) has an additional 10 • tilt, and this rotor is exploited by the controller to give additional force in the x-direction but requiring less pitch. Figures 32 and 33 show, respectively, the position and attitude responses to a 10 ms −1 gust step in the x-direction. Both position responses show an initial lateral displacement downwind. For the untilted case (Υ = 0 • ), the controller responds by pitching the aircraft into the wind until the x-position is recovered and settles with a trimmed pitch of 13.3 • . The maximum displacement is about 1.4 m and the maximum pitch is 21.4 • occurring at 1.35 s. For the tilted case (Υ = 10 • ), the gust causes the aircraft to pitch down into the wind which, after a short downwind displacement, causes an upwind displacement. The controller responds by reducing the pitch so that the x-position is recovered, settling with a trimmed pitch of just under 14 • . The maximum downwind displacement is about 0.3 m and the maximum upwind displacement is just under 0.8 m. The maximum pitch is about 6 • greater than the untilted case, and occurs far earlier after just 0.73 s. This demonstrates the ability of the tilted configuration to withstand wind gusts.   Note that both configurations result in an initial altitude gain due to the additional lift from the lateral wind as seen in Figure 8. With close observation it can be seen that there is a small amount of cross-coupling to the y-direction and the yaw, ψ. It is slightly greater for the tilted configuration resulting from the rotor tilts meaning the force and moments vectors are not aligned with the body axes. However, with a tilt of just 10 • , the amount is small.

Discussion and Conclusions
The problem of gust resilience for multirotor VTOL aircraft is considered. It is proposed that introducing an outward tilt to the rotors improves the gust rejection properties of the vehicle. The effect is analyzed on both a PVTOL aircraft and a quadrotor. The following conclusions can be drawn: • The gust rejection properties of multirotor VTOL aircraft can be analyzed and simulated by modelling the aerodynamics of the rotors using standard methods. The methods are suitable for simulation studies that include the aircraft dynamics, but include a nonlinear implicit equation which must be solved numerically. This can be done rapidly using standard methods. • Trim analysis for the PVTOL shows that for operation in moderate wind (or for cruise) additional lift is generated by the rotors. • The effect of rotor tilt on the dynamics of the PVTOL is explored and analyzed. It is seen that an inwards tilt results in increased open-loop static stability, although dynamic stability can be lost for large tilt angles. However outward tilt results in a loss of static stability and some non-minimum phase. For small tilt angles, the nonminimum phase behaviour is high frequency and hence has little effect on the control performance, but the loss of static stability means that some 'hunting' may result at hover. • Simulations of both the PVTOL and quadrotor show that the outwards tilt causes a pitch down moment when subjected to a lateral wind increase. This anticipates the necessary control action to recover the lateral position by aiming the vehicle thrust into the wind and so improves the gust rejection property beyond that provided by the feedback control.
It should be noted that, in essence, outwards tilt introduces anhedral into the aircraft. The gust disturbance rejection properties of anhedral in fixed-wing aircraft are fairly well-known, and is the reason why VTOL fixed-wing aircraft such as the Harrier have notable anhedral. The dynamics of such aircraft have been explored by Hauser, Sastry and Meyer [40]. It should also be noted that most commercial surveillance and observation quadrotor platforms have the camera located below the centre of gravity, which gives the vehicle some inherent gust resilience because the drag of the camera causes a pitch down into a wind gust. Detailed investigation into this could be the subject of further study.
Apart from loss of static stability, the introduction of fixed rotor tilt has another disadvantage. It is obvious that rotors with non-parallel rotor axes will be less efficient and reduce the endurance of the vehicles in still air. This is not investigated in detail in this paper, and the trade-off between reduced efficiency and improved gust resilience is a matter for further study.
The pitching into the wind that results from outwards rotor tilt could also be achieved by aerodynamic 'tailoring' of the vehicle body shape, a detailed study of aerodynamic design trade-off remains for future work. Finally, validation of the effect of outward tilt requires experimental validation, both by wind tunnel test to add to that of Otsuka, Sasaki and Nagatani [22] and by flight test.

Data Availability Statement:
The MATLAB software that generates the data and a video visualization are available in the Cranfield Online Research Data (CORD) repository, [24].

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