c Distributed Coordination and Control Experiments on a Multi-UAV Testbed by Ellis T. King Bachelor of Engineering The State University of Bu?alo, 2002 Submitted to the Department of Aeronautics and Astronautics in partial ful?llment of the requirements for the degree of Master of Science in Aeronautics and Astronautics at the MASSACHUSETTS INSTITUTE OF TECHNOLOGY September 2004 ? Massachusetts Institute of Technology 2004. Author . .... .... .. .. .... .... .... .... .... .... .... .... . ... ... . .... .... . Department of Aeronautics and Astronautics August 20, 2004 Certi?ed by. ... .... . ... .... .... .... .... .... .... .... ..... .... .. .. .... . Jonathan P. How Associate Professor Thesis Supervisor Accepted by . .... .... .... .... .... .... .... ... . .... .... .... .... .... .... Jaime Peraire Professor of Aeronautics and Astronautics Chair, Committee on Graduate Students Chapter 2 Hardware In the Loop Modeling and Simulation The hardware-in-the-lo op (HWIL) simulation is only useful if it accurately portrays the vehicle dynamics and if the behaviors observed during ?ight tests can be repli- cated on the ground. This chapter focuses on identifying some of the dynamical modes of ?ight for the 60 ARF Trainer aircraft, and verifying that the HWIL sim- ulations re?ect the dynamics expected from the aircraft being employed. Reduced order models for 4 of the 5 dynamical modes are determined for the trainer ARF 60 aircraft using identi?cation techniques on experimental ?ight data and analytical predictions based on aircraft geometry and aerodynamic data. Section 2.1 describes the simulation settings used to create the hardware-in-the-lo op (HWIL) simulations, and Section 2.2 details the procedures used to create models of the aircraft dynam- ics from data collected during ?ight tests and hardware-in-the-lo op simulations. In Section 2.3, the Cloud Cap autopilot is tuned for the trainer ARF 60 aircraft and the closed loop response for several of the modes is measured using the HWIL simulator. 2.1 Hardware in the loop simulations 2.1.1 Aircraft simulation Model Parameters Aerodynamic, inertial and engine calibration information is provided to the Cloud Cap HWIL simulation application to model the aircraft being ?own. For simply 31 (a) The Clark YH airfoil geometry. ?5 0 5 10 15 20 25 30 ?0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 α [deg] C L C D (b) The Clark YH airfoil lift and drag curves. Figure 2-1: The Clark YH airfoil closely resembles the airfoil used on the trainer ARF 60 aircraft and is used to model wing aerodynamics. con?gured aircraft such as the tower trainer 60 ARF used in the testbed, many of the performance characteristics can be obtained using the geometry of the aircraft, such as the data found in Table 2.1. Detailed descriptions of the surface geometry, wing lift curves, and engine performance curves enable simulations of the aircraft under realistic ?ight conditions, providing the input parameters are con?gured accurately. For example, the Clark YH airfoil closely resembles the trainer ARF 60 airfoil and is used to describe the aerodynamic properties of the main wing on the aircraft [25]. Some of the data is shown in Figure 2-1. A more detailed description of the simulator input ?les is given in Ref. [26]. 32 Table 2.1: Trainer 60 ARF measurements, experimentally determined iner- tias as shown in Subsection 2.1.1. Symbolic notation is borrowed from Ref. [27] Measurement Value Units Symbol (SI) Wing Span Wing Area Chord Length Wing Incidence Wing Dihedral Wing Sweep Tail Area Tail Span Tail O?set X Tail Sweep Fin Area Fin Span Fin O?set X Fin O?set Z Fin Sweep Fin Volume Ratio Fuselage CX Area Fuselage Length Gross Mass Empty Mass Roll Inertia* Pitch Inertia* Yaw Inertia* 1.707 0.5200 0.305 1 5 0.0 0.0879 0.606 1.14 9 0.0324 0.216 1.143 0.120 53 0.0189 0.130 1.270 5.267 4.798 0.31 0.46 0.63 m 2 m m deg deg deg 2 m m m deg 2 m m m m deg - 2 m m kg kg kg m 2 · kg m 2 · kg m 2 · b S cˉ Γ Λ Λ l b S t t t t S f b f l f Λ h f f ˉ V f I I I m l S b b m e xx yy zz Aircraft Inertia Experiment The aircraft pitch, roll, and yaw inertias are important parameters for the accurate HWIL simulation of the aircraft dynamics. Fortunately, due to the small scale of the aircraft, experimental measurements can be easily made for each axis of the aircraft. The experimental setup for the roll axis is shown in Figure 2-2. From the aircraft free body diagram, the tension in each cable, T ,is 2T = mg (2.1) 33 Figure 2-2: Torsional pendulum experimental setup to determine roll axis inertia, I xx , for the trainer aircraft. The period of oscillation of a roll angle perturbation, φ, is measured to parameterize the aircraft inertia. The angle α is the small angle deviation of the supporting cables from the vertical position. This experiment was also repeated for the pitch and yaw axes to determine I yy and I zz respectively. For rotational perturbations applied to the airframe, the product of interior angles and distances must be constant Rφ = Lα (2.2) where φ is the aircraft roll angle perturbation and α is the small angle deviation of the supporting cables from the vertical position. The di?erential equation describing the motion of the torsional pendulum is governed by a torsional inertia term and the restoring moment due to tension forces ¨ I xx φ +2TR sin α = 0 (2.3) 34 ? Table 2.2: Experimentally determined aircraft inertias [kg-m 2 ] Aircraft Roll Axis Pitch Axis Yaw Axis No. I xx I yy I zz 1 0.28 2 0.30 3 0.33 Mean 0.30 Std Dev. 0.029 0.46 0.44 0.47 0.42 0.011 0.65 0.61 0.63 0.63 0.021 Using the small angle approximation for α since L ? R, and substituting known values from Eqs. 2.1 and 2.2, Eq. 2.3 reduces to ¨ mgR 2 I xx φ + φ = 0 (2.4) L which is characterized by the undamped natural frequency, ω, and period of oscilla- tion, T p ω = mgR 2 I xx L (2.5) 2π T p = ω (2.6) By ?nding averaged values for the period of oscillation, T p , in each of the pitch, roll, and yaw axes, the inertia about each axis can be approximated. This experiment neglects aerodynamic and other forms of damping as well as the cross-axis inertias (e.g., I xz , I yz ). The experimental results are summarized for each of the axes and three di?erent aircraft in the same con?guration in Table 2.2, showing agreement between di?erent vehicles used in the tests. The largest variation was found in the roll axis due to the di?culty of mounting the aircraft through the center of gravity, which is essential in this experimental setup. 2.1.2 Actuator Models The servos used on the aircraft have saturation limits, limited bandwidth, and limited slew rates which are captured in the actuation models of the Cloud Cap hardware- 35 Figure 2-3: Actuator models used the Cloud Cap Hardware-in-the-lo op sim- ulations in-the-loop simulator. As shown in Ref. [26], the actuator transfer function, G act (s), is given by specifying the bandwidth limit, B W G act (s)= ω 2 n s 2 +2ζω n s + ω 2 n (2.7) ω n =2πB W (2.8) ζ = ζ c =0.707 (2.9) where the damping ratio, ζ, is selected at the critical value to set the actuator band- width equal to the natural frequency (ω b = ω n ) . The aileron, elevator, and rudder channels all respond with approximately the same characteristics (B W = 10 Hz), but the throttle is modeled with less dynamic range (B W = 2 Hz) as the engine RPM requires added time to ramp up to produce thrust. The input/output saturation and slew rate limits are determined as per manufacturer speci?cations (±60 ? ,2 Hz respectively), and applied as shown in Figure 2-3. 2.1.3 Sensor Noises The Cloud Cap hardware-in-the-lo op simulator includes detailed sensor models based on information from the manufacturer to corrupt the simulation measurements. For the purposes of simulation, noises on the pressure, rate gyros and accelerometers onboard the aircraft are modeled using band-limited white noise and speci?ed drift rates [26]. Although the same noise and drift models could be applied to GPS position and velocity measurements, this information is typically assumed to be perfect in the HWIL tests. The values used to parameterize the Piccolo TM pressure sensors, the Crista TM IMU angle-rate sensors, and the accelerometers are shown in Table 2.3. 36 Table 2.3: Crista IMU HWIL Sensor Noise Models Sensor PDynamic PStatic Gyro Accel. [unit] [Pa] [Pa] [deg/s] [m/s/s] Resolution [unit] Min [unit] Max [unit] Noise Gain Butterworth Order BW Cuto? Freq. [Hz] Drift Rate [unit/s] Max Drift value [unit] Drift Update Rate [s] 3.906 -300 4000 20.0 2 11.0 0.05 15.0 5.0 2.00 0.0 110,000 20.0 2 11.0 1.0 100.0 5.0 1.6E-4 6.0E-3 -5.20 -100.0 5.20 100.0 0.10 0.0 2 2 20.0 20.0 1.5E-4 2.0E-3 0.01 0.20 1.0 1.0 2.1.4 Dryden Turbulence Stochastic turbulence disturbances are required for accurate HWIL simulations, as real world experiments are characterized by unpredictable winds acting on the vehicle. The Dryden turbulence model is one of the accepted methods for including turbulence in aircraft simulations [28]. By applying shaped noise with known spectral properties as velocity and angle rate perturbations to the body axes of the vehicle, the e?ect of turbulence is captured during discrete time simulations. The noise spectrum for each of the perturbations is predominantly described by a turbulence scale length parameter, L, the airspeed reference velocity, V o , and the turbulence intensity, σ. The selection of these parameters allows for the turbulence to be modeled according to the prevailing wind conditions. The spectral frequency content for generalized aircraft turbulence have been well studied [29, 28] and are given for each of the aircraft body axes: S ug (ω)= 2σ 2 u L u πV o · 1 1+ ( L u V o ω) 2 (2.10) S vg (ω)= σ 2 v L v πV o · 1+ 3( L v V o ω) 2 ? 1+ ( L v V o ω) 2 ? 2 (2.11) S wg (ω)= σ 2 w L w πV o · 1+ 3( L w V o ω) 2 ? 1+ ( L w V o ω) 2 ? 2 (2.12) 37 S σ 2 0.8 ? πL w ? 1/3 pg (ω)= w 4b ? 2 (2.13) L w V o · 1+ ? 4b ω πV o ?? ω 2 V o S qg (ω)= ? 2 · S wg (ω) (2.14) 1+ ? 4b ω πV o ? ω ? 2 V o S rg (ω)= ? 2 · S vg (ω) (2.15) 1+ ? 3b ω πV o where ω is the spectral frequency of the turbulence and the aircraft wingspan, b,is used as a parameter in the angle rate ?lters to scale the e?ect of rotation on the main lifting surface. The subscripts u, v, w and p, q, r refer to the familiar body frame aircraft wind velocities and angle rates, respectively, thereby allowing independent classi?cation of the turbulence in each axis. These spectral shaping functions are used to form shaping ?lters to give the body axis noise transfer functions [30] ? L u 1 H ug (s)= σ u 2 1+ L u (2.16) πV o · s 1+ √ 3 L v V o s ? L v V o ?? 1+ L v 2 (2.17)H vg (s)= σ v πV o · s V o s ? L w 1+ √ 3 L w V o ?? 1+ L w 2 (2.18)H wg (s)= σ w πV o · s V o ? 1/6 ? 0.8 ? π 4b (2.19) 1/3 H pg (s)= σ w V o L w ? 1+ ? 4b s ?? πV o s H wg (s) (2.20)H qg (s)= ? V o 4b ? · 1+ s πV o s H vg (s) (2.21)H rg (s)= ? V o 3b ? · 1+ s πV o The block diagram for the full 6 DOF Dryden turbulence model is shown in Figure 2-4. Note the cross axis couplings of the angle rate ?lters q g and r g . Example turbulence perturbations values are plotted in Figure 2-5 as a function of the scale lengths and intensities for each of the body axes. The same 4×1 white noise input was used for each trial set. Larger scale lengths, L, increase the time constant of 38 Figure 2-4: Block diagram for the 6 DOF Dryden Turbulence model. The velocity perturbations u g ,v g ,w g are independent outputs of the ?ltered values of the turbulence scale lengths, L, intensity values, σ and the white noise input sources. The principle axis coupling of the aircraft is taken into account through the inputs to the angle rate perturbation ?lters. the turbulence seen for a given airspeed, while larger σ values increase the deviation about zero. L and σ typically vary with altitude in the lower atmosphere as ground e?ects become more prominent, but for HWIL simulations they are usually ?xed. The frequency response of the Dryden ?lters are shown for the same three cases in Figure 2-6. The ?lter cuto? frequency is determined by the ratio of the scale length to airspeed, and this e?ectively produces lower bandwidth ?lters for larger scale lengths. The scale length parameter is chosen according to one of several speci?cations, all of which take into account the variation of L with altitude. The military reference MIL-F-8785C provides one such model of the scale length at low altitudes, h, which 39 0.5 0 ?0.5 Vu Vv 0 1 2 3 4 5 6 7 8 9 1 0 ?1 0 1 2 3 4 5 6 7 8 9 ?1 ?0.5 0 0.5 Vu Set 1 Set 2 Set 3 Qr Qq Qp 0 1 2 3 4 5 6 7 8 9 Time [s] (a) Velocity perturbations 0.4 0.2 0 ?0.2 ?0.4 0 1 2 3 4 5 6 7 8 9 0.2 0 ?0.2 0 1 2 3 4 5 6 7 8 9 0.1 0 ?0.1 ?0.2 0 1 2 3 4 5 6 7 8 9 Time [s] (b) Angle rate perturbations Figure 2-5: The output of the Dryden velocity and angle rate ?lters for di?erent selections of the intensity and scale lengths. Set 1: L = 150,σ =0.5, Set 2: L = 1500, σ =0.5, Set 3: L = 150, σ =1.5 40 is valid up to 1000 feet [29]. L w = h (2.22) h L u = L v = (2.23) (0.177 + 0.000823h) 1.2 The turbulence intensity is a gain factor that scales the magnitude plots in Fig- ure 2-6 to values appropriate for di?erent wind levels (i.e., light, moderate, severe). The intensitylevel has been de?ned for low altitude ?ight according to MIL-F-8785C as σ σ w =0.1W 20 (2.24) w σ u = σ v = (2.25) (0.177 + 0.000823h) 0.4 where W 20 is the wind speed as measured at 20 ft in altitude. According to MIL- F-8785C, W 20 < 15 knots is classi?ed as “light” turbulence, W 20 30 knots is≈ “moderate”, and W 20 > 45 knots is “heavy”. Other military speci?cations such as MIL-HDBK-1797 exist for the low altitude cases [29], and di?erent types of models are used for other regions of the atmosphere. For the purposes of the UAV application, the low altitude models are su?cient. The utility of the Dryden turbulence model is that it allows the expected turbu- lence levels to be described for an aircraft ?ying at a given reference speed for more realistic HWIL simulations. Turbulence is applied to the vehicle body axes consistent with the known parameterized values for scale length and intensity, which e?ectively de?nes the appropriate ?lters with cuto? frequencies and magnitudes needed for sim- ulation. Note that in addition to turbulence, wind is also usually modeled with a static component, W , that represents a prevailing magnitude and direction in an inertial axis. Together these de?ne an arbitrary three-axis wind vector W = W + δW I (2.26) where δW I is the e?ective Dryden wind turbulence in each axis after being rotated through the appropriate body to inertial transformation direction cosine matrix. For small scale aircraft, the static wind component is usually a gross disturbance relative to the aircraft airspeed, and it can have a large e?ect on high level planning 41 Bode Diagram 20 10 0 ?10 ?20 Phase (deg) Magnitude (dB) Phase (deg) Magnitude (dB) ?30 ?40 ?50 ?60 0 ?45 ?90 Set 1 Set 2 Set 3 ?4 ?3 ?2 ?1 0 1 10 10 10 10 10 10 Frequency (rad/sec) (a) Velocity ?lter for w g Bode Diagram 0 ?20 ?40 ?60 ?80 ?100 ?120 ?140 ?160 ?180 ?200 0 ?45 ?90 ?135 ?180 ?3 ?2 ?1 0 1 2 3 10 10 10 10 10 10 10 Frequency (rad/sec) (b) Angle rate ?lter for q g Figure 2-6: The Bode plots of the Dryden velocity and angle rate ?lters, given white noise inputs to H vg (s) in (a) and H qg (s) H wg (s)· in (b). Various selections of the intensity and scale lengths are shown in di?erent sets. Set 1: L = 150,σ =0.5, Set 2: L = 1500, σ =0.5, Set 3: L = 150, σ =1.5 42 algorithms. The e?ect of this type of disturbance on the planning system and aircraft dynamics is discussed in more detail in Chapters 4 and 5. 2.2 Open Loop Aircraft Modeling 2.2.1 Longitudinal Dynamics A common model for the longitudinal motion of the aircraft is [27, 31] ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? u˙ x x x x u x θ δ e u w q w˙? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? z z z z w z θ δ e u w q δ e (2.27)+ = q˙ θ ˙ m m m m q θ m θ δ e u w q ? ? ? ?? ? ? ? 0 0 1 0 0 x˙ = Ax + Bu (2.28) where the state variables x =[uwqθ] T refer to the longitudinal velocities, u and w, the pitch rate, q and the angle of inclination, θ. The elements of the A matrix in Eq. 2.27 represent the concise form aerodynamic stability derivatives referring to the airplane body axis. Tables of values relating the concise form derivatives to the dimensionless or dimensional derivatives are available in numerous sources [27, 32]. The control input u = δ e is the elevator defection angle with the engine thrust ?xed and is input to the dynamics through the aerodynamic control derivative matrix B. The Longitudinal Dynamics in Eq. 2.27 are typically resolved into two distinct phugoid and short period modes, which represent dynamics of the aircraft on di?erent timescales. The short period is characterized by high frequency pitch rate oscillations and can have high or low damping, depending on the dynamic stability of the aircraft. In contrast, the phugoid mode is characterized by lightly damped, low frequency oscillations in altitude and airspeed with pitch angle rates, q, remaining small. Short Period Mode A simple approximation for the short period mode of the aircraft can be obtained by assuming the speed of the aircraft is constant over the timescale of the short period dynamics (˙u = 0), that the aircraft is initially in steady level ?ight and that the 43 derivatives refer to a wind-axis system (θ = α = 0). The equations of motion then reduce to ? w˙ ? ? z w z q ?? w ? ? z δ e ? ? δ e (2.29) ? ? = ? ?? ? + ? q˙ m w m q q m δ e Following further approximations shown in [27] which make assumptions about the relative size of the m q ,z q and z w derivatives, the transfer functions for the two short term equations describing the response to elevator are: w(s) z n ? s+ V o m δ e ? ? k w (s+1/T α ) (2.30) z δ e = δ e (s) (s 2 ?(m q + z w )s+(m q z w ?m w V o )) s 2 +2ζ s ω s s+ ω s 2 q(s) m n (s?z w ) ? k q (s+1/T θ ) (2.31)= δ e (s) (s 2 ?(m q + z w )s+(m q z w ?m w V o )) s 2 +2ζ s ω s s+ ω s 2 where k q ,k w ,T θ ,T α ,ζ s , and ω s represent approximate values for the short period mode and V o is the vehicle reference speed. One of the most accurate ways to obtain models for the aircraft data is to use actual ?ight data. Identi?cation algorithms such as those in the Matlab System Iden- ti?cation Toolbox [33] can be used to used to obtain open loop models of the system dynamics from ?ight data collected during experiments. These models can then be used to validate the HWIL simulation environment as well as to help determine the gain settings for the autopilot control loops as shown in Section 2.3. Input-output data was collected by disengaging all of the autopilot loops and performing a series of maneuvers to measure the aircraft response to de?ections from the elevator and aileron control surfaces. Example data from two experiments are shown in Figures 2-7(a)-(b) depicting the longitudinal and lateral modes, respectively. To capture the longitudinal dynamics, the bank angle was held ?xed at zero degrees, while a series of pitch oscillations were commanded using the elevator. Figure 2-7(a) shows a sample of data that was collected on one run of the pitch test. From the plot it is clear that the longitudinal modes are being excited due to input from the elevator, while the lateral motions in the roll and yaw axes are essentially ?xed. Sample data from a roll excitation run is plotted in Figure 2-7(b). This plot also clearly shows coupling in the yaw axis due to the dihedral angle of the wing. 44 (a) Pitch test data induced with elevator de?ections. 100 1055 1060 1065 1070 1075 1080 1085 ?100 ?50 0 50 Roll Axis [deg/s] Aile Rudd Elev Yaw Axis [deg/s] 40 20 0 ?20 ?40 1055 1060 1065 1070 1075 1080 1085 0 0 100 ?100 Pitch Axis [deg/s] ?100 1055 1060 1065 1070 1075 1080 1085 time [s] 306 308 310 312 314 316 318 ?100 ?50 0 50 100 Roll Axis [deg/s] ?100 ?50 0 50 Yaw Axis [deg/s] Aile Rudd Elev 306 308 310 312 314 316 318 ?20 ?10 0 10 20 Pitch Axis [deg/s] 306 308 310 312 314 316 318 Time [s] (b) Roll test data induced with aileron de?ections. Note the yaw axis coupling due to the large dihedral angle of the wing Figure 2-7: Sample ?ight test data used in the estimation algorithms. Dash- dot lines represent rate gyro data output for each axis. Control surface inputs are plotted for each corresponding axis in degrees of de?ection 45 Parametric models for the input-output transfer functions can be formed using experimentally collected data and used to determine the unknown coe?cients in Eq. 2.31. Then from the characteristic equation, the longitudinal short period dy- namics can be inferred from models of the transfer function from elevator angle to pitch rate. Also note that qualitative predictions about the values of the parameters in Eq. 2.31 are available, since they depend on the concise stability derivatives which all have physical signi?cance. Once models are obtained the approximate parameters can then be veri?ed against these qualitative predictions. Figure 2-8 shows the output of a parametric subspace model based on experimental data such as that shown in Figure 2-7(a). The model output (solid line) tracks the actual measurements of pitch rate (dashed line) quite well and was validated on data sets from di?erent test days and aircraft. Model residuals within the 99% con?dence bounds for the auto- and cross correlation functions are plotted in Figure 2-9 and indicate that the 3 rd order model is su?cient to describe the input-output dynamics. This model is represented by the continuous transfer function in Eq. 2.32 and has zero-pole pairs as indicated in Figure 2-10(a). The short period is well represented by the high frequency, oscillatory mode of the system, while there is one low frequency pole located near the imaginary axis. The transfer function is: q(s) 9.539s 2 ?1440s +60.52 T qδ e (s)= = δ e (s) s 3 +21.77s 2 + 325.8s +29.94 9.539(s ?150.9)(s ?0.0420) (2.32)= (s 2 +21.7s + 323.7)(s +0.0925) ζ The third order model T qδ e (s) was selected because it provided the best ?ts to a large number of data sets and residuals that remained below the 99% con?dence intervals in Figure 2-9. Approximate values for the short period dynamics can be obtained by resolving the oscillatory mode in Eq. 2.32 to determine the corresponding values for ω s and ζ s . A step response for this mode is plotted in Figure 2-10(b), indicating a reasonable short period response time with settling time 0.4 seconds, and s =0.6. Second order models produced from the same data set were found to have di?culty reproducing the outputs of the experiment, and exceeded the con?dence bounds as shown in Figure 2-9. Higher order models (4 th and higher) tended to 46 place additional pole-zero pairs near the origin without achieving better tracking or residual bound performance, therefore indicating the 3 rd order model as the best representation for the system. The zero locations in Eq. 2.32 are not consistent with the expected values for the stability derivatives in Eq. 2.31, which is an indication of the delay acting on this input-output channel and consistent with the servo response time in an actual physical system. In addition, it should be noted that the zero location for these types of models is generally more uncertain than the pole location. As a result, the zeros in Eq. 2.31, are di?cult to identify without more sophisticated validation techniques. In order to validate the HWIL model the same set of tests was performed in simulation and the results are compared in Figure 2-11. As shown in Figure 2-11(a), both models reproduce the experimental outputs with excellent tracking. Since the data set used to validate these two models is di?erent from either of the sets used to generate them, this model can be relied upon with much higher con?dence. The transient responses are shown in Figure 2-11(b) and the fast acting short period mode is shown to agree well, however there is variation on the longer timescales. The HWIL model for pitch response is given by q(s) 2.232s 2 ?1265s +2.449 T qδ e (s)= = δ e (s) s 3 +19.2s 2 + 283.7s ?2.3 0.0022(s ?566.92)(s ?0.0019) (2.33)= (s 2 +19.209s + 283.8924)(s +0.0081) which has a short period characteristic equation similar to Eq. 2.32, identi?ed using experimental ?ight test data. The discrete P-Z plot in Figure 2-12 shows the pole and zero locations for the longitudinal models. The short period mode is shown to agree well, however the slower frequency dynamics are not as well modeled. Phugoid Mode The data from the preceding experiments captures the short period mode of the system well, however the phugoid mode is not well represented. This is physically consistent because the phugoid mode excites the airspeed and pitch response over longer timescales, speci?cally when q˙ =˙w ≈ 0. Setting the corresponding terms to 47 ?1 ?0.5 0 0.5 1 [rad/s] Sim Output Exp Data 236 238 240 242 244 246 248 250 252 254 256 Time Figure 2-8: Simulated and measured outputs from elevator to pitch rate. Val- idation on a di?erent data set than that used to create the model. Autocorrelation of residuals for output q 0.1 0.05 0 ?0.05 ?0.1 ?20 ?15 ?10 ?5 0 5 10 15 20 Cross correlation for input δ and output q residuals e ?0.2 ?0.1 0 0.1 2 nd Order 3 rd Order ?15 ?10 ?5 0 5 10 15 20 Samples Figure 2-9: Auto- and cross-correlations with 99% con?dence intervals (dashed lines) for the elevator to pitch rate models obtained from experimental data. 2 nd order model has high cross correlations for positive sample delays, indicating unmodeled higher order dynamics, but 3 rd order models stay within the 99% bounds. 48 10 15 System: untitled1 Pole : ?10.8 + 14.4i Damping: 0.602 Overshoot (%): 9.34 Frequency (rad/sec): 18 System: untitled1 Zero : 0.042 5 Damping: ?1 Overshoot (%): 0 Frequency (rad/sec): 0.042 0 System: untitled1 Pole : ?0.0925 System: untitled1 Zero : 151 ?5 Damping: 1 Overshoot (%): 0 Frequency (rad/sec): 0.0925 Damping: ?1 Overshoot (%): 0 Frequency (rad/sec): 151 ?10 ?15 ?20 0 20 40 60 80 100 120 140 160 (a) The open loop pole and zero locations indicate the short period dynamics identi?ed by the estimation algorithms. Step Response 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Amplitude Time (sec) (b) Short period approximation step input Figure 2-10: Estimated longitudinal model, T qδ e (s), from elevator input to pitch rate. 49 Measured and simulated model output ?? Pitch Axis ?1 ?0.5 0 0.5 1 HWIL Sim Model Flight Data Exp Flight Model Pitch Axis [rad/s] ?1.5 188 190 192 194 196 198 Time [s] (a) Both models created from the actual ?ight data and the HWIL simulator are able to track the experimental data well. ?5 ?4.5 ?4 ?3.5 ?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0 Pitch Rate Response [rad/s] HWIL Model Flight Model 0 0.5 1 1.5 2 Time [s] (b) Transient response for the longitudinal models. Figure 2-11: Transient response for the 3 rd order models are agree in time constant and damping, showing agreement in the short period dynamics of the HWIL and ?ight test models. 50 0.6 0.5 0.4 0.3 0.2 0.1 0 ?0.1 0.3π/T 0.2π/T 0.1π/T 0.9 0.3π/T 0.2π/T 0.1π/T 0.4π/T 0.3 0.1 0.4π/T 0.2 0.8 0.4 0.5 0.6 0.7 HWIL Poles HWIL Zeros Flight Poles Flight Zeros 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Figure 2-12: Discrete P-Z plot for the experimental and HWIL sim models. Damping and time constant of short period mode agree to a certain margin, but the slow dynamics are not captured well in this experiment. zero in Eq. 2.27, the dynamics can be simpli?ed to a second order system involving only the body x-acceleration, u˙, and the pitch rate, θ ˙ [27] ? u˙ ? ? ? m u V o ?m q z u ? ?g ?? ? ? ? m δ e V o ?m q z δ e ? ? u ? θ ˙ ? = ? x u ?x w m w V o ?m q z w ?? ? + ? x δ e ? m w V o ?m q z w ? δ e ? m u z w ?m w z u ? ? m δ e z w ?m w z δ e ? 0 θ m w V o ?m q z w m w V o ?m q z w x˙ p = A p x p + Bu (2.34) The phugoid mode can then be approximated by ?nding the poles of Eq. 2.34, Δ(s)= sI ?A p = s 2 +2ζ p ω p + ω 2 p | ? | ? m u V o ?m q z u ?? ? m u z w ?m u z u ? 2 = s ? x u ?x w s + g m w V o ?m q z w (2.35) m u z w ?m u z u m For conventional aircraft in subsonic ?ight several approximations can be made [27] u ≈ 0; m w z u |; m w (2.36)|m u z w |?| | V o |?|m q z w | 51 These simpli?cations yield closed-form expressions for the damping and natural fre- quency of the mode x u (2.37)2ζ p ω p ≈? gz u ? ? (2.38)ω p ≈ V o As outlined in [27], the concise stability derivatives can be transformed into their dimensionless equivalents and related to the aerodynamic parameters of C L and C D . The phugoid frequency and damping can then be expressed in terms of the aerody- namics and airspeed, assuming lift is equal to the weight in the trimmed condition g √ 2 (2.39)ω p ≈ V o C D total C D w + C D b + C D p (2.40)ζ p ≈ √ 2C L w ≈ √ 2C L w The subscripts w, b and p denote the contributions due to the wing, body and par- asitic forms of the drag for the aircraft, respectively. Eq. 2.39 is known to provide reasonable approximations for the phugoid natural frequency, but due to the simpli- fying assumptions and uncertainty in the lift and drag coe?cients, Eq. 2.40 is less accurate in determining the damping ratio, ζ p . V In order to observe the phugoid dynamics, an experiment under trimmed condi- tions with airspeed and/or pitch angle dynamics is required. Unfortunately due to space constraints, an experiment requiring long term “hands o?” dynamics on the actual aircraft is not feasible, however this test can easily be performed on the HWIL simulator. The results of a 3 m/s airspeed impulse response with initial airspeed, o = 28 m/s and zero angle of attack are shown in Figure 2-13. The best ?t curve for the oscillation was found to have damping and natural frequency ω p(?t) =0.0732 rad/s and ζ p(?t) =0.20 (2.41) which are close to the predictions from the phugoid approximations in Eqs. 2.39 and 2.40 with predicable uncertainty in the damping term ω p(theory) =0.0786 rad/s and ζ p(theory) =0.14 (2.42) 52 ?2.5 ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 Δ Airspeed [m/s] HWIL Exp data Best Fit Theory 14.3 s 28.6 s 41.0 s ω P(theory) = 0.0786ζ P(theory) = 0.14 ω P(fit) = 0.0732 ζ P(fit) = 0.20 ?10 0 10 20 30 40 50 60 70 Time [s] Figure 2-13: HWIL simulation of the phugoid mode using an airspeed per- turbation con?rms analytical predictions from the initial speed, V o = 28 m/s and Lift/Drag curves. C L w and C D w are deter- mined from the wing lift curve at α = 0. The values calculated for the body and parasitic components of drag are C D b =0.021 and C D p =0.01, respectively. For these values, the wing lift and drag, C L w and C D w terms were determined from the wing lift curves with α = 0 in Figure 2-1(b), while the body drag coe?cient can be estimated from the power curve at the operating speed and the fuselage geometry. 2.2.2 Lateral Dynamics Similar to the longitudinal model, the dynamics describing lateral perturbations about an equilibrium trim condition can be written in concise state space form [27] ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? v˙ p˙ y v y p y r y φ l v l p l r l φ v p y δ a y δ r l δ a l δ r n δ a n δ r ? δ a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? (2.43) + = ? ? ? r˙ φ ˙ ? ? ? ? ? ? n v n p n r n φ 0 1 0 0 ? ? ? ? ? ? r φ ? ? ? ? ? ? δ r 0 0 53 where v is the sideslip velocity, p and r represent the roll and yaw rates, and φ is the roll angle. The aileron control de?ection angle is denoted by δ a , and since the e?ect of the rudder cannot be neglected, it provides an additional control input δ r . Roll Subsidence Mode While the two longitudinal modes decouple rather easily, this is not the case in lateral aircraft dynamics. The lateral motion is characterized by one oscillatory mode known as the Dutch roll mode, and two ?rst-order lags known as the roll subsidence and spiral modes. Because all of the lateral motion is coupled to some extent, it makes identifying individual modes more di?cult. Simplifying approximations can still be made for small roll angles, since there is little yaw or sideslip induced, particularly for fast dynamics. In this case the lateral-directional model can be reduced to [27] ? ? ? ?? ? ? ?? ? p˙ l p l φ p l δ a l δ r δ e (2.44) ? = ? ?? ? + ? ?? ?? φ ˙ 10 φ 0 0 δ a If aircraft wind axes are assumed (l φ = 0), then Eq. 2.44 further simpli?es to a ?rst- order equation representing the roll subsidence mode p(s) k p T pδ a (s)= = (2.45) δ a (s) s + T 1 r Similar to the longitudinal modes, data can be sampled in the roll and yaw axes and transfer function models can be formed to represent the input-output relations. The ?rst-order dynamics in Eq. 2.45 represent the initial response of the system after the ailerons are actuated. Figure 2-14 presents the experimental and simulated output for the roll axis – the ?rst-order model (solid line) is shown to track the output (dashed line) well. The output of the model in Figure 2-14 is represented by the ?rst order system 106.48 T pδ a (s)= (2.46) s +14.36 where the time constant T r = 14 1 .36 =0.070 sec. identi?es the dominant physical properties for the roll mode. Figure 2-15(a) demonstrates the tracking of two models for the roll subsidence mode created from experimental ?ight test data and the same test performed in HWIL 54 1 [rad/s] 0.5 0 ?0.5 ?1 ?1.5 100 102 104 106 108 110 112 114 116 Sim Output Exp Data Time Figure 2-14: Simulated and measured outputs from aileron to roll rate. tests. The data sets show strong agreement indicating that the HWIL simulation provides an accurate representation of the testbed aircraft in roll response. As shown in Figure 2-15(b), the time constant value for both models was found to be roughly 0.1 seconds, which is reasonable for aircraft of this size. Subtle model di?erences created a discrepancy in the overall gain factor of the roll-rate transfer function, but this was corrected in the HWIL simulator by adjusting the aileron chord length from its true physical value until the transfer functions responses were in better agreement. The roll subsidence transfer function for the HWIL simulator is 108.10 T pδ a (s)= (2.47) s +14.27 which has time constant of T r H W = 14 1 .27 =0.07 sec. This model agrees well with Eq. 2.45, indicating strong correlation with the true roll response of the aircraft. Spiral Mode The roll subsidence mode can be simply decoupled from the lateral dynamics, how- ever the spiral and Dutch roll modes are more di?cult to identify. The spiral mode 55 Measured and simulated model output ?? Roll Axis ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 Roll Axis [rad/s] Experimental Model Flight Data Simulated HWIL Model ?2.5 ?3 14 16 18 20 22 24 26 28 30 Time [s] (a) Models created from both the actual ?ight data and HWIL simulator data can track the experimental data. Step response roll subsidence mode 0 1 2 3 4 5 6 7 8 9 10 Roll Rate [deg/s] Flight Model HWIL Model (aile chord 18%) HWIL Model (aile chord 100%) 0 0.2 0.4 0.6 0.8 1 Time [s] (b) Transient response for roll subsidence models. Time constant is ≈0.1 sec. Aileron chord length was varied from true value to obtain better agreement with model obtained from ?ight test data. Figure 2-15: Agreement between the models from aileron input to roll rate output for the HWIL and ?ight data models. 56 is a non-oscillatory, slowly acting mode that captures complex motions in roll, yaw, and sideslip. It is characterized by the interaction between directional stability (also known as weathercock or ?n stability) and the lateral stability (dihedral) of the air- craft. These two e?ects counterbalance one another resulting in the slowly acting spiral mode, which can be stable, neutrally stable or unstable depending on the rel- ative strengths of the lateral directional e?ects. Over long periods, the motion variables v, q, and p can be assumed steady, and Eq. 2.43 can be reduced using the approximations v˙ =˙q =˙ 0. This allows thep ≈ reduced order model for the spiral mode to be formed as shown in [27] ? ? ? (l v n p ?l p n v ) ?? p ? 0 (2.48) ? φ ˙ ? = ? y r l r n v 1 ?l v n r ) y 0 φ ?? φ ? Since φ ˙ = p, Eq. 2.48 can be reduced to a ?rst order di?erential equation describing the unforced spiral mode dynamics. With reference to Appendix 1 of [27], the time constant of the mode, T s can be expressed in terms of the dimensionless derivatives T s = V o (L v N p ?L p N v ) (2.49)? g(L r N v ?L v N r ) Due to the timescales over which the spiral mode acts, it typically is not possible to identify the mode using input-output estimation techniques. Instead, a series of analytical approximations based on aircraft geometry and aerodynamic data can be used to identify the parameters of the spiral mode, and this prediction is used to validate the HWIL settings. Approximations for the derivatives in Eq. 2.49 were found in terms of the aerodynamic and geometric properties of the ARF 60 aircraft [27] ˉ N v = a 1F V f (2.50) 1 ? N p = (2.51)? 12 C L ? ? ? α=0 1 ? l f N r =(N r ) wing +(N r ) fin = ? 6 C D ? ? ? α=0 ?N v (2.52) b 1 L p = (C L α + C D | α=0 ) (2.53)? 12 1 ˉ h f l L v =(L v ) wing +(L v ) fin = ? 4 C L α Γ ?a 1F V f (2.54) f 57 Table 2.4: Dimensionless derivatives for the ARF 60 aircraft contributing to the lateral spiral mode Derivative Approx. ARF 60 N v N p N r L p L v L r Yaw Sti?ness Yaw Moment due to Roll Rate Yaw Damping Roll Damping Dihedral E?ect Roll Moment due to Yaw Rate 0.0792 -0.0167 -0.0547 -0.3550 -0.1103 0.0389 1 l f L r = C L α=0 ?(L v ) fin (2.55) a 6 | b 1F denotes the wing lift curve slope for the tail, which is analytically approximated as a thin wing in subsonic ?ow [32] 2πAR f a 1F =(C L α ) fin = (2.56) 2+ AR f ? 1 + tan 2 Λ f +4 where AR f is the aspect ratio of the ?n. Using the measurements of the aircraft geometry from Table 2.1, the derivatives for the ARF 60 aircraft are computed and summarized in Table 2.4. The time constant of the spiral mode can then be computed from Eq. 2.49 as T s =28. 95. From this rough analysis it is con?rmed that the spiral mode will be a very slow, but stable mode with the positive root given by Δ s (s )=(s +0. 0345) (2.57) This mode was di?cult to observe in ?ight tests or HWIL simulations, but it was con?rmed that the trimmed open loop ARF 60 model does not exhibit divergent behavior over extended periods of time (> 10 min), showing that the low frequency modes are stable as predicted. 2.3 Autopilot Tuning The Cloud Cap autopilot has a prede?ned set of controllers that need to be tuned to suit the 60 ARF aircraft modeled in the previous sections. The HWIL simulator is a 58 useful tool that can be used to tune the autopilot gains and control settings before ?ight testing. Previous sections have veri?ed the accuracy of the HWIL simulations, it is expected these predictions will map well to the actual aircraft being ?own. This section deals with the tuning of the autopilot to suit the requirements of the MILP tra- jectory planner, which pushes the limits of the waypoint tracking strategy employed. By uploading closely spaced waypoints and making frequent changes to the plan the autopilot is tracking, the closed loop system is required to be responsive, however as with any control system, care needs to be taken to protect against instability. In addition to path following, the closed loop tracker is required to reject any disturbances acting on the system in order to ?y the tra jectory plans designed. Later sections incorporate methods for estimating the disturbance levels and accounting for them on higher planning levels, however these estimates will be imperfect, requiring that the autopilot be successful in rejecting bounded levels of wind and measurement noise error. Subsection 2.3.1 addresses the lateral autopilot controllers, the waypoint tracker is described in Subsection 2.3.2, and Subsection 2.3.3 describes the process in tuning the airspeed/altitude controllers. All of the tunings completed using the HWIL simulator are also validated during actual ?ight testing of the ARF 60 aircraft. 2.3.1 Lateral Autopilot E?ective strategies for autonomous control are often accomplished by nesting numer- ous control loops to build up to a complicated desired level of functionality. However the overall performance of the system is then limited by the bandwidth and perfor- mance of the innermost loops. Because the autopilot will eventually be used for tight tra jectory following, the tuning of the inner loops is important to gain adequate sys- tem performance. Quick, tight response is needed on the turn rate loops on the order of the roll dynamics response (from Subsection 2.2.2, T r =14.27), however settings that are too aggressive quickly lead to roll oscillations and eventual instability. 59 Figure 2-16: Lateral Autopilot block diagram Turn Rate The turn rate of the aircraft is regulated by maintaining a desired bank angle, giving an acceleration in the radial direction as shown in Figure 4-3. Assuming a coordinated turn, the desired bank angle, φ d ,is φ d = (2.58)? v g a ψ ˙ d where the desired turn rate, ψ ˙ d , is the output of the waypoint tracker controller, v a is the aircraft airspeed, and φ is the bank angle (positive for right wing down). The bank angle is regulated by generating a control signal to the ailerons, u aile , using PI control on the bank angle error, e φ (t)= φ ? (t)?φ d (t), where φ ? is the bank angle estimate obtained from ?ltered roll rate measurements with bias correction, t ? 1 ? ? T u aile (t)= K φp e φ (t)+ e φ (τ )dτ (2.59) φ 0 and K p and T φ are the controller gain and reset rate, respectively. This PI control strategy is appropriate for this application, as integral action is needed to ensure the tracking of the desired turn rate command from the waypoint tracker, while the separate roll damper loop handles upsets caused by turbulence. Note that another 60 (a) Tail view, with the roll angle, (b) Overhead view, with the air- φ, and lift forces, L, shown. speed vector, v a , and yaw an- gle, ψ, shown. Figure 2-17: Aircraft in a coordinated turn undergoing lateral acceleration, a. The ?gure assumes that the airspeed, v a , is constant. strategy using strictly yaw rate measurements is available, if the roll angle estimates cannot be used due to high noise levels. Using the HWIL simulator the turn rate controller can be tuned to provide good performance before ?ight testing. Figure 2-18 shows the response to a bank angle command as given by Eq. 2.58 for varying K p gains and constant reset rate, 1/T φ = 0.1. The plot shows the increased rise time and settling time for increasing K P gains, however at the cost of increased turn rate oscillations. The most appropriate gains found for the HWIL simulation model are shown with the thicker lines, indicating reasonable settling time of about 5 seconds, and only slight oscillations in the turn rate. Roll damper K The roll damper is used to counteract roll turbulence by direct feedback of the roll rate measurement, φ ˙ to the ailerons before it integrate to the e?ect the turn rate controller. A low pass ?lter (10 rad/s bandwidth) is applied to the roll rate measurement, and GPS is used to correct for the drift bias before feeding it through the roll damper, φ . Figure 2-19 shows the e?ect of applying the added derivative action to the turn rate loop, with increasing damper gain, K φd . The derivative action softens the roll response of the autopilot and increases the bank angle oscillations. Since the 60 ARF trainer aircraft already has a high degree of roll stability, very little, if any, damper 61 Response to Bank Angle Cmd [deg/s] [deg] Turn Rate Controller Tuning 30 25 20 15 10 5 0 Roll angle φ Response [deg] Yaw rate ψ Response [deg/s] . φ d [deg] K p =?0.1, K i =?.01 K p =?0.1, K i =?.01 K p =?0.3, K i =?.03 K p =?0.3, K i =?.03 K p =?0.5, K i =?.05 K p =?0.5, K i =?.05 K p =?0.7, K i =?.07 K p =?0.7, K i =?.07 ?5 0 2 4 6 8 10 12 14 16 Time [s] Figure 2-18: The performance variation of the turn rate loop with varia- tion of the proportional and integral gains, K φp , and K φi = K p T φ respectively. The bank angle estimate, φ ? , and a ?ltered yaw ? ˙ rate measurement, ψ, is shown in response to a step change in desired bank angle. The thick line provides the best tradeo? between roll response and undamped oscillations in turn rate. gain is needed to prevent turbulence upsets. The thicker line in Figure 2-19 shows the response for the slight damper settings that might be used in ?ight. Note that from the time constant of the graph the bank angle dynamics of the aircraft may be determined. Modeling the bank angle as a ?rst order lag of the desired command, the dynamics can be written φ(s) 1 = (2.60) φ d (s) T c s +1 with T c 1.25 seconds from Figure 2-19. This measurement will be later used in ≈ Section 4.2 to modify the MILP dynamics formulation to account for the lag in accel- eration that would be seen as a result of the step change. 62 Roll Damper Tuning 30 Response to Bank Angle Cmd [deg/s] [deg] 25 20 15 10 5 0 Roll angle φ Response [deg] Yaw rate ψ Response [deg/s] . K d = 0.0 [deg] K d = 0.0 [deg/s] φ d [deg] K d = 0.02 [deg] K d = 0.02 [deg/s] K d = 0.05 [deg] K d = 0.05 [deg/s] K d = 0.1 [deg] K d = 0.1 [deg/s] ?5 0 2 4 6 8 10 12 14 16 Time [s] Figure 2-19: Variation of Turn Rate performance with selection of the appro- priate roll damper gains, K φd with constant proportional and integral gains, K φp =0.3, K φi =0.01 respectively Rudder Mixing Aileron-Rudder mixing rules are employed to assist in making coordinated turns for the aircraft, reducing the sideslip motion while banking the aircraft. For this trainer aircraft this is achieved by setting the aileron to rudder mixing ratio, K ar =0.15, which helps to reduce the errors in GPS heading when making turns. The e?ective- ness of the mixing ratio could not be determined without additional measurements on the aircraft true heading, however by including an onboard magnetometer the sideslip motion can be measured through comparison of the GPS and true heading estimates. Added true heading measurements would facilitate the selection of the optimum rudder mixing ratio for a variety of ?ight conditions and airspeeds, thereby reducing the GPS errors due to sideslip motion in bank turns. 63 (a) Waypoint tracker geometry. (b) Tracker Response Depicted Figure 2-20: Lateral track control law for the Cloud Cap autopilot. The yaw damper loop as shown in Figure 2-16 is typically needed on aircraft with less lateral directional stability, however the trainer ARF 60 aircraft have a large enough dihedral angle, Γ = 5 ? , and vertical ?n to reject perturbations in yaw (note the relative stability of the derivatives in Table 2.4). Yaw damper augmentation of the plant dynamics is therefore not needed and disabled for the trainer ARF 60 aircraft. 2.3.2 Waypoint Tracker The Piccolo TM autopilot uses a robust nonlinear waypoint tracker, originally designed for implementation on the Aerosonde UAV [34], which is capable of tracking the vehicle through series of inertially ?xed waypoints using GPS position and velocity 64 measurements. The strategy employs PD control on the error signal e ψ = ψ ? ?ψ d , where ψ ? is the aircraft heading estimate and ψ d is the desired heading vector as shown in Figure 2-20(a). The heading estimate, ψ ? , is approximated using the GPS velocity vector, which is a valid assumption for low wind conditions and small aircraft bank angles (i.e., coordinated turns), ? v ay ? v ψ ? = arctan (2.61) ax The output of the waypoint tracker is a turn rate command, which is subsequently related to the bank angle of the aircraft as shown in Eq. 2.58 [35]. Convergence Parameter Scaling As shown in Figure 4-5, the desired heading vector, ψ d , is determined in the intrack reference frame by selecting a point L d meters ahead of the current along track po- sition. From Figures 2-20(a) and 2-20(b), the distance L d , (also referred to as the tracker convergence parameter ), determines how sensitive the waypoint tracker will be to cross-track errors and is selected to be approximately equal to the vehicle turn radius. For an aircraft in a coordinated turn at a speci?ed maximum bank angle, φ max =30 ? , with airspeed v a = 25 m/s the vehicle turn radius is given by 2 v L d ≈ρ min = g tan a φ max ≈ 120m (2.62) Simulations con?rm that selecting L d < 120 leads to instability in the waypoint tracker, as noted by highly oscillatory ?ight patterns of the vehicle in Figure 2-21. With the selection of the convergence parameter, the desired heading vector is given by ? y IT ? L ψ d = ?arctan (2.63) d where y IT is the cross track distance of the vehicle and the negative sign is required for correct operation. The convergence parameter, L d , is in fact not constant for all operating conditions, as there are certain situations which require adjustment to the closed loop dynamics. The control law described in Figure 2-20(a) works well when the GPS vector provides 65 Latitude [deg] 41.967 L p L = 100 p 41.966 41.965 41.964 = 125 Wind 41.963 5 m/s 41.962 41.961 41.96 ?70.892 ?70.89 ?70.888 ?70.886 ?70.884 Longitude [deg] Figure 2-21: Groundtrack of a ?xed counter-clockwise waypoint plan for dif- ferent selections of the intrack convergence parameter, L d . Dif- ferent selections determine the autopilot sensitivity to cross- track error, creating unstable oscillations for small values of L d . Headwinds tend to excite oscillations for the marginally stable system, as the turn rate command provides too much authority to the aileron actuators. an accurate estimate of the vehicle wind-relative heading. However, due to the e?ects of wind acting on the vehicle, the aircraft will experience changes in the e?ectiveness of the turn rate command with changes in the wind relative to the vehicle body axis. By scheduling L d with varying airspeed to increase the tracker performance in the case of winds, the sensitivity the tracking algorithm is mitigated over the variation in the no-wind case. The convergence parameter is also scaled according to the vehicle cross-track error, in order to provide higher performance for small cross-track errors, y <L d . In order IT to achieve this one of several proprietary methods is employed in the Cloud Cap 66 autopilot including a nonlinear scaling of L d with cross-track error [36]. This scale factor smoothly transitions the convergence parameter to small values for y <L d , IT providing a means to keep the waypoint tracker sensitive for all regions of operation. While detailed simulations of the system in Figure 2-16 can be performed using the aircraft models developed in earlier sections, the nonlinear scalings and ambi- guities in the autopilot internal control loops make accurate predictive simulations di?cult to achieve. The simulations lead to mismatches between the predicted and actual gains needed for HWIL or ?ight testing, although general trends and limits can be con?rmed. In practice these tests do not capture all the dynamics required to accurately simulate the closed loop for controller tuning purposes. As a result, the most e?ective way to select the control gains is to test on the HWIL simulator and subsequently verify on actual ?ight tests. Tracker Tuning The performance of the tracker loop depends primarily on the inner turn rate loop performance as well as the airspeed of the aircraft being ?own. Gains that are set at low airspeeds will typically not have as much authority as the same gains at higher airspeeds, while gains set at the high end of the airspeed range can lead to oscillations if left unchecked. For this reason, the HWIL simulator is exercised at several points of the ?ight envelope to ensure proper operation. K One suitable method for ?nding autopilot gains is to increase the proportional term, K p , until the critical point for sustained oscillations is reached. Once this critical point is found, the proportional term can be reduced by 20-30% in magnitude, depending on the airspeeds the aircraft will be ?ying. The derivative gain, K d , is then used to soften the autopilot response and further reduce oscillations as needed. The results of this process are shown in Figure 2-22, for an aircraft being ?own at 26 m/s, close to the upper limit of the speed range. With the proportional gain set at p = ?0.25, the derivative gain of K d = ?0.5 is shown to signi?cantly reduce the undamped response, with a settling time of 17 seconds and 10% overshoot. 67 Tracker Performance for derivative settings, K = 0.25, L d =120 m p 120 Crosstrack Error [m] 100 80 60 40 20 0 ?20 K d = 0.0 K d = ?0.1 K d = ?0.2 K d = ?0.3 K d = ?0.4 K d = ?0.5 0 5 10 15 20 25 30 Time [s] Figure 2-22: The response of the closed loop system is shown for varying derivative tracker gains, reducing the undamped oscillations of the system. Closed loop settling time from a step disturbance of 100 m in cross-track error is shown to be reduced to 17 seconds with 10% overshoot. Dynamic MILP Trajectory Conversion The output of the MILP planner is a series of ordered waypoint lists (?ight plans ) that are uploaded to the vehicle in real-time as each optimization completes. These waypoint lists represent the physical locations in (x, y) space that the MILP planner assumes the vehicle will be passing through at a given time, and they need to be converted into a form that is consistent with the closed loop response of the waypoint tracker. The waypoint lists are updated at a rate of roughly 1 plan every 4-5 seconds, depending on the wind conditions. As shown in Ref. [8], the process is repeated when the vehicle reaches the horizon point of the current plan. 68 Figure 2-23: Conversion of MILP waypoints to vehicle ?ight plans, utilizing the waypoint tracker control law to keep bank angle constant throughout a turn. With a small enough waypoint spacing, Γ ≈ 100m, the desired heading vector, H, will point as shown for the tra jectories with large changes in heading angle, Δψ. Figure 2-23 shows the tra jectory points {A,B,C,D} with spacing Γ that were designed by the MILP solver. To prevent the vehicle from experiencing undesirable transients when new plans are uploaded mid-turn, the modi?ed plan {A ? ,C,D} is uploaded to the vehicle. A ? is selected as the straight line projection of the line CB of distance Γ, as shown in Figure 2-23. The point B is not needed and is therefore removed from the plan upload to the vehicle, however B is used as the reference to the start of the next optimization (i.e., the Horizon Point ). The modi?ed plan, {A ? ,C,D}, closely resembles the situation in Figure 4-5, and allows the vehicle to maintain a constant bank angle as each point in the plan is reached. Note that in the current implementation the vehicle receives new plans approximately every 4 seconds (Γ = 100m), meaning that large cross-track errors will not have time to accumulate before new plans are uploaded to the vehicle. As a result, the convergence parameter is scaled using the nonlinear scaling law for small cross-track, and for most of the MILP tra jectory design the waypoint tracker will be operating in the scaled nonlinear region. This provides bene?ts to the performance of the system, but makes the closed loop response more di?cult to predict. 69 Figure 2-24: The longitudinal altitude, h, and airspeed, v a , control loops shown for the Cloud Cap autopilot. The total and static pres- sure, pitch rate and air density P t ,P 0 ,θ ˙ , and ρ, respectively, are measured and used to determine the throttle and elevator control inputs δ throt and δ e . There is also an optional PD loop from altitude to elevator not shown here. 2.3.3 Airspeed/Altitude Control Due to the inherent coupling between airspeed and altitude, the Cloud Cap autopilot utilizes several coupled air, pitch and altitude control loops to maintain the desired reference values. Altitude is not controlled directly, rather using an energy expres- sion to help damp phugoid oscillations. Due to the coupling between altitude and airspeed, the tuning of the various controller settings in Figure 2-24 is an involved process requiring several iterations to ensure satisfactory performance in both vari- ables. Later sections will utilize the reference speed commands to perform timing control of the vehicle, however the altitude will remain essentially ?xed for the most ?ight experiments. The true airspeed of the aircraft is estimated by measuring the total and static 70 pressures, P t and P 0 , of the incident air?ow ? 2P t P 0 v a = and ρ ≈ ρ 287T where ρ is the air density estimated using onboard measurements of the air tem- perature, T . The airspeed is regulated using PID control on the dynamic pressure error signal, e(t)= P t (t) ?P tD , where P tD =1/2ρv 2 is the desired total pressure a corresponding to the target airspeed of the aircraft. The tuning of the airspeed loop is shown in Figure 2-25, which shows the airspeed response to a step input of 6 m/s for varying cases of control selections, I, II, III, and IV. Case I shows the e?ect of a high proportional gain in the airspeed loop, as the elevator is excited too heavily producing undamped oscillations. The airspeed proportional gain is shown reduced by 25% in II, and setting the altitude derivative gain 25% higher in Step III reduces some of the oscillation due to throttle. Finally in step IV, the pitch damper term, K θp , is used to further decrease the oscillations, and the closed loop response is shown to have a settling time of approximately 25 seconds (note that the initial step at t = 5 seconds). In Figure 2-25(b), the altitude variation for the same series of tests is shown to lie within ±10 m of the nominal for these control settings. These control settings provide adequate performance for the tests required, however the entire space of options for control selection in this axis was not fully explored. It is therefore possible that higher performance solutions exist for the loops shown, and could be found with more detailed modeling of the closed loop dynamics. The closed loop response of the airspeed loop will be required in later chapters to develop an additional control loop on the relative timing of the mission plan being ?own. 2.4 Conclusions This chapter presented dynamics models for the trainer 60 aircraft which were used to validate the Cloud Cap HWIL simulation environment. By analyzing the response to speci?c modes of lateral and longitudinal aircraft dynamics, the reduced order 71 Airspeed Step Response 28 27 26 Altitude [m] Airspeed [m/s] 25 24 23 22 21 CMD I II 20 III IV 19 0 10 20 30 40 50 (a) Airspeed step response. Altitude Responses to Step in Airspeed Cmd 165 160 155 150 145 140 CMD I II III IV 0 10 20 30 40 50 Time [s] (b) The altitude response to a step change in desired airspeed, showing the coupling between the airspeed and altitude loops. Figure 2-25: Airspeed step response input and output of the airspeed con- troller. Note the timescale indicates the initial step at t =5 seconds. 72 models developed here for the ARF 60 aircraft were shown to have similar responses and characteristics to the models created from HWIL simulations. While a complete model of the ?ight dynamics could not be formed, it was shown that the HWIL simu- lator does represent certain key aspects of the dynamics in both the fast (short period, roll subsidence modes) and slow (phugoid and spiral modes) regions, indicating its suitability to represent the ARF 60 aircraft for autopilot controller tuning and high level HWIL simulations. The dutch roll mode was not modeled due to the di?cultyin collecting experimen- tal data isolating this response. The inherent roll-yaw axis couplings are notorious for being di?cult to experimentally verify, and reduced order analytical models require making gross approximations in the dynamics and aerodynamic derivatives. Had dif- ?culties arisen in making the transitions from HWIL simulation to ?ight testing more detailed models could have been considered. Section 2.3 described the Piccolo TM autopilot tunings set up to be utilized in the trainer ARF 60 aircraft. The closed loop system performance has been identi?ed for the lateral track waypoint controller, as well as the airspeed controller shown in Subsection 2.3.3, which will allow planning control loops to be closed in later sections with greater accuracy. The performance of the inner loops has been veri?ed on repeated occasions during actual ?ight tests. 73