《航空器的稳定与控制》（英文版）Distributed Coordination and Control

Distributed Coordination and control Experiments on a Multi-UAV Testbed Ellis T King elor of Engineerin The State University of buffalo, 2002 Submitted to the Department of Aeronautics and Astronautics in partial fulfillment of the requirements for the degree of Master of Science in Aeronautics and Astronautics t the MASSACHUSETTS INSTITUTE OF TECHNOLOGY September 2004 C Massachusetts Institute of Technology 2004 Authe Department of Aeronautics and astronautics August 20. 2004 Certified b Associate professor Thesis Supervisor Accepted by Jaime Peraire Professor of aeronautics and astronautics Graduate Stude

c Distributed Coordination and Control Experiments on a Multi­UAV Testbed by Ellis T. King Bachelor of Engineering The State University of Buffalo, 2002 Submitted to the Department of Aeronautics and Astronautics in partial fulfillment 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 Certified 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-loop(HWIL) simulation is only useful if it accurately portrays the vehicle dynamics and if the behaviors observed during fight tests can be repli cated on the ground. This chapter focuses on identifying some of the dynamical modes of fight for the 60 ARF Trainer aircraft, and verifying that the HWIL sim- ulations reflect 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 identification techniques on experimental fight 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-loop(HWIL) simulations and Section 2.2 details the procedures used to create models of the aircraft dynam ics from data collected during flight tests and hardware-in-the-loop simulations. I 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 fown. For simply

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 flight tests can be repli￾cated on the ground. This chapter focuses on identifying some of the dynamical modes of flight for the 60 ARF Trainer aircraft, and verifying that the HWIL sim￾ulations reflect 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 identification techniques on experimental flight 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 flight 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 flown. For simply 31

(a)The Clark YH airfoil geometry △△△△A会 deg (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 configured 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 fight conditions, providing the input parameters are configured 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 files is given in Ref. 26

(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] CL CD (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. configured 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 flight conditions, providing the input parameters are configured 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 files 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 Spa 1.707 b Wing Area 0.5200 S Chord length 0.305 Wing Incidence 1 de Wing dihedral deg Wing Sweep 0.0de Tail area 0.0879m2 Tail Spa 0.606 b t Fin area 0.0324m2 Fin Span 0.216m Fin Offset X 1.143 Fin offset Z 0.120 h Fin Sy deg A Fin volume 0.0189 0.130 Fuselage Length.270 lb Gross mass 5267kg Empty Mass 4.798 ne Roll iner 0.31 kg.m2I Pitch inertia 0.46kg·m2 Yaw inertia* 0.63kg·m22 Aircraft Inertia ex periment 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 rcraft, 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 ng

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 Offset X Tail Sweep Fin Area Fin Span Fin Offset X Fin Offset 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 m2 m deg deg deg m2 m m deg m2 m m m deg ­ m2 m kg kg kg m2 · kg m2 · kg m2 · b S c¯ Γ Λ Λ l b St t t t Sf bf lf Λ hf f V¯f I I I m l Sb 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

L R一 L sin(a)=R sin(g) ng Figure 2-2: Torsional pendulum experimental setup to determine roll axis inertia, Ixc, for the trainer aircraft. The period of oscillation of a roll angle perturbation, is measured to parameterize the aircraft inertia. The angle a 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 Iyy and Izz respectivel For rotational perturbations applied to the airframe, the product of interior angles and distances must be constant 2 where o is the aircraft roll angle perturbation and a is the small angle deviation of the supporting cables from the vertical position. The differential equation describin the motion of the torsional pendulum is governed by a torsional inertia term and the restoring moment due to tension forces Irao+ 2TRsina=0

Figure 2­2: Torsional pendulum experimental setup to determine roll axis inertia, Ixx, 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 Iyy and Izz 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 differential equation describing the motion of the torsional pendulum is governed by a torsional inertia term and the restoring moment due to tension forces ¨ Ixxφ + 2T R sin α = 0 (2.3) 34

Table 2.2: Experimentally determined aircraft inertias [kg-m21 Aircraft Roll Axis I Pitch Axis Yaw Axis N 0.46 0.65 0.61 0.63 Mean 0.30 0.42 0.63 Std Dev.I 0.029 0.011 0.021 Using the small angle approximation for a since L R, and substituting known values from Eqs. 2.1 and 2.2, Eq 2.3 reduces to Irao L which is characterized by the undamped natural frequency, w, and period of oscilla Tp R2 (2.5) 2 By finding averaged values for the period of oscillation, Tp, 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,, Ixz, Iy2). The experimental results are summarized for each of the axes and three different aircraft in the same configuration in Table 2.2, showing agreement between different vehicles used in the tests. The largest variation was found in the roll axis due to the difficulty 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-

� Table 2.2: Experimentally determined aircraft inertias [kg­m2] Aircraft Roll Axis Pitch Axis Yaw Axis No. Ixx Iyy Izz 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 ¨ mgR2 Ixxφ + φ = 0 (2.4) L which is characterized by the undamped natural frequency, ω, and period of oscilla￾tion, Tp ω = mgR2 IxxL (2.5) 2π Tp = ω (2.6) By finding averaged values for the period of oscillation, Tp, 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., Ixz, Iyz). The experimental results are summarized for each of the axes and three different aircraft in the same configuration in Table 2.2, showing agreement between different vehicles used in the tests. The largest variation was found in the roll axis due to the difficulty 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

wn.2 Actar h s2 +2*z'wnstwn'2 Actuator Actator Transfer Fcn Rate Limiter Figure 2-3: Actuator models used the Cloud Cap Hardware-in-the-loop sim- ulation n-the-loop simulator. As shown in Ref. 26, the actuator transfer function, Gact(s) is given by specifying the bandwidth limit, Bw 2丌B1w 8) =c=0.707 (29) where the damping ratio, S, is selected at the critical value to set the actuator band width equal to the natural frequency (wb Wn). The aileron, elevator, and rudder channels all respond with approximately the same characteristics(Bw= 10 Hz),but the throttle is modeled with less dynamic range(Bw= 2 Hz) as the engine RPM equires added time to ramp up to produce thrust. The input/output saturation and slew rate limits are determined as per manufacturer specifications(++60, 2 Hz respectively), and applied as shown in Figure 2-3 2.1. 3 Sensor noises The Cloud Cap hardware-in-the-loop 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 specified drift rates26. 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 pressure sensors, the CristaTM IMu angle-rate sensors. and the accelerometers are shown in Table 2.3

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, Gact(s), is given by specifying the bandwidth limit, BW Gact(s)= ω2 n s2 + 2ζωns + ω2 n (2.7) ωn = 2πBW (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 (BW = 10 Hz), but the throttle is modeled with less dynamic range (BW = 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 specifications (±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 specified 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 PiccoloTM pressure sensors, the CristaTM IMU angle­rate sensors, and the accelerometers are shown in Table 2.3. 36

Table 2.3: Crista IMU HWIL Sensor noise models PDynamic PStatic G unit P Pa] [deg/s] [m/s/s) Resolution unit 3.906 200|1.6E-46.0E-3 Min unit 300 0.0-5.20|-1000 4000110.0005.20 100.0 20.0 Butterworth Order 2 2 2 2 BW Cutoff Freq. Hzl 11.0 11020.0 20.0 Drift Rate unit 0.05 1.0 1.5E-42.0E-3 Max Drift value [unit] 15.0 100.00.01 0.20 Drift Update Rate [ s 5.0 5.0 1.0 1.0 2.1.4 Dryden Turbulence Stochastic turbulence disturbances are required for accurate HWIL simulations, as eal 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 effect 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 e airspeed reference velocity, Vo, and the turbulence intensity,o The selection of these parameters allows for the turbulence to be modeled accordin 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 r1+(u)2 Suq(w) L。1+3(÷)2 (211) +( n()=9L.1+3({) (1+(÷=u)2

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 Cutoff 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 effect 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, Vo, 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: Sug(ω)= 2σ2 uLu πVo · 1 1 + ( Lu Vo ω)2 (2.10) Svg(ω)= σ2 vLv πVo · 1 + 3( Lv Vo ω)2 � 1 + ( Lv Vo ω)2 �2 (2.11) Swg(ω)= σ2 wLw πVo · 1 + 3( Lw Vo ω)2 � 1 + ( Lw Vo ω)2 �2 (2.12) 37

0.8(m) Lm1+( S9() (214) 15 + where w is the spectral frequency of the turbulence and the aircraft wingspan, b, is used as a parameter in the angle rate filters to scale the effect of rotation on the main lifting surface. The subscripts u, U, w and p, g, r refer to the familiar body frame aircraft wind velocities and angle rates, respectively, thereby allowing independent classification of the turbulence in each axis. These spectral shaping functions are used to form shaping filters to give the body axis noise transfer functions 30 1 H 丌V1+-ms (216) H9(s) 1+√3{m (218) + H (1+(#) (220) Hrg(s) H (221) (#) 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 filters gg and ro 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 x l white noise input was used for each trial set. Larger scale lengths, L, increase the time constant of

S σ2 0.8 � πLw �1/3 pg(ω)= w 4b �2 (2.13) LwVo · 1+ � 4b ω πVo � � ω 2 Vo Sqg(ω)= �2 · Swg(ω) (2.14) 1+ � 4b ω πVo � ω �2 Vo Srg(ω)= �2 · Svg(ω) (2.15) 1+ � 3b ω πVo where ω is the spectral frequency of the turbulence and the aircraft wingspan, b, is used as a parameter in the angle rate filters to scale the effect 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 classification of the turbulence in each axis. These spectral shaping functions are used to form shaping filters to give the body axis noise transfer functions [30] � Lu 1 Hug(s)= σu 2 1+ Lu (2.16) πVo · s 1+ √3 Lv Vo s � Lv Vo � � 1+ Lv 2 Hvg(s)= σv (2.17) πVo · s Vo s � Lw 1+ √3 Lw Vo � � 1+ Lw 2 Hwg(s)= σw (2.18) πVo · s Vo � �1/6 0.8 � π 4b (2.19) 1/3 Hpg(s)= σw Vo Lw � 1+ � 4b s �� πVo s Hqg(s)= � Hwg(s) (2.20) Vo 4b � · 1+ s πVo s Hrg(s)= � Hvg(s) (2.21) Vo 3b � · 1+ s πVo 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 filters qg and rg. 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

Filters an velocities ngular Rates Hg→ Hos Figure 2-4: Block diagram for the 6 DOF Dryden Turbulence model. The velocity perturbations ug, Ug, wg are independent outputs of the filtered values of the turbulence scale lengths, L, intensity values g 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 filters. the turbulence seen for a given airspeed, while larger o values increase the deviation about zero. L and o typically vary with altitude in the lower atmosphere as ground effects become more prominent, but for HWIL simulations they are usually fixed The frequency response of the Dryden filters are shown for the same three cases in Figure 2-6. The filter cutoff frequency is determined by the ratio of the scale length to airspeed, and this effectively produces lower bandwidth filters for larger scale lengths The scale length parameter is chosen according to one of several specifications, 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

Figure 2­4: Block diagram for the 6 DOF Dryden Turbulence model. The velocity perturbations ug, vg, wg are independent outputs of the filtered 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 filters. 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 effects become more prominent, but for HWIL simulations they are usually fixed. The frequency response of the Dryden filters are shown for the same three cases in Figure 2­6. The filter cutoff frequency is determined by the ratio of the scale length to airspeed, and this effectively produces lower bandwidth filters for larger scale lengths. The scale length parameter is chosen according to one of several specifications, 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