Modeling and Closed Loop Flight Testing of a Fixed Wing Micro Air Vehicle

This paper presents the nonlinear six degrees of freedom dynamic modeling of a fixed wing micro air vehicle. The static derivatives of the micro air vehicle are obtained through the wind tunnel testing. The propeller effects on the lift, drag, pitching moment and side force are quantified through wind tunnel testing. The dynamic derivatives are obtained through empirical relations available in the literature. The trim conditions are computed for a straight and constant altitude flight condition. The linearized longitudinal and lateral state space models are obtained about trim conditions. The variations in short period mode, phugoid mode, Dutch roll mode, roll subsidence mode and spiral mode with respect to different trim operating conditions is presented. A stabilizing static output feedback controller is designed using the obtained model. Successful closed loop flight trials are conducted with the static output feedback controller.


Introduction
The micro air vehicles (MAV) are small aircraft with a wingspan of 150 mm and a maximum takeoff weight less than 200 g [1,2]. MAV can be deployed in confined air spaces where it is arduous for bigger UAV (unmanned air vehicle) to operate. The three categories of MAVs are the fixed wing, flapping wing and rotary wing [3][4][5]. Among the three types of MAVs, fixed wing MAVs have higher flight velocity and longer endurance when compared to flapping wing and rotary wing counterparts. Fixed wing MAVs face many unique challenges that make their design and development difficult. The development of a fully autonomous fixed wing MAV is a very challenging task. The fixed wing MAVs operate in low Reynolds number (Re) regime where many complex flow phenomena take place within the boundary layer. The Separation, transition, and reattachment of air flow can all occur within a short distance along the chord line of the wing [1].
Since MAVs are lightweight, low moment of inertia air vehicles, they are very much susceptible to atmospheric disturbances [6]. A robust feedback control is required for MAVs to operate in such adverse atmospheric conditions. A feedback controller is designed based on the mathematical model describing the dynamics of the MAV. To construct the mathematical model of the MAV, the forces and moments acting on MAV are estimated through the wind tunnel testing and computational fluid dynamic (CFD) analysis. MAVs have low aspect ratio wing (AR < 3). The forces generated by propeller wake is significant for MAVs [7]. Hence, a conventional aircraft modeling techniques that do not take into account of propeller wake are inadequate to model the aerodynamic forces acting on MAV. The aerodynamic effects of propeller flow on MAV and small UAV are studied in [7][8][9][10][11][12][13][14]. The influence of propeller wake on the lift, drag and wing stall characteristics is reported in [7][8][9][10][11][12][13][14]. There is an increase in lift force generated by rotating propeller when compared to stationary propeller conditions. The drag force with a rotating propeller is considerably higher than that in the case of the wing without the rotating propeller. The stall of the wing is delayed by 10 • of an angle of attack (α) with the aerodynamic effects of rotating propeller. In [15], a six degree of freedom nonlinear model for MAV is developed without including the propeller effects. Thereafter, a linear model which includes terms which are often neglected in the linear model of a bigger aircraft is constructed by linearizing the six degrees of freedom nonlinear model around a trim point. A specific nonlinear model of a gun launched MAV (GLMAV) for hover and near-hover flight conditions are presented in [16]. In [17], the aerodynamic modeling of the longitudinal aerodynamics of MAV at the high angle of attack under unsteady conditions is carried out. For a 75 mm wingspan fixed wing nano air vehicle (NAV), in [18], the longitudinal, lateral and coupled models excluding propeller effects computed using CFD tools are presented.
This paper explains in detail the nonlinear six degrees of freedom (6DOF) dynamic modeling of a fixed wing MAV. A very brief explanation of the dynamics of the MAV considered in this paper is presented in [19]. The lift, drag, pitching moment and side force is computed from wind tunnel testing with rotating propeller. The effect of propeller wake in the lift, drag, pitching moment and side force is quantified with explicit equations obtained through curve fitting the wind tunnel test data. The dynamic thrust of motor-propeller is estimated by measuring the axial force with rotating propeller and stationary propeller. Some of the static derivatives and the dynamic derivatives are obtained using empirical relations mentioned in [20]. The trim conditions for straight and constant altitude flight are presented for the entire flight envelope. The linear longitudinal and lateral state space model and corresponding transfer function model are presented. A static output feedback controller is designed for the obtained model and closed loop flight testing is conducted in the presence of wind disturbances. The flight test data shows satisfactory stability characteristics of the MAV flight.
The paper is organized as follows. The nonlinear dynamic model of the MAV is presented in Section 2. Wind tunnel test results are discussed in detail in Section 2. The flight envelope of MAV and trim conditions for straight and constant altitude flight is given in Section 3. Linear longitudinal and lateral state space model of the MAV is presented in Section 4. Section 5 gives controller design details and closed loop flight test results followed by conclusions in Section 6.

Nonlinear Dynamic Model of MAV
This section develops the nonlinear dynamic model of 150 mm MAV, named KH2013A shown in Figure 1 with specifications given in Table 1 [19]. The details of the force and moment coefficients used in this section are given in Appendix A.  The MAV is a flying wing aircraft with autopilot hardware and servo motors for actuating control surfaces placed inside the wing. The airfoil used is modified version of Eppler-387 (E387), suitable for low Reynolds number flyers. The thickness of the conventional E387 airfoil of about 9% is increased to 25% to accommodate the autopilot hardware and other components. The profile of the airfoil used in MAV is shown in Figure 2. The specifications of the airfoil are given in Table 2.  The force equations, moment equations and kinematic equations for position and orientation constitutes the six degree of freedom nonlinear model of the MAV. The equations described in this section is explained in more detail in [21]. LetV = [u, v, w] T denote the velocity components of the MAV along the body axis andP = [x, y, z] T be the position of MAV in inertial axis. Let ω = [p, q, r] T denote the angular velocity along the body axis andĒ = [φ, θ, ψ] T be the Euler angles. The vector form of force equation is given in (1).˙V where F g is the gravitational force, F a is the aerodynamic force and F t is the propulsive force. The gravitational force is given in (2).
where m is the mass of MAV and g is the acceleration due to gravity. Lift and drag forces contribute to the aerodynamic forces acting on MAV and is given in (3).
where X a , Y a and Z a denotes the aerodynamic forces along body X axis, Y axis and Z axis of the MAV respectively. The expression for X a is given in (4).
where ρ is the density of air, V a is the MAV velocity magnitude, S is the wing area andc is the chord length. Similarly the expression for Y a and Z a are given in (5) and (6) respectively.
The propulsive force acts only along body X axis as the thrust line of the propeller is aligned with body X axis. Hence the propulsive force of the MAV is written as given in (7). In this, T denotes the dynamic thrust and is estimated from wind tunnel tests.
The moment balance equations along body X axis, Y axis and Z axis of the MAV is given in where L, M and N are rolling moment, pitching moment and yawing moment respectively and is given in (11)-(13) respectively.
The parameters t 1 to t 8 mentioned in (8)-(10) are constants depending on the moment of inertia of the MAV. The parameters t 1 to t 8 are combinedly represented by a vector t m as given in (14). The expression for t m is given in (15).
The relation between body angular rates (ω) and Euler angles (Ē) is given beloẇĒ where the rotation matrix R 1 (Ē) is given in (18).
The velocity component along the body axis is multiplied by the rotation matrix to yield velocity components along the inertial axis. The relation is given in (19).P = R 2 (Ē)V (19) where the rotation matrix R 2 (Ē) is given in (20).

Estimation of Static Derivatives and Control Derivatives from Wind Tunnel Tests
The wind tunnel tests for the MAV was conducted at National Aerospace Laboratories (NAL). The specifications of the open circuit wind tunnel used for testing is given in Table 3. The MAV undergoing wind tunnel testing is shown in Figure 3. The lift force (L), side force (Y) and aerodynamic pitching moment (M) are measured in wind tunnel with rotating propeller. The propeller used here is GWS 5 × 3 (5-inch diameter and 3-inch pitch) along with Hobbyking AP05 BLDC motor weighing 5.4 g. The specifications of the propeller is given in Table 4 [22]. The static thrust of the propeller with AP05 BLDC motor at an input voltage of 8.4 V is measured using a load cell and is given in Table 5. The thrust force (T) and drag force (D) are estimated from the axial forces measured with rotating propeller and with the stationary propeller. The maximum RPM (rotations per minute) value is 12,000 with airflow. The stall angle of attack (α stall ) is not captured in the wind tunnel tests as the there were provisions for varying the angle of attack from −5 • to +25 • only. The air velocity is varied from 5 m/s to 13 m/s and measurements are taken. This corresponds to the Reynolds number range of 37,500 to 97,500. The lift force obtained with zero elevator deflection and zero propeller RPM is shown in Figure 4. The measurements are taken for an interval of 1 • of α and are marked in all the plots. From Figure 4, we can see that the total lift force generated increases with increase in Reynolds number. The plot of lift generated when the elevator is deflected 25 • upwards is shown in Figure 5. From Figures 4 and 5, we can see that there is a reduction in the lift force when the elevator is deflected upwards as it reduces the effective camber of the airfoil. Figure 6 shows the lift generated for 15 • downward elevator deflection. From Figures 4 and 6, we can see that there is an increase in the lift force when the elevator is deflected downwards as it increases the effective camber of the airfoil. The amount of lift force generated is significantly influenced by the propeller wake as the ratio of the propeller diameter to the wingspan is very high (0.68) when compared to the other UAVs. The plot of the lift force generated for a propeller RPM of 8000 and 12,000 is given in Figures 7 and 8 respectively. While comparing Figures 4 and 8, we can see that there is at least 10% increase in the lift force generated under the influence of propeller wake.     The relation between coefficient of lift (C L ) and total lift force is given in Equation (21).
The coefficient of lift has contributions from wing, elevator and propeller induced air flow as given in (22). In (22), C Lw , C Lδ e and C Lt are coefficient of lift contributed by wing, elevator and propeller air flow, respectively. The coefficient C Lt is scaled for a maximum propeller RPM value of 12,000. The coefficients of C Lδ e are different for upward and downward elevator deflection (δ e >= 0) and upward elevator deflection (δ e < 0) respectively. The coefficients C Lw , C Lδ e and C Lt are given in Table 6. The relation between coefficient of drag (C D ) and total drag force is given in (23).
The coefficient of drag (C D ) has contributions of wing (C Dw ), elevator (C Dδ e ) and due to propeller wake (C Dt ). The total axial force measured with rotating propeller includes drag force and thrust force. The contribution of wing and elevator to drag force is measured with the propeller in a stationary condition. The contribution of propeller wake to drag force is estimated through empirical relations. The plot of C D for zero elevator deflection and the stationary propeller is given in Figure 9. The total drag force increases with increase in Reynolds number. However, C D decreases with increase in Reynolds number as C D ∝ 1 V 2 a for a given drag force. The plot of C D for −25 • elevator deflection and +15 • elevator deflection with the stationary propeller is given in Figures 10 and 11 respectively.   The coefficient of drag of the MAV is given in (24). The coefficient of drag due to wing (C Dw ) is modeled as a function of airspeed, and angle of attack. The expression for C Dw is given in (25) and for C Dδ e is given in (26). The coefficients of C Dw and C Dδ e are given in Table 7 for positive and negative elevator deflection. The contribution of propeller wake to drag, C Dt is taken from [23] and is given in (27).
where S prop is the area swept by one propeller rotation. While taking the measurements with rotating propeller, the total drag force measured is the difference between the MAV drag and the dynamic thrust generated by the propeller. The estimated dynamic thrust, T (δ th ) is given in (28).
where n is the propeller rotation per second (RPS) and d p is the propeller diameter.
The relation between the coefficient of pitching moment C m and pitching moment (M) is given in (29). where In (30), c mw is due to wing (31), C mδ e due to propeller (32) and C mt due to propeller wake (33).
C y is expressed as a function of sideslip angle (β) and rudder deflection (δ r ) and is given in (35).
The derivatives C yβ and C yδ r are given in (36) and (37) respectively. and The other static derivatives C lβ , C nβ and control derivatives C lδ r , C nδ r are obtained through empirical relations given in [20].

Estimation of Dynamic Derivatives
Estimation of dynamic derivatives through wind tunnel tests requires sophisticated instruments like rotating rigs. The wind tunnel does not have the required instrumentation for measuring dynamic derivatives. So the empirical relations provided in [20] is used here for estimating the dynamic derivatives. The dynamic derivatives are a function of α, V a , β and the MAV configuration parameters like specifications of the vertical tail, winglets etc. The estimated static derivatives from wind tunnel tests and dynamic derivatives from empirical relation is given for a velocity of 8 m/s, α = 13.1 • and β = −3.05 • is shown in Table 8. This corresponds to trim conditions as explained in the next section.

MAV Flight Envelope and Trim
The MAV is said to be flying in a trimmed condition when the net forces and moments acting on MAV equals to zero. The upper limit on velocity in the flight envelope of the MAV is determined by the thrust available (T a ) and thrust required for different flight velocities. The lower value of velocity is determined by the maximum safe operating angle of attack. Flying at an angle of attack near stall value is not recommendable as a small disturbance may lead to a stall. From the wind tunnel testing, it is observed that α stall of the MAV is higher than 25 • . Hence a maximum angle of attack of 25 • is taken as the limit for normal flying conditions. The trim conditions for straight and constant altitude flight are obtained by equatingu = 0,v = 0,ẇ = 0,ṗ = 0,q = 0,ṙ = 0,φ = 0,θ = 0 andż = 0. The equations are solved using the fsolve routine in MATLAB. The trim conditions for straight and constant altitude flight are given in Table 9. As seen in Table 9, at V a = 6 m/s, cruise α is slightly lower than 25 • . At V a = 13 m/s, thrust required is almost equals thrust available. So the flight envelope of MAV is determined as in the velocity range of 6-13 m/s. Table 9. Trim conditions of MAV. Note that unlike bigger aircraft, the value of β, φ and δ r is non-zero for MAV for straight and constant altitude flight conditions. This is due to the counter torque exerted by the rotating propeller.The counter torque cannot be neglected for MAVs as the inertia J xx is small when compared to bigger UAV. The effects of counter torque on the dynamics of a small fixed wing MAV is explained in [24].

Linear State Space Model of MAV
The linear state space model is computed by linearzing the nonlinear dynamic model equations of the MAV about trim conditions. The variables used in this section are linearized, for exampleũ denotes the linearized variable for u and so on. The linear state space model equations can be separated into linear longitudinal state space model and linear lateral state space model [25]. The linear longitudinal state space model is given in (38).Ẋ long = A long X long + B long U long where X long = [ũ,w,q,θ] T and U long = [δ e ,δ t ] T . The matrices A long and B long is given in (39) and (40) respectively, for a trim condition corresponding to V a = 8 m/s.
The phugoid mode natural frequency (ω ph ) and damping ratio (ζ ph ) is 1.94 rad/s and 0.283 respectively. The short period mode natural frequency (ω sp ) and damping ratio (ζ sp ) is 35.7 rad/s and 0.246 respectively. The variation of short period mode and phugoid mode natural frequency and damping ratio with respect to flight velocity is shown in Table 10. The ω sp increases with increase in velocity where as ω ph decreases. The ζ sp increases initially with increase in velocity, but later decreases. Whereas, ζ ph decreases initially with increase in velocity, but later increases. The transfer functions  The lateral dynamics consists of Dutch roll mode, roll subsidence mode and spiral mode. The variation of Dutch roll mode natural frequency (ω dr ) and damping ratio (ζ dr ) is given in Table 11. Similar to ω sp , ω dr also increases with increase in flight velocity. The damping ratio, ζ dr decreases with increase in velocity. The roll subsidence mode pole (P rl ) remains nearly constant as shown in Table 11.
The spiral mode pole, (P sl ) is stable for all flight velocity. The transfer functionsṽ In (53) to (56), ∆ lat (s) is the characteristic polynomial for linear lateral dynamics and is given in (57).

Closed Loop Flight Test Results
A static output feedback controller is designed for improving the damping ratios of short period mode, phugoid mode and Dutch roll mode. The controller synthesis methodology is similar to the algorithm explained in [26]. The controller is synthesized in discrete time domain with a sampling time of 20 ms. The MAV linear state space model corresponding to a velocity of 8 m/s is used as the plant model. The structure of the controller is given by where U c (k) = [δ e (k),δ th (k)] T , Y(k) = [q(k),θ(k)] T for longitudinal control and F is a constant feedback gain matrix. For the control of lateral dynamics, U c (k) = [δ r (k)] and Y(k) = [p(k),r(k),φ(k)] T . The procedure adopted in controller synthesis is explained below.
Let the state space model of the open loop plant be given as below.
The control design objective is to find a static output feedback controller of the form given in (58) that minimizes the H ∞ norm of the following mixed sensitivity function.
In (61), T zw d represents the closed loop transfer function from disturbance input to performance output. The sensitivity function, complementary sensitivity function and control input sensitivity function is denoted by S, T and KS respectively. The open loop plant given in (59) and (60) is augmented with frequency dependent weighting functions W 1 , W 2 and W 3 to obtain the generalized plant. The discrete time equivalent of the closed loop generalized plant is given in (62)-(65).
In (63), Z 1 (k) is the output of weighted sensitivity function (W 1 S). Similarly, Z 2 (k) and Z 3 (k) given in (64) and (65) are the outputs of weighted complimentary sensitivity function (W 2 T) and weighted control input sensitivity function (W 3 KS) respectively. The disturbance input to the plant is denoted by W d . The controller is synthesized by solving the following linear matrix inequality (LMI) stated in (66), [27].
In (66), µ < 1 is the radius of the circle in discrete time complex plane, inside which poles of the closed loop system need to placed. The unknowns are the matrix P > 0 and feedback gain matrix F. The elements of the matrix R are selected initially to solve for P and F. The elements of the matrix R are the input variables of genetic algorithm (GA), that minimizes the performance index given in (61). The feedback gain matrix is given in (67) for control of longitudinal dynamics and in (68) for control of lateral dynamics respectively. Please note that the gain from q and θ to δ th is high when compared to δ e . This is because δ e falls in a range of −0.2618 radians (−15 • ) to 0.4363 radians (+25 • ), where as δ th falls between 0-120 RPS (rotations per second).
A comparison between open loop and closed loop damping ratios and natural frequencies of short period mode and phugoid mode for continuous time system are given in Table 12. The short period mode damping ratio is increased from open loop system value of 0.246 to 0.509 for closed loop system. The short period mode frequency of closed loop system is slightly increased when compared to open loop system. The damping ratio ζ sp is improved by feedback of q to δ e . Due to weight and power restrictions in MAV, a low bandwidth actuator is used for elevator and rudder control surfaces. So a very high ζ sp is not practically realizable. The phugoid mode frequency remains almost same for open loop and closed loop system. There is a slight improvement in ζ ph for closed loop system. For achieving a high ζ ph , velocity to δ th feedback is required. Since MAV does not have any light weight and accurate airspeed sensor, a very high ζ ph is not achievable. The open loop and closed loop modes of lateral dynamics for continuous time system are given in Table 13. The Duthch roll mode damping ratio is improved for closed loop system. The ζ dr is improved by r feedback to δ r . Due to limitations in actuator bandwidth, a very high ζ dr is not practically feasible. The flight testing is conducted using a customized 7 g autopilot hardware shown in Figure 12. The autopilot consists of 3-axis accelerometer, 3-axis rate gyro, 3-axis magnetometer, altimeter, GPS and a bi-directional communication module. Servomotors are used for elevator and rudder control surfaces. The thrust is generated by a BLDC motor. The sensors and actuators are interfaced to a 32-bit ARM Cortex-M3 processor, where all the computations are done. Control input is given at a sampling time of 20 ms. The MAV is hand launched at a trim condition corresponding to V a = 8 m/s as in Table 9. Apart from the feedback control inputs, the inputs from the pilot are also given. The flight tests are conducted under mild wind conditions to evaluate the flight stability characteristics. The plot of outputs q, θ and inputs δ e , δ th for longitudinal control are shown in Figure 13. When δ th is high close to 2 seconds, the MAV exhibits pitch-up motion as seen from θ plot. The outputs p, r, φ and input δ r for lateral control is given in Figure 14. The magnitude of p is higher than q and r due to counter torque and cross wind effects. The control inputs shown in Figures 13 and 14 are combined input from pilot and feedback control. The inputs from feedback control and those from pilot are not recorded separately. A photograph of MAV flight is shown in Figure 15. The flight test indicate that the static output feedback controller designed based on the longitudinal and lateral model developed in this paper stabilizes the MAV flight.

Comparison between Predicted p, q, r Data from Model and Flight Test
A comparison between data obtained from flight tests and those predicted by the developed nonlinear model for the angular rates p, q, r is shown in Figures 16-18 respectively. The predicted data is obtained by giving the elevator, rudder and thrust inputs shown in Figures 13 and 14, to the nonlinear equations given in (1) to (18). The discrepancy between flight data and predicted data is due to the following reasons. (1) The wind disturbances acting on the MAV during flight is not measured. So the response of MAV due to wind disturbances cannot be predicted. (2) Due to manufacturing errors, the MAV used for flight testing is not exactly identical to the one used for wind tunnel testing.   The error between flight data and predicted data is given in Table 14. The percentage error P E (v) is computed using (69).
where A(.) denotes the average, v p denotes the predicted value of variable (p, q, r), v f d denotes the flight data of the same variable and M(.) denotes the maximum. The angular rate q depends upon elevator input and thrust input. The angular rate p depends upon rudder input and counter torque generated due to propeller rotation. The angular rate rate r depends upon mainly on rudder input. So the prediction error is higher for p and q when compared to r. This is also due to the fact that the thrust and counter torque model used in the nonlinear model is determined at constant voltage of 8.4 V input to the motor. However, in actual flight, the battery voltage is varying and is not measured. This accounts for an error in predicting dynamic thrust and counter torque during flight. Table 14. Percentage error in prediction.

Conclusions
This paper presented the nonlinear 6DOF dynamic modeling of a 150 mm wingspan fixed wing MAV. The effects of propeller wake on lift, drag and pitching moment are quantified. There is an increase of 10% in lift considering the propeller effects. The propeller wake had increased drag substantially. The longitudinal static stability improves with propeller effects. The dynamic thrust is also estimated from wind tunnel test results.
Unlike a bigger aircraft, the value of sideslip angle, roll angle and rudder deflection is non-zero to maintain a straight and constant altitude flight conditions. This is due to the effect of significant counter torque generated by the rotating propeller. The linear longitudinal and lateral state space models are presented. The H ∞ static output feedback controller improves the damping ratios of various modes of longitudinal and lateral dynamics. The closed loop flight test results with the static output feedback controller indicate satisfactory stability characteristics of the MAV flight. The obtained flight data is compared with the flight data generated from the nonlinear model to verify the modeling accuracy.

Acknowledgments:
The authors acknowledge ARDB-NPMICAV program for providing partial funding for this study. The authors are thankful to National Aerospace Laboratories, Bangalore for helping in conducting the wind tunnel testing.
Author Contributions: Harikumar Kandath designed the MAV and controller. Jinraj. V. Pushpangathan, Titas Bera and Sidhant Dhall performed wind tunnel testing. M. Seetharama Bhat contributed to flight control design, MAV configuration design and analyzing wind tunnel experimental data.

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

Abbreviations
The following abbreviations are used in this manuscript: Appendix A C L -Coefficient of lift C D -Coefficient of drag C m -Pitching moment coefficient of the wing C Lq -Coefficient for change in lift with respect to q C Dq -Coefficient for change in drag with respect to q C y -Side force coefficient C yp -Coefficient for change in side force with respect to p C yr -Coefficient for change in side force with respect to r C yβ -Slope of the side force coefficient of MAV versus β curve C yδ r -Slope of the side force coefficient of MAV versus δ r curve C mq -Coefficient for change in pitching moment with respect to q C l p -Coefficient for change in rolling moment with respect to p C lr -Coefficient for change in rolling moment with respect to r C lβ -Coefficient for change in rolling moment with respect to β C lδ r -Slope of the rolling moment coefficient of MAV versus δ r curve C np -Coefficient for change in yawing moment with respect to p C nr -Coefficient for change in yawing moment with respect to r C nβ -Coefficient for change in yawing moment with respect to β C nδ r -Coefficient for change in yawing moment with respect to δ r