
small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


Navigation and Control Technologies for Autonomous Micro Vehicles Sponsorship: Air Force Research Laboratory / Munitions Directorate. During this project we focused on four primary objectives which are listed below. 1. Magnetometer Integration. Integrate and flight test magnetometers with the current version of the autopilot (Spiral 1). 2. HeightAboveGround Sensor Integration. Integrate and flight test heightaboveground sensors with the current version of the autopilot (Spiral 1). 3. Automatic Gain Adjustment. Develop and flight test automatic gain adjustment algorithms that automatically tune the servo loops of the autopilot. 4. Automatic Trim Seeking. Develop and flight test automatic trim seeking algorithms to recursively estimate the trim values of the UAV. Final Report Award No. F08630 03 1 0017 Navigation and Control Technologies for Autonomous Micro Vehicles August 25, 2005 Principal Investigator/ Technical Point of Contact Randal W. Beard, Associate Professor Department of Electrical and Computer Engineering Brigham Young University Provo, Utah 84604 USA voice: ( 801) 422 8392 fax: ( 801) 422 0201 beard@ ee. byu. edu Co Principal Investigator Timothy W. McLain, Associate Professor Department of Mechanical Engineering Brigham Young University Provo, Utah 84604 USA mclain@ byu. edu AFRL/ MN Point of Contact Carrie Fowler AFRL/ MNAV 101 West Eglin Blvd, Suite 346B Eglin AFB, FL 32542 ( 850) 882 8876 x3383 carrie. fowler@ eglin. af. mil CONTENTS Contents Cover Page i Table of Contents iii 1 Objectives 1 2 Magnetometer Integration 1 2.1 Magnetometer Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2.2 Motor Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3.1 Ground Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3.2 Airborne Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4 Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4.1 2 Axis Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4.2 3 Axis Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5 2 Axis Tilt Compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.6 Magnetometer Wind Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Height Above Ground Sensor 9 3.1 Landing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Orbit Routine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.1 Heading Field Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.2 Correcting Steady State Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2.3 Glide Slope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 Barometric Altimeter Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4 Ultrasonic Range Finder Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Optic Flow Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Automatic Trim and Gain Tuning 23 4.1 UAV Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 Scheme A: nonlinear MRAC using backstepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3 Scheme B: nonlinear MRAC without backstepping . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Scheme C: nonlinear MRAC grouping terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.4.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.4.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.5 Scheme D: nonlinear MRAC simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 ii CONTENTS 4.6 Flight Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.6.1 Platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.6.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5 Students Funded Under Project 73 A Appendix: PNI Micromag Specifications 74 iii 1 Objectives During this project we focused on four primary objectives which are listed below. 1. Magnetometer Integration. Integrate and flight test magnetometers with the current version of the autopilot ( Spiral 1). 2. Height Above Ground Sensor Integration. Integrate and flight test height above ground sensors with the current version of the autopilot ( Spiral 1). 3. Automatic Gain Adjustment. Develop and flight test automatic gain adjustment algorithms that automatically tune the servo loops of the autopilot. 4. Automatic Trim Seeking. Develop and flight test automatic trim seeking algorithms to recursively estimate the trim values of the UAV. The last two objectives were handled within one uniform framework and will therefore be described together. 2 Magnetometer Integration Adding a magnetic compassing device to the sensing suite of the Kestrel Autopilot is very desirable for several reasons. The BYU autopilot uses GPS for an estimate of ground track. With the use of a magnetometer, true heading ( ) can also be estimated. The use of true heading instead of ground track could significantly improve the accuracy of our gimbaled camera tracking system. Without true heading the gimbaled tracker will be in error by at least the value of the crab angle and the crab angle can be significant in high winds. In addition, magnetic heading allows for the computation of the wind vector which can easily be estimated in real time by using magnetometer heading, GPS velocity, airspeed, and ground track heading. 2.1 Magnetometer Selection When first selecting a magnetometer to interface with the autopilot several models were considered including the Honeywell HMC 1022 2 axis analog magnetometer and the PNI Micromag 2 axis digital magnetometer. Both sensors have a footprint of less than 0.6 x 0.5 inches. However the PNI Micromag magnetometer was selected because it incorporates built in temperature compensation and uses a SPI bus protocol that was straight forward to implement on the autopilot. The PNI Micromag also has selectable resolution and sample rate and the technical support at PNI is excellent. PNI developed a 3 axis version of the Micromag within several months which is currently being used on the autopilot. A table of Specifications and pin outs is included in Appendix A. Figure 1: PNI Micromag two axis digital magnetometer 1 2.2 Motor Noise 2.2 Motor Noise It was theorized that the magnetometer heading would be adversely influenced by induced magnetic fields from large brushless motors that are used for propulsion onboard Unmanned Air Vehicles ( UAVs). Fortunately, it was found that motor noise could be mitigated by twisting and shielding the high current carrying wiring. The alternating nature of brushless motors causes high frequency low amplitude noise and can also introduce a bias on the magnetometer readings. The high frequency component of the noise can be filtered out digitally. The biasing component of the motor noise is correlated with motor current. Work is currently being done to remove the effects of the electric motor on the magnetometer readings. The magnetometer is mounted on the autopilot which is located about 4 to 6 inches from the motor. If possible this distance should be increased to further reduce the effects of noise. All current carrying wiring should also be placed as far from the magnetometer as possible in addition to being twisted and/ or shielded. 2.3 Calibration In order to use magnetometers effectively a good calibration must be employed. This is easier said than done. First the magnetometer must be placed in the UAV in its completed form. Every time that something metallic is added to the UAV or every time the magnetometer is relocated on the UAV a new calibration should be performed. The magnetometer should be mounted as far as feasibly possible from current carrying wires, ferrous metals, and magnetic components. Then a magnetically clean environment should be located for calibration. A magnetically clean environment is one that is distant from steel and iron structures, power lines, and magnetic items. As a rule of thumb the calibration should be performed at least 5 times the diameter away from the nearest ferrous and magnetic objects. For example, if a 10 foot automobile is nearby, move at least 50 feet away from the automobile before calibrating the magnetometers. Beware that underground power lines and nearby power transformers could significantly affect the calibration. 2.3.1 Ground Calibration If it is intended to calibrate on the ground, the person performing the calibration should be sure to remove any metal objects and powered devices, such as keys and cell phones, from the area. Next, the ground station operator should set Calibrate Magnetometers to a value of 1. This starts logging data for the X and Y axes and stores the maximum and minimum values ( X m a x , X m i n , Y m a x , Y m i n ) for each axis. Then the person calibrating should hold the plane as level as possible while slowly turning in a circle in the horizontal plane through two complete rotations. Next the ground station operator should set Calibrate Magnetometers to a value of 2. This stops logging magnetometer data and calculates the scale factors ( X scale, Y scale) and offset values ( X offset, Y offset) for each axis. To ensure good calibration it is recommended that the X and Y magnetometer data be plotted on the x y plane and displayed for the operator. The X and Y magnetometer data should appear to be an ellipsoid. Applying the scale factors and offsets should circularize the ellipsoid and shift its center toward zero on the plot ( see Figure 2). For a 3 axis magnetometer the Z  axis should be calibrated separately. To do this the ground station operator should set Calibrate Magnetometers to 3. This starts searching for maximum and minimum Z axis values ( Z m a x , Z m i n ). The airplane must then be rotated about either the X or Y  axis for a full rotation. It is also valid to only turn the aircraft completely upside down, gathering the data points in both orientations. Once the rotation is completed Calibrate Magnetometers should be set to 4 to calculate the Z scale factor and offset values ( Z scale and Z offset respectively). The scale factors and offsets are calculated using the following equations: X offset = ( X m a x + X m i n ) = 2 ; Y offset = ( Y m a x + Y m i n ) = 2 ; Z offset = ( Z m a x + Z m i n ) = 2 ; ( 1) X scale = ( X m a x X m i n ) = 2 ; Y scale = ( Y m a x Y m i n ) = 2 ; Z scale = ( Z m a x Z m i n ) = 2 : ( 2) Once the calibration is complete be sure to save the calibration parameters in the autopilot flash memory and to a file. 2 2.4 Heading Estimation  2000  1500  1000  500 0 500 1000 1500 2000  2000  1500  1000  500 0 500 1000 1500 2000 X axis counts Y axis counts raw data calibrated data Figure 2: Y axis data plotted against X axis data for a typical two axis calibration. The uncalibrated data is in red while the data with the applied scaling and offset is blue. 2.3.2 Airborne Calibration The most effect method for calibrating the magnetometer was found to be an airborne calibration. An airborne calibration places the UAV in an environment that is far from magnetic disturbances on the ground. To calibrate the magnetometers in the air the first step is to set up an orbit point with a large orbit radius ( between 200 and 300 meters). Next, get the UAV up in the air and flying in its commanded orbit. Then start the magnetometer X Y calibration by setting Calibrate Magnetometers to 1. This starts logging data for the X and Y axes and stores the maximum and minimum values ( X m a x , X m i n , Y m a x , Y m i n ) for each axis. Once the plane has completed at least one complete orbit, set Calibrate Magnetometer to 2 which stops logging magnetometer data and calculates the scale factors ( X scale, Y scale) and offset values ( X offset, Y offset) for each axis. To calibrate the Z  axis, fly on an east west flight path. Then set Calibrate Magnetometers to 3. Next, invert the UAV momentarily before returning to normal flight. Barrel rolls could also be performed to achieve the same goal. Finally, set Calibrate Magnetometer to 4. This completes the Z  axis calibration and calculates the Z scale factor and offset values ( Z scale and Z offset respectively). As an alternative method for calibrating the Z  axis, set up a north south flight path. Then set Calibrate Magnetometers to 3. Next perform a complete loop or split S maneuver. Then set Calibrate Magnetometers 4. In either case the goal of the calibration routine is to capture the maximum and minimum values of all three axes. Be sure to save the calibration parameters in the autopilot flash memory and to a file. Once the rotation is completed the scale factors and offsets are calculated using Equations ( 1). 2.4 Heading Estimation 2.4.1 2 Axis Heading Estimation To calculate heading, raw magnetometer readings ( B x and B y ) are adjusted using the scale factors and offsets calculated during calibration using the following equations: B 0 y = B y Y 0 Y g B 0 x = B x X 0 X g : 3 2.4 Heading Estimation Figure 3: 3 axis heading ( red) compared to 2 axis heading ( green) along with corresponding roll ( ) and pitch ( ) angles taken during a clockwise circular orbit. To estimate heading the following formula is used: = atan2 ( B 0 y ; B 0 x ) + declination The heading estimate on a 2 axis magnetometer is only valid when the magnetometer is oriented level. A 2 axis magnetometer will be in error by as much as 20 degrees with roll and pitch angles of up to 15 degrees ( see Figure 3). 2.4.2 3 Axis Heading Estimation To calculate heading, raw magnetometer readings ( B x , B y , and B z ) are adjusted using the scale factors and offsets calculated during calibration, where B 0 x , B 0 y , and B 0 z are the now calibrated x , y , and z components of the magnetic field, i. e., B 0 x = B x X 0 X g B 0 y = B y Y 0 Y g B 0 z = B z Z 0 Z g : Then the adjusted magnetometer readings are un rolled and un pitched using standard rotation matrices as follows: 0 @ X h Y h Z h 1 A = 0 @ c o s 0 s i n 0 1 0 s i n 0 c o s 1 A 0 @ 1 0 0 0 c o s s i n 0 s i n c o s 1 A 0 @ B 0 x B 0 y B 0 z 1 A : ( 3) The un pitched and un rolled horizontal components of the magnetic field then become: X h = B 0 x c o s B 0 y s i n s i n B 0 z c o s s i n Y h = B 0 y c o s B 0 z s i n : From these equations heading is computed as: = atan2 ( Y h ; X h ) + declination : 4 2.5 2 Axis Tilt Compensation Figure 4: 3 axis estimation of heading as the UAV flew in a clockwise orbital flight path. Corresponding roll and pitch angles are given in the lower plot This method been used successfully to estimated the actual heading of the UAV even while maneuvering at pitch and roll angles of over 35 degrees ( see figure 5). 2.5 2 Axis Tilt Compensation A 2 axis magnetometer can be pitch and roll compensated to achieve results comparable to those of a 3 axis magnetometer by estimated the z component of the magnetometer. This is done by using additional information that is known about the Earth s magnetic field such as inclination and magnitude of the Earth s field. Inclination is the angle between the Earth s surface and the total magnetic field vector. This angle, often referred to as the dip angle, mainly varies depending on latitude. Inclination is not to be confused with declination which is the angular difference between true north and magnetic north. The magnitude of the Earth s magnetic field is relatively constant, however, aperiodic variation of earth s magnetic field over time and place does occur. The Z  axis component can be estimated using the following expression: B 0 z = s i n B Total + s i n B 0 x c o s s i n B 0 y c o s c o s ; where is the inclination and B Total is is the magnitude of the Earth s magnetic field. Both of these values are only valid in the region where the magnetometer is to be used. In the case of the PNI magnetometer which outputs readings in terms of ADC counts rather than in Tesla or Gauss, The value for B Total must also be input in terms of ADC counts. After estimating the z component, the adjusted magnetometer readings are un rolled and un pitched using Equation ( 3) with B 0 z replaced by B 0 z . This method has been used successfully to estimate the actual heading of the UAV even while maneuvering at pitch and roll angles of over 35 degrees. In simulation this method of estimating the Z  axis is equally as good as using a 3 axis magnetometer ( see Figure 5). 2.6 MagnetometerWind Estimation Knowledge of true heading and ground track heading as well as airspeed and ground speed allows for the computation of the wind vector which can easily be estimated in real time. Using what is called the wind triangle, the wind speed and wind direction can both be found using simple trigonometry ( see figure 6). Using the wind triangle to subtract the air velocity vector from the ground velocity vector yields the wind vector. This method has accurately measured the wind speed and wind direction in flight on a UAV ( see figure 7). 5 2.6 Magnetometer Wind Estimation 0 10 20 30 40 50 60 70 80 90  200  100 0 100 200 Comparison Between Three magnetic heading estimation methods degrees 0 10 20 30 40 50 60 70 80 90  5 0 5 10 15 Time ( sec) degrees psi 3D Mag 2D Mag 2D Mag Est 3D Mag f q Constant roll of 10 Degrees Figure 5: A simulation comparison of all three methods discussed for estimating psi including using a 3 axis mag ( 3D Mag), 2 axis mag ( 2D Mag), and 2 axis tilt compensation by estimating the Z axis ( 2D Mag Est 3D Mag). Corresponding roll and pitch angles are given in the lower plot. Wind speed / Direction Airspeed / Heading ( y ) Ground speed / Ground track ( c ) Figure 6: Wind triangle showing that the addition of the air velocity vector and wind vectors should equal the ground velocity vector. Additional magnetometer data is shown in Figure 8 and 9. 6 2.6 Magnetometer Wind Estimation 0 20 40 60 80 100 120 300 320 340 Wind Estimation ( Actual Wind: from North West at 1 to 2 m/ s) Wind Direction ( degrees) 0 20 40 60 80 100 120 0 2 4 6 Wind Speed ( m/ s) 0 20 40 60 80 100 120 240 260 280 300 320 degrees sample number Estimated Reported Estimated Reported GPS Course Mag Heading Figure 7: Top: Estimated wind direction. Middle: Estimated wind speed. Bottom: GPS Course and Magnetometer heading. The actual wind speed and direction as reported on day of flight were 2 m/ s from the north west ( 315 degrees). 0 200 400 600 800 1000 1200 1400 1600 1800  4  2 0 2 4 x 10 4 x y z 0 200 400 600 800 1000 1200 1400 1600 1800  4  2 0 2 4 x 10 4 total field magnitude 0 200 400 600 800 1000 1200 1400 1600 1800  50 0 50 sample num degrees phi theta Figure 8: Top: Magnetometer flight data of X , Y , and Z axes. Middle: Total magnetic field vector calculated from X , Y and Z data. It is relatively constant which indicates that the magnetometer did not enter the presence of any significant magnetic disturbances. Bottom: Roll ( ) and Pitch ( ) data corresponding to the collected magnetometer data in flight. 7 2.6 Magnetometer Wind Estimation Figure 9: Heading estimates in a circular orbit using both 2 and 3 axis magnetometers. 8 3 Height Above Ground Sensor The usefulness of mini and micro unmanned air vehicles ( UAVs) depends largely on the ability to hand launch and recover the UAV in uncontrolled environments. Efforts to implement effective autonomous landing of mini UAVs in uncontrolled environments have been most hindered by lack of small, lightweight, effective sensors to determine height above ground ( HAG). Because of the lack of a suitable sensor for determining HAG, small UAVs have traditionally relied solely on the barometric altimeters with which they are already equipped for HAG estimation. This approach works well for controlled environments in which the altitude of the landing environment relative to the altitude of the takeoff environment is known, and the terrain is relatively smooth. However, in uncontrolled environments where this offset is not known, or where the terrain is not acceptably smooth, other methods must be utilized to effectively determine HAG. Two methods with acceptably small sensor payloads have been proposed as effective means of HAG estimation. These methods are: 1. Estimation of HAG via Ultrasonic Ranging 2. Estimation of HAG via Optic Flow It is proposed that these sensors can be combined with information from the barometric altimeter to reliably and accurately determine HAG. This accurate estimate of HAG could then be combined with robust control algorithms to effectively and safely land small UAVs in uncontrolled environments. This report presents the experimental results of efforts to estimate HAG and land small fixed wing UAVs in uncontrolled environments using each of the sensors outlined above. Because each of the sensors tested will couple with the barometric altimeter in determining HAG, the experimental results will be presented as follows: 1. Overview of the landing algorithm, 2. Landing results using only the barometric altimeter, 3. Landing results using an ultrasonic sensor combined with the barometric altimeter, 4. Landing results using an optic flow sensor combined with the barometric altimeter, 5. Comparison of results using each of the sensors. Each of these sections will address the specific sensor tested and the algorithms used to determine HAG from the data provided by the sensor. 3.1 Landing Algorithm 3.1.1 Overview The landing algorithm developed and tested for this project consists of two user defined waypoints: the approach point and the landing point ( Figure 10). LANDING POINT APPROACH POINT Figure 10: Landing Algorithm The approach point serves as a point about which the airplane may orbit until it reaches a specified altitude, at which point it will break out of its orbit and follow a glide slope to the landing point. The approach point is fully defined by the following parameters: 9 3.2 Orbit Routine 1. Airspeed The desired airspeed which will be maintained throughout the circle down portion of the landing. 2. Radius The radius of the orbit which will be maintained throughout the circle down portion of the landing 3. Start Altitude The altitude to which the plane will attempt to climb or descend while en route to the approach point. The true start altitude for the descent will be the altitude to which the plane has successfully climbed or descended by the time it is within 1.5 orbit radii of the landing point. 4. End Altitude  The altitude at which the airplane will break out of its orbit and begin its glide slope in to the landing point 5. Relative East Coordinate  The distance east of home ( in meters) around which the plane will orbit while descending from its start altitude to its end altitude. 6. Relative North Coordinate  The distance north of home ( in meters) around which the plane will orbit while descending from it start altitude to its end altitude. 7. Descent Rate  The rate ( in meters/ second) at which the plane will descend while orbiting the approach point. The landing point is fully defined by the following parameters: 1. Airspeed  The airspeed which will be maintained throughout the glide slope from approach point to landing point. 2. Flare Height  The height above ground at which the airplane will cut throttle and attempt to hold zero roll while maintaining its glide slope by commanding pitch to control around desired altitude. 3. Relative East Coordinate  The distance east of home ( in meters) at which the plane will touch down. 4. Relative North Coordinate  The distance north of home ( in meters) at which the plane will touch down. Altitude control throughout the landing routine is accomplished through activation of various PID loops. For the portion of the landing sequence in which the plane is en route from its last waypoint to the approach point altitude is maintained as follows: 1. Actual altitude is within a preset distance ( usually around 15 m) of desired altitude  Pitch from Altitude, Throttle from Airspeed. 2. Actual altitude is greater than the same preset distance from desired altitude and is less than desired altitude  Pitch from Airspeed, Full Throttle. 3. Actual altitude is greater than the same preset distance from desired altitude and is greater than desired altitude  Pitch from Airspeed, Zero Throttle. For the portion of the landing sequence in which the plane is circling down to its end altitude and then following its glide slope to the landing point altitude is maintained using the pitch from altitude and throttle from airspeed loops. Throughout the entire landing sequence, until the flare at which desired roll is held at zero, heading is controlled using the roll from heading PID loop. Roll rate and pitch rate are also controlled to provide damping. Desired headings are commanded based on the generation of a heading field which forms a limit cycle around the desired orbit radius. A full description of this algorithm is provided in the following section. 3.2 Orbit Routine 3.2.1 Heading Field Generation The ability to implement the landing routine described in the previous section required the development and implementation of a control algorithm which would allow the airplane to fly accurate orbits, at minimal orbit radii, even in high wind conditions. The orbit algorithm implemented generates a desired GPS heading based on the position of the airplane such that the corresponding heading field forms a limit cycle around the desired orbit radius ( Figure 11). Such a heading field is generated using Equations ( 4). 10 3.2 Orbit Routine Figure 11: Typical Heading Field for a 50 Meter Radius Orbit ( Desired Orbit Radius Indicated by Solid Circle). d = 8 > < > : center t a n 1 r c d ; d 2 r c center r ; r c d 2 r c center + + + r ; d < r c ( 4) In the above equations d is the distance from the center of the orbit, d is the desired course, center is the compass heading to the center of the orbit, r c is the desired orbit radius, is the included angle between d and center for d = 2 r c , is the complimentary angle to , and is a gain corresponding to the field strength of the orbit. Increasing makes the field more reactive around the desired orbit radius and decreasing makes the field less reactive as shown in Figure 12. r is a parameterization of the distance from the center of the orbit based on r c . The formulation of r is such that 0 r 1 and r = 1 when d = r c . The value of r is calculated according to Equations ( 5). r = ( 1 d r c r c ; r c d 2 r c d r c ; d < r c ( 5) Figure 12: High Gain Heading Field ( Left) and Low Gain Heading Field ( Right) for a 50 Meter Radius Orbit ( Desired Orbit Radius Indicated by Solid Circle). The desired heading from the heading field described by equations ( 4) is used to generate a roll command for the airplane according to d = k ( d a ) ( 6) 11 3.2 Orbit Routine 3.2.2 Correcting Steady State Error Because the heading control represents a type one system and flying an orbit based on desired heading is essentially tracking a ramp in heading, the actual course of the airplane tracks the desired course with a steady state error. This steady state error in course represents a lag between desired and actual course which translates to actual orbits that are larger than the commanded orbit. The magnitude of the steady state error in orbit radius is a function of the desired radius of the orbit, the groundspeed of the airplane, and the heading field gain . Two different approaches were used to overcome this effect. The first approach involved adapting the radius command to remove the steady state error. The second approach, which was more effective in high winds, involved calculating a nominal roll angle based on the desired heading rate and the coordinated turn equation. Adaptive Radius Commands Because the actual orbit radius of the airplane is always exterior to the commanded radius, for every dynamically feasible orbit there exists some smaller radius command that will result in the actual desired orbit. The adaptive radius command progressively increases or decreases an offset between the desired orbit radius entered by the operator and the desired orbit radius passed to the orbit algorithm. This offset is continuously incremented or decremented until the plane is orbiting at the radius requested by the operator. Orbit results with and without the adaptive radius command are shown below in Figure 13. Telemetry from actual flight tests of the adaptive radius method are shown in Figure 14. Figure 13: Circular Orbits ( From Simulation) with and without an adaptive radius command. Nominal Roll Angle Based on Desired Heading Rate The adaptive radius command approach works extremely well for situations where wind speed does not exceed 50% of airspeed. However, high winds cause the groundspeed of the airplane to vary significantly based on the course of the airplane relative to the wind direction. This variation in groundspeed causes a corresponding variation in the radius offset necessary to maintain the desired orbit. A better approach for high wind conditions was developed which involved calculating a nominal roll angle that resulted in the heading rate necessary to fly the desired orbit. The heading field was then used to give roll commands about this nominal roll value. Calculating the nominal roll angle first required taking the derivative of the heading field with respect to time. This resulted in the formation of a heading rate field over the same space as the heading field. The equations for the heading rate field are given by _ d = V r c r 2 c + 2 d 2 d ( r 2 c + d 2 ) 3 2 ( 7) = V d s i n ( + r ) + r ( 1 ) V c o s ( + r ) r c ( 8) = V d ( s i n ( + r ) r c o s ( + r ) ) ( 9) An example of a heading rate field for a given heading field strength and constant groundspeed is shown in Figure 15. Equations ( 4) combined with Equations ( 7) ( 9) allow the calculation of both a desired heading and a desired heading rate for any location relative to the center of the orbit. The desired heading rate can then be converted to a 12 3.2 Orbit Routine Figure 14: Actual Flight Results for 50, 70, 100, and 150 m Radii Orbits Flown in 2 3 m/ s Wind . nominal roll angle using the standard coordinated turn equation nom = _ g V g ( 10) and the control law described in Equation ( 6) can be adapted as nom = _ g V g ( 11) d = nom + k ( d a ) ( 12) The control law shown in Equation ( 12) is more effective at tracking orbits in high winds, but somewhat less steady at tracking orbits in low wind conditions than is the adaptive radius command method. A comparison of the two methods flying a 100 meter radius orbit in simulation with a constant wind speed of 8 m/ s is shown in Figure 16. 3.2.3 Glide Slope The glide slope implemented in this landing routine was a simple linear glide slope calculated by decreasing desired altitude as a function of the distance from which the glide slope began. This formulation allows the desired altitude to be driven negative for distances from the initiation of the glide slope which are farther than the desired landing point. This makes the glide slope robust with respect to negative offsets between measured HAG and actual HAG. Expected overshoot or undershoot of the landing point can approximated as Expected Over / Undershoot 1 t a n a + b ; ( 13) where is the glide slope, a is the altitude control error ( difference between desired and actual HAG), and b is the measurement error. The relationship shown in Equation ( 13) indicates that larger glide slopes will tend to reduce the overall overshoot or undershoot error for the landing provided the plane is capable of matching the commanded glide 13 3.3 Barometric Altimeter Based Landing Figure 15: Typical Heading Rate Field. slope. However, larger glide slopes also cause increased impact velocities at touch down. Testing showed a glide slope of anywhere from 7 to 12 degrees to provide a good balance between accuracy and impact at touch down for the MAGICC Lab airplanes. 3.3 Barometric Altimeter Based Landing 3.3.1 Theory The Kestrel autopilot comes equipped with a barometric pressure sensor as part of its standard sensor suite. Before a given flight, the barometric pressure sensor is zeroed such that it references the barometric pressure at the takeoff location as the pressure corresponding to zero altitude. Pressure decreases with altitude according to P = P ref Z h h ref g h d h : ( 14) If the density of air ( ) is assumed to be sufficiently constant over the range of altitudes being measured, then altitude above the takeoff point can be determined by h = P ref P actual g = ( P ref P actual ) ; ( 15) where the constant of proportionality ( ) can be determined by taking readings of the pressure sensor at a number of different known reference altitudes and fitting a straight line to the data set. Once the constant of proportionality has been determined, readings of the pressure sensor from the analog to digital converter can be converted directly into relative altitudes above the takeoff point. This sort of calibration allows the barometric pressure sensor to provide valid and accurate estimates of HAG where the terrain is sufficiently flat. Where this is not true, barometric pressure can still provide a rough estimate of HAG until other sensors are able to collect valid readings of true HAG. After another sensor has reported a valid reading of true HAG variations in barometric pressure from the pressure at which true HAG was determined can be used to further estimate the true HAG. 3.3.2 Application Temperature Compensation The pressures reported by the pressure sensor on the Kestrel autopilot proved to be somewhat sensitive to the temperature of the sensor. Particularly this would become a problem when the autopilot was allowed to heat up as it sat static during calibration and programming, and then cooled as airflow over the airfoil dissipated heat during flight. This calibration at high temperature and subsequent cooling would cause the autopilot to consistently read approximately 20 meters high. The increase in pressure readings with decreased sensor temperature was empirically determined to be linear over the range of operating temperatures. The constant of proportionality was likewise empirically determined and implemented as P actual = P measured + ( T measured T ref ) ( 16) 14 3.3 Barometric Altimeter Based Landing Figure 16: A Comparison of Nominal Roll Angle and Adaptive Radius Based Methods ( Simulation). This temperature compensation makes the pressure readings from the autopilot robust to the variations in temperature that are encountered in standard autopilot operation. Compensation for Airflow Over the Airfoil Another source of error in the readings from the autopilot s pressure sensor stems from the fact that the MACICC Lab airplanes are not normally equipped with static pressure tubes. While in flight a considerable offset in the absolute pressure measured by the pressure sensor is introduced by airflow over the airfoil in which the autopilot is housed. The increased speed of the airflow over the airfoil results in a Bernoulli effect that causes the pressure measured along the airfoil to be less than the true absolute pressure. This low pressure region in which the actual pressure measurements are taken causes the altitude measurement computed by the autopilot to consistently read high. The relationship between the velocity of the airflow over the airfoil and the induced offset in the pressure reading was determined empirically through wind tunnel testing to be a second order function of airspeed. This functional relationship can be combined with the temperature correction term in Equation ( 16) to calculate pressure as P actual = P measured + C 1 v 2 diffpres + C 2 v diffpres + ( T measured T ref ) ( 17) where v diffpres is the airspeed past the airplane as measured by the airplane s pitot probe attached to a differential pressure port. Once the constants ( C 1 , C 2 , ) have been empirically determined, the offset can be computed each time through the loop, and the barometric altitude reading can be adjusted accordingly. Testing has shown the constants determined for one airplane to be valid for other airplanes of a similar type regardless of slight variations in the airfoil and in autopilot location. This allows the constants for a fleet of similar airplanes to be determined through a rigorous wind tunnel calibration of a single airplane. Further testing is required to determine how robust these constants are with respect to major changes in airframe or autopilot location. This method of compensating for airflow over the airfoil allows a large number of similar airplanes to be flown, using the exact same constants, without each being equipped with static tubes which can add unnecessary weight and complication. 3.3.3 Results Using the calibration described above and testing over relatively short flights, the MAGICC Lab airplanes consistently read within 1  1.5 meters of the actual altitude relative to home. This allows relatively accurate landings where the altitude of the landing point can be safely assumed to be the same as the altitude of home. This works especially well when trying to land the airplane at or near its home location which is most often the case. 15 3.4 Ultrasonic Range Finder Based Landing If the actual altitude of the airplane matched the desired altitude perfectly, then overshoot or undershoot would be solely a function of error in HAG estimation. For glide slopes of 7 to 12 degrees and a measurement error of 1  1.5 meters, a corresponding overshoot or undershoot error of around 4  12 meters can be expected ( see Equation ( 13)). Actual results for 26 consecutive landing runs are shown below in Figure 17. The first 22 autoland runs were each from the main approach point. These were followed by one run apiece from each of the four alternate approach points. Figure 17: Results of Twenty six Consecutive Autolands from Five Different Approach Points. An average overshoot or undershoot error of 4  12 meters was estimated. Actual results from these 26 landing show an average over/ undershoot error of 7.6 meters with a standard deviation of 5.4 meters and a median value of 6 meters. 3.4 Ultrasonic Range Finder Based Landing 3.4.1 Theory Ultrasonic sonar transducers determine distance from an object by measuring the time of flight of an ultrasonic signal which is emitted by the transducer and then reflected by the object. The speed of sound in an ideal gas is given by c = p R T k ; ( 18) where is the adiabatic constant ( 1.402 for air), R is the universal gas constant ( 297.05 kJ/( kg K) for air), and T is the absolute temperature in Kelvin. Linearizing Equation ( 18) around 0 degrees C gives a convenient approximation for the speed of sound as c air = 3 3 1 : 5 + 0 : 6 T c : ( 19) 3.4.2 Application Sensor The MUST01 Ultrasonic Sonar Transducer manufactured by Mekatronix was selected and used as the ultrasonic ranging module and transducer for these tests. This sensor provided both the ranging board and transducer in one integrated package as well as allowing a sensing range of up to 12.2 meters with a resolution of 4.8 cm. The entire unit has an overall diameter of 4.3 centimeters and a total mass of 17 grams. The sensor emits a 40  50 kHz sound wave and has a nominal beam angle of 15 degrees. The MUST01 unit emits a signal 10 times a second and returns an analog voltage proportional to each signal s time of flight. This value is alpha filtered on the ranging board according to V out = 0 : 9 V new + 0 : 1 V t e x t l a s t : ( 20) 16 3.4 Ultrasonic Range Finder Based Landing Figure 18: MUST01 Ultrasonic Sonar Transducer. The voltages from the ultrasonic sensor are read into the autopilot via the autopilot s analog to digital converter. The ultrasonic sensor was calibrated by holding the sensor at a number of fixed distances from a surface with good acoustical properties and then recording the analog to digital converter values corresponding to that distance. The A/ D values were then related directly to their corresponding distances through a constant of proportionality. This constant was determined through a least squares linear fit of the data. Compensation for Airplane Attitude Because the ultrasonic sensor has a nominal beam angle of 15 degrees and reports the closest item from which it receives a valid return signal, the reading from the ultrasonic sensor may not always reflect HAG depending on the attitude of the airplane. If the airplane is either rolled or pitched sharply the signal emitted by the ultrasonic sensor will reflect off the ground at some distance significantly displaced from the point located directly beneath the airplane. The distance from the airplane to this displaced point will be significantly greater than the distance from the airplane to the point directly beneath it. Thus this reading will be significantly greater than the airplane s actual HAG. However, if the attitude of the airplane is known, and the earth can be assumed to be sufficiently flat over the span from the point directly below the airplane to the point at which the signal is reflected, then the actual HAG can still be determined by h actual = h measured c o s ( j j ) c o s ( j j ) + d s i n ; ( 21) where is one half the nominal beam angle, and refer to pitch angle and roll angle respectively, and d is the lateral distance from the center of the airplane at which the sensor is mounted. Equation ( 21) is only valid, and only necessary, when either , , or both are greater than . It is important to note that when this equation must be used to determine HAG, the maximum distance at which true HAG can be accurately determined becomes a dynamic function of airplane attitude rather than a constant ( 12.2 meters for the sensor tested). The airplane attitude dependant value of maximum distance at which true HAG can be accurately determined can be found setting h measured equal to the maximum distance at which true HAG can be determined for level flight, and then solving Equation ( 21) for h actual. Another important consideration when using an ultrasonic sensor to determine HAG is that there will be an angle of incidence associated with every surface at which the ultrasonic signal will not be reflected back toward the sensor. Equation ( 21) is not valid if the angle of incidence is greater than this angle. Instead there is a discontinuity at which the maximum value of true HAG that can be determined goes immediately to zero because no return signal is received. Flight tests have shown that for all outdoor environments tested this angle of incidence is acceptably small such that this discontinuity is not encountered under normal flight conditions. Limitations Due to Propeller Noise A major limitation of the ultrasonic sensor is its tendency to pick up false readings due to noise from the propeller. For the MAGICC Lab airplanes with the motor running at full throttle to near full throttle any reading over 8 meters would be rejected in favor of noise from the propeller. This makes it impossible to trust readings from the ultrasonic sensor at a range greater than 8 meters. For our testing we set the cutoff at which we would trust readings from the ultrasonic sensor at 7.9 meters. This severely limits the usefulness 17 3.5 Optic Flow Based Landing of the ultrasonic sensor by severely limiting its maximum attainable range. This also introduces uncertainty as to whether the sensor will be more or less susceptible to noise from a different combination of motor and propeller. The sensor s maximum range then becomes very platform dependant rather than strictly sensor dependant. The magnitude of the noise which the sensor picks up has been minimized by placing the sensor on the far wingtip of the airplane, as far from the motor and propeller as possible. This, however, introduces a significant payload to the extremity of the airplane causing balancing issues. Also, if the sensor were to be used on an airplane where it could not be mounted a significant distance from the propeller, a significant decrease in maximum range should be expected. Environmental Factors In order for the ultrasonic sensor to determine distances it must transmit and then receive back a signal of acceptable magnitude. For the case of the airplane this signal must be of acceptable magnitude such that it will be accepted by the sensor as the significant signal over noise from the propeller. Testing has shown that this is extremely reliable for flights over surfaces that are sufficiently acoustically reflective such as streets or well trimmed yards. However, for flights over more dense foliage ( fields with tall grass for example) the signal can be sufficiently muffled such that the sensor does not receive an acceptably large return signal until it is 1  2 meters off the ground. Thus maximum attainable range for the ultrasonic sensor is not only dependant on the particular airframe it is flown in ( especially sensor placement on that airframe) and that airplane s attitude, but it is also highly dependant on the acoustical properties of its flight environment. 3.4.3 Results In order for the ultrasonic sensor to facilitate autonomous landing in truly uncontrolled environments where landing altitude relative to takeoff altitude is not known, it must allow the airplane to circle down at a steady rate until true HAG is detected, and then break out of its orbit and follow a glide slope to the landing point. The ultrasonic sensor was not able to reliably perform this task. The sensor s extremely limited range coupled with loss of ultrasonic readings over dense foliage and loss of ultrasonic readings due to airplane attitude were the main contributing factors to the sensor s unreliability. 3.5 Optic Flow Based Landing 3.5.1 Theory Optic flow sensors can be used to determine HAG by relating the flow of features across an imaging array to the speed the imaging array is moving past the surface it is imaging and its distance from that surface. The number of pixels a given object moves in the imaging plane can be combined with data from the airplane s IMU and GPS to determine HAG according to [ 1] h = x 2 t a n fov 2 p n _ t s 2 c o s c o s ; ( 22) where h represents HAG, x represents the displacement of the sensor parallel to the plane being imaged over the sample time t s , represents the average number of pixels of displacement of features across the imaging plane over the same sample time, fov represents the field of view of the sensor, p n is the number of pixels in the imaging array in the direction of motion of the sensor, _ is the average pitch rate of the sensor over the sample period, is the average pitch of the sensor over the sample period, and is the average roll of the sensor over the sample period. 3.5.2 Application Sensor Testing of autonomous landing using optic flow to determine HAG was performed using Agilent s ADNS 2610 optic flow sensor. The ADNS 2610 has a small form factor, measuring only 10 x 12.5 mm and runs at 1500 frames per second. It requires a light intensity of at least 80 mW/ m2 at a wavelength of 639 nm or 100 mW/ m2 at a wavelength of 875 nm. The ADNS 2610 measures the flow of features across an 18 x 18 CMOS imager. It outputs two values, d x and d y , representing the total optic flow across the sensor s field of view in both the x and y directions respectively. These values are stored in a buffer which can store accumulated optic flow up to 128 pixels. Once the buffer is full no more optic flow readings are accumulated until the buffer has been read. Reading the buffer clears it and allows it to begin reaccumulating optic flow measurements. The flow data in the camera y direction corresponds to lateral motion of the airplane and for purposes of establishing HAG can be ignored. The flow data in the camera 18 3.5 Optic Flow Based Landing x direction can be combined with data from the IMU, GPS, and the field of view of the sensor to determine HAG according to Equation ( 22). Optics A critical consideration in outfitting an airplane with an optic flow sensor is the setup of the optics. Narrow angle lenses increase the distance at which the airplane will be able to pick up appreciable optic flow. This increases the distance at which the sensor will be able to accurately measure HAG for a given groundspeed. This comes at the expense of increasing the longitudinal size of the sensor. Narrow angle lenses can also decrease the low end range of the sensor by causing overflow in the optic flow sensor s dx register or by causing the angular flow rate seen by the sensor to become too fast for the sensor to register. Wide angle lenses allow for smaller longitudinal sensor size, but require that larger features be available in the environment. They also decrease the distance at which appreciable optic flow will be recorded for a given groundspeed. A particular advantage of wider angle lenses is that they are less susceptible to noise introduced by angular oscillations of the airplane in pitch and roll. In addition to focal length, lens diameter is an important consideration. This is especially important with regard to the corresponding f stop value of the lens. The f stop is defined as the ratio of a lens focal length to its diameter and is a measure of the amount of light that the optics allow to pass through to the imaging array. Lenses with larger f stops allow less light to pass through. Keeping in mind that the imaging array requires a light intensity of 80  100 W/ m2 in the correct spectral range it is important to select a lens with a low f stop in order to increase the sensor s operational range. Because the sensor automatically adjusts shutter speed to keep the light intensity at an ideal level, there is no disadvantage to selecting a lens with a lower f stop other than increased size. Another important implication of the f stop is that selecting a lens with an increased focal length ( a narrower angle lens) will not only increase the longitudinal size of the sensor, but also increase the lateral size of the sensor if the f stop is to be maintained. For the purposes of these tests lenses were selected which yielded fields of view of 6.5, 2.5, and 1.2 degrees when mounted on the sensor. These lenses had accompanying f stops of 2.0, 2.0, and 2.5 respectively. Each of the configurations is shown in Figure 19. Left FOV: 6.5o Long. Size: 25 mm Lat. Size: 30 mm Mass: 15 g F stop: 2.0 Center FOV: 2.6o Long. Size: 35 mm Lat. Size: 30 mm Mass: 23 g F stop: 2.0 Right FOV: 1.2o Long. Size: 50 mm Lat. Size: 30 mm Mass: 23 g F stop: 2.5 Figure 19: 1.2, 2.5, and 6.5 degree Field of View Optics Configurations. Comparable results were achieved using both the 2.5 degree and 1.2 degree field of view setups. However, the 1.2 degree field of view lens offered the slightly greater range while performing better over feature poor surfaces such as asphalt. The one disadvantage of the 1.2 degree field of view setup up was its increased susceptibility to noise caused by pitch and roll oscillations. Sample Rate Another critical consideration when equipping an airplane with an optic flow sensor is the selection of sample rate. The resolution of the sensor ( measured in units of distance per sensor count) is given by r = l i m h ! 0 0 @ h fov 2 p n t a n 1 V t s 2 ( h h ) t a n 1 V t s 2 h 1 A : ( 23) 19 3.5 Optic Flow Based Landing For a given sensor with a fixed field of view the sensor s resolution is a function of sample rate, groundspeed, and HAG. Plots of sensor resolution for a fixed field of view and varying sample rates and groundspeeds are shown in Figure 20. Figure 20: Plots of Resolution as a function of HAG and Groundspeed and HAG and Sample Rate. For high values of r each sensor counts corresponds to a relatively large distance measure. Correspondingly the noise on the sensor measured in sensor counts is amplified by r into higher magnitude noise in the HAG measurement. This makes large values r undesirable because they can significantly decrease the sensor s signal to noise ratio. It is possible to decrease r by decreasing the sample rate as shown in Figure 20. This, however, can be undesirable because it decreases the sensor s response time and can lead to overflow in the sensor s d x register. Furthermore, for a fixed sample rate, maintaining a low value of r for relatively large values of HAG results in unnecessarily low values of r for lower values of HAG. This means unnecessarily low sample rates are being maintained where significantly faster sample rates could be used without significant loss in signal to noise ratio. Such an approach also results in premature sensor overflow and corresponding loss of low end range. To solve these problems two approaches were developed using a dynamic rather than static sample rate. The first approach maintains an r value along a specified resolution curve regardless of groundspeed of the plane. The second approach uses the current velocity estimate and HAG estimate to maintain a constant value of r regardless of groundspeed and HAG. Constant Curve Resolution In general, flight control systems for UAVs utilize a feedback loop to control around a commanded airspeed rather than groundspeed. Even when the feedback loop is closed around groundspeed, the commanded groundspeed may or may not be realizable due to wind conditions. In either case in the presence of wind it is reasonable to assume that groundspeed may vary significantly depending on the course of the airplane relative to the wind direction. As shown in Figure 20 variations in groundspeed can alter the resolution curve significantly. However, if the sample rate is set according to t s = V ; ( 24) where represents the sample rate divisor, the value of r will follow a fixed resolution curve as a univariate function of HAG. This is useful because it eliminates the variations in resolution due to velocity and provides reliable repeatable resolutions regardless of wind conditions or commanded airspeeds. Constant Value Resolution When an estimate of the current HAG is available, the correct sample rate can be calculated such that the sensor resolution can be maintained constant regardless of HAG or groundspeed. This sample rate is found by solving Equation ( 23) for t s . The result is t s = l i m h ! 0 0 B B @ 1 V 2 ( h h ) V 2 h 1 t a n h fov 2 p n r V 2 4 h ( h h ) 1 C C A ; ( 25) where r now represents a user defined parameter that defines the resolution of the sensor for all combinations of HAG and groundspeed. In practice, the actual value of t s is capped at both the high and low end, but allowed to change in order to maintain a fixed resolution inside of some saturation block. In flight tests using a dynamic sample rate to fix sensor resolution cut noise drastically for high values of HAG while significantly increasing the sample rate and thus 20 3.5 Optic Flow Based Landing improving the response time for low values of HAG. In practice Equation ( 25) is implemented by dropping the limit expression and fixing h to some reasonably small number. Flight test results using both constant curve resolution and constant value resolution are shown in Figure 21. The spike in the optic flow measurement at sample 400 in the constant value resolution plot and 275 in the constant curve resolution plot corresponds to flight over the asphalt road. In practice such spikes are ignored by ignoring values for which the optic flow sensors reports poor surface quality. Figure 21: Optic HAG Results from Glideslope Descents Using Fixed Value Resolution ( Left) and Fixed Curve Resolution ( Right). 3.5.3 Results Ranging The accuracy of the distance measurements calculated from optic flow was tested by comparison with a laser rangefinder. Testing was performed by mounting the optic flow sensor and laser rangefinder perpendicular to the motion of a vehicle driving along a freeway at varying distances from the freeway sound barrier wall. Range values were recorded for both the optic flow calculations and the laser rangefinder. The results of this test are shown in Figure 22. The results show very close agreement and validate the readings obtained from the optic flow calculations. Figure 22: Range Data from Optic Flow Calculations Compared to Laser Rangefinder Data. Range data from the optic flow sensor in flight tests showed reliable, repeatable readings for any height above ground less than 40 meters. For the purposes of landing based on HAG calculations from optic flow, values above 40 meters were discarded as noise. 21 3.6 Conclusion Landing Using the fixed value resolution sample rate and the 2.65 degree field of view optics, HAG measurements based on optic flow were used to perform twenty seven consecutive autonomous landings from a single approach point. The results are shown in Figure 19. For each of these autonomous landings HAG was calculated according to h = h bar + h offset ( 26) h offset = ( h bar h optic ) + ( 1 ) h offset : ( 27) The value of h was updated every time through the loop and the value of h offset was updated every t s seconds where t s varied with the current HAG estimate and groundspeed. The readings from the barometric pressure sensor were allowed to drift over the twenty seven landings in order to simulate an unknown landing point altitude. The barometric altimeter drifted to 15 meters of inaccuracy over the first 22 landings and was then rezeroed to introduce negative offsets between true HAG and barometric altitude. The value of h offset ( corresponding to the drift in the barometric altimeter) for each of the landings is shown in the third subplot of Figure 19. The second subplot of Figure 23 shows the distance from the landing point for each of the twenty seven landings. Figure 23: Results of 27 Consecutive Optic Flow Based Autonomous Landings. The average distance between the desired and actual touchdown points for autonomous landings using optic flow data was 4.3 meters. This represents a 56% improvement over the average distance from the desired touchdown point using only the barometric altimeter. Landings based on HAG from optic flow also had a much narrower spread than those based on barometric altitude alone with a standard deviation of 2.2 for optic flow based landings compared to 5.4 for barometric. Because the circle down portion of the landing routine is based on barometric altitude until HAG readings fall under a threshold at which they are trusted, it is important to note that the distance from the desired touchdown point does not increase with increasing error in barometric altitude. This demonstrates the ability of the optic flow sensor to land the airplane accurately even when the altitude of the landing point relative to the takeoff point is not known. 3.6 Conclusion The first step necessary to land a vehicle autonomously is the design of control algorithms that will allow the vehicle to accurately perform a landing maneuver. A method for controlling the landing maneuver based on heading fields has been presented. This method has been tested both in simulation and in hardware and has proven effective, repeatable, and robust to wind. The method has been used effectively to perform hundreds of autonomous landings in a wide variety of wind conditions. 22 In addition to effectively tracking a landing maneuver, in order to accomplish true autonomous landing in uncontrolled environments, the UAV must be able to detect HAG reliably. In environments where the terrain is sufficiently flat and the altitude of the landing point relative to the takeoff point is known, a well calibrated barometric altimeter ( with both airflow and temperature compensation) is sufficient to accurately determine HAG and land the UAV. In environments where this is not the case fusion of data from the barometric altimeter and the optic flow sensor provides a feasible and robust method for estimating HAG. This estimate of HAG can be combined with robust landing algorithms to successfully land mini and micro fixed wing UAVs within meters of their desired landing point. No other sensor combination to date with a similar payload and computational cost has demonstrated the same capability. 4 Automatic Trim and Gain Tuning Traditional autopilot design involves taking the equations of flight and linearizing them about an operating point. Controllers are then designed with these linearized models. A popular choice for such controllers are PID controllers. Although such designs are easy to implement and can work well, this method requires that the plane remain within the design flight conditions and are typically not very failure tolerant. To allow for a larger flight envelope, gain scheduling can be used to bridge different PID controllers designed around different operating points. However such designs still require tuning at the different operating points, do not address fault tolerance, and must be tuned for different platforms. Adaptive control was developed to address some of these problems. Examples of different adaptive control approaches are neural network, least squares estimation, and Lyapunov based methods. One adaptive approach is to use neural networks. This often entails training a network offline which is used in the autopilot. However, to compensate for training errors and to adapt to changes in the aircraft ( i. e. battle damage, system failures, etc.), a separate online neural network is also used. This online network is then constantly adapting to changes on the plane and training error. [ 2], [ 3], [ 4], [ 5], and [ 6] discuss this approach. Another approach identifies the plane parameters using recursive least squares techniques. These parameters are then used to adjust the controller. Such controllers have the ability to quickly converge on the plane parameters [ 7], and therefore adapt quickly to failures. [ 8], [ 9], and [ 5] use least squares estimation. Also, a Lyapunov based approach can be used. In this approach parameters are adjusted according to a Lyapunov function to assure stability. These algorithms typically have a small code footprint and require little computing power. [ 10], [ 11], and [ 12] are examples of Lyapunov based algorithms. In our research, we deal with Miniature Aerial Vehicles ( MAVs). We have found that PID values are different for each MAV, even if they have the same geometry. This is mostly due to the manufacturing process. The desire to have one autopilot that works on any of our MAVs with only one set of gains has lead us into adaptive control. Because of the size of our MAVs, we need to use a small autopilot. This means we have code space, memory, and computational restrictions on our autopilot. Also, because of the differing platforms and an imperfect manufacturing process, we do not have an accurate model of our MAVs. This means we are not able to use neural networks because they require a good model of the MAV, require some computational power to implement, and are not well suited for being placed on MAVs of differing geometries. We also cannot use the least squares algorithms which use inverse matrices that require memory and computer power to handle. Such resources come at a premium. This has lead us to use the Lyapunov based approach. Using model reference design and Lyapunov theory, we will show that a stable, adaptive algorithm can be derived that is computationally simple. This new autopilot design will automatically adapt itself to whatever MAV it is flying. Its inputs will be desired pitch and roll angles. We assume that both heading and altitude looped are closed outside of the current control design and feeding the autopilot with the desired pitch and roll angles. It should be noted that adaptive control algorithms are primarily used to compensate for actuator and control surface failure. To the authors knowledge it is not being used to make a universal autopilot that can be used on a variety of different platforms. This makes our approach somewhat unique. Further research should be done to see how many different platforms our autopilot can fly on. 23 4.1 UAV Model In this paper we will show four different Model Reference Adaptive Control ( MRAC) schemes. All four schemes will use the same reference model, but will differ on how much is lumped into the regressor s bias term. The first scheme will employ backstepping to handle derivatives on angular rates. The second does not use backstepping; instead we assume that the angular rate derivatives change slowly enough that we can assume they are slow varying parameters. The third and forth schemes further simplify by lumping more slowly changing signals into the bias term. Flight data from three different MAV designs will also be shown. 4.1 UAV Model If we assume that the inertia matrix has the form J = 0 @ J x 0 J x z 0 J y 0 J x z 0 J z 1 A ; and that the aerodynamic moments are linear, then the rotational dynamics of the UAV are given by _ = p + q s i n t a n + r c o s t a n ( 28) _ = q c o s r s i n ( 29) _ = q s i n s e c + r c o s s e c ( 30) p _ = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V + C p a a + C p r r ( 31) q _ = J x z J y ( r 2 p 2 ) + J z J x J y p r + 1 2 J y V 2 c S C m 0 + C m + C m q c q V + C m e e ( 32) r _ = 3 p q 4 q r + 1 2 V 2 S b 2 C r 0 + C r + C r p b p 2 V + C r r b r 2 V + C r a a + C r r r ; ( 33) where , , roll, pitch, yaw angles d , d desired roll and pitch angles m , m reference model roll and pitch angles k m reference model gain p , q , r roll, pitch, yaw rates , angle of attack, sideslip angle , flight path heading angle, flight path angle e , a , r elevator, aileron, rudder commands V airspeed S , b , c wing area, wing span, average cord length J i j inertial moments about i and j axies Table 1: Definition of terms 24 4.1 UAV Model = ( J x J z J 2 x z ) 1 = J x z ( J x J y + J z ) = ; 2 = ( J z ( J z J y ) + J 2 x z ) = 3 = ( ( J x J y ) J x + J 2 x z ) = 4 = ( J x z ( J x J y + J z ) ) = C p 0 = J z C 0 + J x z C n 0 C p = J z C + J x z C n C p p = J z C p + J x z C n p C p r = J z C r + J x z C n r C p a = J z C a + J x z C n a C p r = J z C r + J x z C n r C r 0 = J x z C 0 + J x C n 0 C r = J x z C + J x C n C r p = J x z C p + J x C n p C r r = J x z C r + J x C n r C r a = J x z C a + J x C n a C r r = J x z C r + J x C n r : For all schemes we use the following first order models, _ m = k m d m ( 34) for pitch and _ m = k m d m ( 35) for roll. We define tacking error as = m and ( 36) = m ( 37) for pitch and roll, respectively. The derivatives are _ = _ _ m = q c o s r s i n _ m and ( 38) _ = _ _ m = p + q s i n t a n + r c o s t a n _ m : ( 39) We also make the following assumptions: 25 4.2 Scheme A: nonlinear MRAC using backstepping A1: The following states can be measured: , , p , r , q . A2: The pitch angle is limited to where < = 2 . A3: The roll angle is limited to where < = 2 . A4: The pitch angle = + where is the flight path climb angle. A5: The pitch angle = where is the flight path heading angle. A6: For the pitch controllers is constant. A7: For the roll controllers is constant. 4.2 Scheme A: nonlinear MRAC using backstepping 4.2.1 Pitch Attitude Hold For pitch attitude hold, we will work with Equations ( 29) and ( 32). Our objective is to design a controller so that the pitch angle follows the reference model ( 34). Our desire is to pick e such that Equation ( 38) becomes _ = 1 ; where 1 > 0 . If we define q d e s 4 = _ m 1 + r s i n c o s ; and let q = q d e s in Equation ( 38), we get the desired result. Using the actual pitch rate but adding and subtracting q d e s we obtain _ = q c o s r s i n _ m = q + q d e s q d e s c o s r s i n _ m = q d e s c o s r s i n _ m + q q d e s c o s = 1 + q c o s ; ( 40) where q 4 = q q d e s . Differentiating q we get _ q = _ q _ q d e s = J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y C m 0 + C m ( ) + C m q c q V + C m e e q _ d e s : Define the Lyapunov function candidate V = 1 2 2 + 1 2 q 2 : 26 4.2 Scheme A: nonlinear MRAC using backstepping Differentiating we get _ V = _ + q _ q = 1 2 + q c o s + q _ q = 1 2 + q q _ d e s + J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y h C m 0 + C m ( ) + C m q c q V + C m e e i : ( 41) If we knew the coefficients, we would like to pick e to satisfy 2 q = q _ d e s + J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y h C m 0 + C m ( ) + C m q c q V + C m e d e s e i : Solving for d e s e gives d e s e = J y 1 2 V 2 S c C m e 2 q c o s + q _ d e s + 1 1 2 V 2 S c C m e J x z ( p 2 r 2 ) ( J z J x ) p r 1 C m e C m 0 + C m ( ) + C m q c q V : ( 42) If we define K 4 = 0 B B B B B B B B B B B @ J y 1 2 S c C m e J x z 1 2 S c C m e J z J x 1 2 S c C m e C m 0 C m C m e C m C m e C m q c C m e 1 C C C C C C C C C C C A and 4 = 0 B B B B B B @ 2 q c o s + q _ d e s V 2 p 2 r 2 V 2 p r V 2 1 q V 1 C C C C C C A ; ( 43) then d e s e = T K : Let ^ K be the current estimate of K , then we will use the control e = T ^ K : Therefore, plugging into Equation ( 41) we get _ V = 1 2 + q q _ d e s + J x z J y ( p 2 r 2 ) + J z I x J y p r + 1 2 V 2 S c J y h C m 0 + C m + C m q c q V + C m e e + d e s e d e s e i = 1 2 2 q 2 + q 1 2 V 2 S c J y C m e ( e d e s e ) = 1 2 2 q 2 + 1 2 S c C m e J y ! V 2 q T ( ^ K K ) = 1 2 2 q 2 1 2 S c C m e J y ! V 2 q T K ; 27 4.2 Scheme A: nonlinear MRAC using backstepping where K 4 = K ^ K . To cancel the last term, let the modified Lyapunov function candidate be V 2 = 1 2 2 + 1 2 q 2 + 1 2 V S c C m e I y ! 1 2 K T 1 K : ( 44) where is a diagonal matrix. Then _ V 2 = 1 2 2 q 2 1 2 S c C m e I y ! V 2 q T K + 1 2 V S c C m e I y ! _ K T 1 K : Observing that _ K = _ K _ ^ K = _ ^ K since K is constant, gives _ V 2 = 1 2 2 q 2 1 2 V S c C m e I y ! h 1 _ ^ K + V q i T K : ( 45) Therefore letting _ ^ K = V q ; ( 46) gives _ V 2 = 1 2 2 q 2 : The algorithm can be summarized as follows: Algorithm 1 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: q _ d e s 4: Compute: = 2 q c o s + q _ d e s V 2 p 2 r 2 V 2 p r V 2 1 q V T , 5: Update the gain estimate according to: _ ^ K = V q 6: Compute the elevator command: e = T K . 28 4.2 Scheme A: nonlinear MRAC using backstepping 4.2.2 Roll Attitude Hold Using Equation ( 35) as our roll model reference, we want Equation ( 39) to be _ = , so we define p d e s to be p d e s 4 = q s i n t a n r c o s t a n + _ m : We define p as p 4 = p p d e s ( 47) and take its derivative _ p = _ p _ p d e s : ( 48) Substituting equation ( 31) into equation ( 48) we get the following. _ p = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V + C p a a + C p r r p _ d e s ( 49) We define our Lyapunov candidate function V = 1 2 2 + 1 2 p 2 ( 50) and take it derivative _ V = _ + p _ p = p + q s i n t a n + r c o s t a n _ m + p _ p : ( 51) We add and subtract p d e s to get _ V = p + p d e s + q s i n t a n + r c o s t a n _ m + p _ p = 1 + p + p _ p = 1 2 + p + p _ p = 1 2 + p + _ p : ( 52) To arrive at the Lyapunov equation we want, we set + _ p = 2 p . If we now plug this into equation ( 49) and solve for a , calling it d e s a , gives us d e s a = 4 V 2 S b C p a 2 p + p _ d e s 4 1 V 2 S b C p a p q 4 2 V 2 S b C p a q r C p 0 C p a C p C p a C p p C p a b p 2 V C p r C p a b r 2 V C p r C p a r : ( 53) Using the assumption that = where is the flight path heading angle gives us d e s a = 4 V 2 S b C p a 2 p + p _ d e s 4 1 V 2 S b C p a p q 4 2 V 2 S b C p a q r C p 0 C p a C p C p a + C p C p a C p p C p a b p 2 V C p r C p a b r 2 V C p r C p a r : ( 54) which can be represented as d e s a = T K ( 55) 29 4.2 Scheme A: nonlinear MRAC using backstepping where = 0 B B B B B B B B B B @ 2 p + p _ d e s V 2 p q V 2 q r V 2 1 p V r V r 1 C C C C C C C C C C A K = 0 B B B B B B B B B B B B B B B B B @ 4 b V 2 S C a 4 1 b V 2 S C p a 4 2 b V 2 S C p a C p 0 C p a + C p C p a C p C p a b C p p 2 C p a b C p r 2 C p a C p r C p a 1 C C C C C C C C C C C C C C C C C A : ( 56) However, we do not know the terms in K and instead only have an estimate ^ K so really we have a = T ^ K . Adding and subtracting d e s a to equation ( 49) we get _ p = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V ( 57) + C p a a d e s a + d e s a + C p r r p _ d e s _ p = V 2 S b C p a 4 a d e s a 2 p _ p = V 2 S b C p a 4 T ^ K K 2 p _ p = V 2 S b C p a 4 T K 2 p + _ p = V 2 S b C p a 4 T K 2 p ( 58) Defining a new Lyapunov function candidate V 2 = 1 2 2 + 1 2 p 2 + V S b C p a 8 K T 1 K ( 59) where is a diagonal matrix. We differentiate the Lyapunov function to get V _ 2 = _ + p _ p + V S b C p a 4 _ K T 1 K = 1 2 + p + p _ p + V S b C p a 4 _ K T 1 K ( 60) Substituting the results from ( 57) we get V _ 2 = 1 2 + p + _ p + V S b C p a 4 _ K T 1 K = 1 2 + p V 2 S b C p a 4 T K 2 p + V S b C p a 4 _ K T 1 K = 1 2 2 p 2 + p V 2 S b C p a 4 T K + V S b C p a 4 _ K T 1 K = 1 2 2 p 2 + V S b C p a 4 p V T + _ K T 1 K : ( 61) Setting p V T + _ K T 1 = 0 we arrive at the update law _ ^ K = p V : ( 62) The algorithm can be summarized as follows. 30 4.3 Scheme B: nonlinear MRAC without backstepping Algorithm 2 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: p _ d e s 4: Compute: = 2 p + p _ d e s V 2 r _ + p q V 2 q r V 2 1 p V r V r T , 5: Update the gain estimate according to: _ ^ K = p V 6: Compute the aileron command: a = T ^ K . 4.3 Scheme B: nonlinear MRAC without backstepping 4.3.1 Pitch Attitude Hold Using the same definition of tracking error for pitch as above, we design controller without using backstepping. Again, our desire is to pick e such that Equation ( 38) becomes _ = 1 ; ( 63) where 1 > 0 . Solving for q in equation ( 32), and assuming q _ is slow varying parameters, we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 64) If we knew the parameters we could we could command d e s e = T K where T = 0 B B B B @ 1 k m ( d m ) r s i n V c o s 1 p 2 r 2 p r 1 C C C C A K = 0 B B B B B B B B @ c C m p C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e J x z 1 2 V 2 S c C m e J z J x 1 2 V 2 S c C m e C m C m e 1 C C C C C C C C A ( 65) Adding and subtracting this in equation ( 64) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 66) 31 4.3 Scheme B: nonlinear MRAC without backstepping Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 67) We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 68) We summarize the algorithm. Algorithm 3 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: T = 1 k m ( d m ) r s i n V c o s 1 p 2 r 2 p r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the elevator command: e = T ^ K . 32 4.3 Scheme B: nonlinear MRAC without backstepping 4.3.2 Roll Attitude Hold Assuming p _ and r _ are slow varying parameters and solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 69) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 70) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 71) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 72) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us = 0 B B B B B B B B B B @ 1 + q s i n t a n + r c o s t a n k m ( d m ) V 1 V 2 p q V 2 q r V 2 1 r V r 1 C C C C C C C C C C A and K = 0 B B B B B B B B B B B B B B B B @ b C p p 2 C p a 4 S b p _ 4 1 S b 4 2 S b C p 0 C p a C p C p a C p C p a C p r 2 C p a C p r C p a 1 C C C C C C C C C C C C C C C C A : ( 73) 33 4.4 Scheme C: nonlinear MRAC grouping terms Adding and subtracting d e s a from Equation ( 70) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 74) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 75) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 76) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 4 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: = 1 + q s i n t a n + r c o s t a n k m ( d m ) V 1 V 2 p q V 2 q r V 2 1 r V r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 4.4 Scheme C: nonlinear MRAC grouping terms 4.4.1 Pitch Attitude Hold We assume that q _ , 1 c o s , p 2 r 2 , and p r terms are slow varying parameters. Our desire is to pick e such that Equation ( 38) is _ = 1 ; ( 77) 34 4.4 Scheme C: nonlinear MRAC grouping terms where 1 > 0 . Solving for q in equation ( 32) we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 78) If we knew the parameters we could we could command d e s e = T K where T = 0 @ 1 k m ( d m ) r s i n V 1 1 A and K = 0 B B @ c C m p c o s C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e + J x z 1 2 V 2 S c C m e p 2 r 2 J z J x 1 2 V 2 S c C m e p r C m C m e 1 C C A : ( 79) Adding and subtracting this in equation ( 78) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 80) Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 81) We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 82) We summarize the algorithm. 35 4.4 Scheme C: nonlinear MRAC grouping terms Algorithm 5 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d m ) V c o s 1 T , 4: Update the gain estimate according to: _ ^ K = V c o s 5: Compute the elevator command: e = T ^ K . 4.4.2 Roll Attitude Hold Let the derivative of the tracking error be _ = p + q s i n t a n + r c o s t a n _ m : ( 83) We make the assumption here that the terms p _ , r _ , q s i n t a n + r c o s t a n , p q , q r , and r V are slow varying parameters Solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 84) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 85) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 86) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 87) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us = 0 B B @ 1 k m ( d m ) V 1 r 1 C C A and K = 0 B B B B B @ b C p p 2 C p a C p b i a s C p C p a C p r C p a 1 C C C C C A ; ( 88) 36 4.4 Scheme C: nonlinear MRAC grouping terms where C p b i a s = 4 V 2 S b p _ + C p 0 C p a + b C p p 2 C p a ( q s i n t a n + r c o s t a n ) 4 1 S b p q V 2 + 4 2 S b q r V 2 C p C p a C p r 2 C p a r V Adding and subtracting d e s a from Equation ( 89) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 89) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 90) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 91) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 6 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: = 1 k m ( d ) 2 V 1 r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 37 4.5 Scheme D: nonlinear MRAC simple 4.5 Scheme D: nonlinear MRAC simple 4.5.1 Pitch Attitude Hold We assume that q _ , 1 c o s , p 2 r 2 , and p r terms are slow varying parameters. We want to pick e such that _ = 1 ; ( 92) where 1 > 0 . Solving for q in equation ( 32) we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 93) If we knew the parameters we could we could command d e s e = T K where T = 1 k m ( d m ) V 1 ! T and K = 0 B B @ c C m p c o s C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e ( ) + J x z 1 2 V 2 S c C m e p 2 r 2 J z J x 1 2 V 2 S c C m e p r C m C m e 1 C C A : ( 94) Adding and subtracting this in equation ( 95) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 95) Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 96) 38 4.5 Scheme D: nonlinear MRAC simple We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 97) We summarize the algorithm. Algorithm 7 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d m ) V 1 T , 4: Update the gain estimate according to: _ ^ K = V c o s 5: Compute the elevator command: e = T ^ K . 4.5.2 Roll Attitude Hold We make the assumption here that the terms p _ , r _ , q s i n t a n + r c o s t a n , p q , q r , and r V are slow varying parameters Solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 98) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 99) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 100) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 101) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us T = 0 @ 1 k m ( d ) V 1 r 1 A T and K = 0 B B @ b C p p 2 C p a C p b i a s C p r C p a 1 C C A ; ( 102) 39 4.5 Scheme D: nonlinear MRAC simple where C p b i a s = 4 V 2 S b p _ + C p 0 C p a + b C p p 2 C p a ( q s i n t a n + r c o s t a n ) 4 1 S b p q V 2 + 4 2 S b q r V 2 C p C p a ( ) C p r 2 C p a r V Adding and subtracting d e s a from Equation ( 99) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 103) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 104) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 105) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 8 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d ) V 1 T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 40 4.6 Flight Data 4.6 Flight Data 4.6.1 Platform The algorithms were flown on MAVs equipped with the Kestrel Autopilot running a 29 MHz Rabbit microcontroller with 512K Flash and 512K RAM. The autopilot is equipped with a variety of sensors, including rate gyros, accelerometers, an absolute pressure sensor for measuring altitude, a differential pressure sensor for measuring airspeed, and a GPS receiver. The autopilot is programmed in C. The adaptive control algorithms run at about 80 90 hertz with this platform. The airframes flown were a Zagi style flying wing design ( Deliverance), a larger and heavier flying wing design ( Ray), and a fixed wing glider design ( Phidipides). In the initial phase of testing, we used Deliverance to tune adaptive control gains and to obtain a comparison between the four algorithms presented in this paper. Afterward, we narrowed our testing to one algorithm ( D) and flew with the other airframes. Deliverance Weight: 0.734 k g Wingspan: 1.207 m Mean geometric chord: 0.254 m Ray Weight: 1.321 k g Wingspan: 1.511 m Mean geometric chord: 0.305 m Phidipides Weight: 0.969 k g Wingspan: 1.524 m Mean geometric chord: 0.168 m Nose to tail: 0.914 m 41 4.6 Flight Data Pitch k m 3.0 Pitch 1 45.0 Pitch 1 0.06 Pitch 1 0.01 Pitch Leakage gain 1 0.001 Pitch Leakage gain 2 0.0005 Roll k m 4.0 Roll 1 140.0 Roll 1 0.005 Roll 1 0.001 Roll Leakage gain 1 0.001 Roll Leakage gain 2 0.0001 Table 2: Adaptive control gains used in Algorithm D 4.6.2 Implementation While Lyapunov theory provides for stability, adaptive control engineers have several techniques available to ensure good stability in practice. Many times, flight conditions do not provide sufficient excitation for all of the adaptive parameters. Sensor noise and other disturbances can also cause parameter drift. [ 13] Some of the parameters may diverge such that when flight conditions change quickly ( i. e., we switch from straight, level flight into a climbing turn), the plane can go unstable. We use three techniques in our implementation to improve stability: dead zone, leakage, and normalization. With dead zone, adaption is stopped once the actual value gets within a certain region around the desired value. Once we get into a flight pattern with little excitation, the ^ K parameters will adapt until the plane is close to the desired value. Then adaption freezes until flight conditions change such that we are no longer within the no adaption ( dead ) zone. Leakage subtracts off a small percentage of each ^ K parameter each time the code is executed. This prevents parameters from ever becoming too large. So as not to affect performance, the leakage gain is kept small, about 0.0001 for most parameters. Normalization is similar to leakage. Each time through the control loop, all of the ^ K parameters are normalized with respect to themselves. This prevents one parameter from become much larger than the others, but preserves the relative weighting between each parameter. Because all of our MAVs do not have a rudder we have implemented the controller with r = 1 . 4.6.3 Results We wanted to compare the four algorithms against each other and against the usual PD controller used in our lab to control pitch and roll on the MAVs. To obtain this comparison, we flew an hourglass shaped pattern with legs between 200 and 300 meters long on the same airframe ( Deliverance) with each of the four algorithms. The flights incorporated auto takeoff and auto land functions, giving the adaptive algorithms full control for the entire flight. We also obtained data for the flight plan flown with more conventional PD control, with gains tuned by hand. For all of the adaptive algorithms, transients die out within about 3 seconds for pitch and about 10 seconds for roll. Along with the adaptive control plots, we show the values of the ^ K parameters adapting with time. In most cases, these parameters settle to a steady value after adaption has been running for a while. Note that while the parameters represent certain combinations of aerodynamic factors and other estimated quantities, they do not, in most cases, converge to these true values. The Lyapunov stability proof used to derive the algorithms does not guarantee that these parameters will converge to the true values. Rather, the proof shows that the parameters stabilize to values that will cause the tracking error to go to zero. The results show that each adaptive control algorithm s performance is comparable to PD control. Before each algorithm could be flown, several adaptive control gains ( namely k m , 1 , and ) had to be tuned. Algorithm D, with the fewest number of gains, was easiest to tune. Algorithm D also is the simplest algorithm to implement and takes up the smallest amount of autopilot code space. 42 4.7 Conclusion Algorithm Pitch Roll A 2.0 2.5 B 2.0 3.5 C 2.0 5.2 D 2.0 3.9 PD control 1.5 3.5 Table 3: Average error in degrees for each algorithm  from last third of hourglass flight tests on Deliverance We used average error as a metric to compare between the algorithms. We define average error as the average of p ( d e s i r e d a c t u a l ) 2 over the sample. Each figure included in this paper displays average error over the total sample as well as the average error over each third of the flight time. This shows the gradual improvement in performance characteristic of adaptive control as the parameters become more finely tuned. In most of the figures, the last third of the flight has less average error than the first third. The longer the adaptive algorithms run, the better the performance will be. The adaptive gains k m , 1 , 2 , and were very difficult to tune for most of the algorithms, especially algorithm C. Algorithms A and B were tuned to give good performance, but did not have very good stability. Based on low average error, good stability, and simplicity of implementation, we chose to focus on algorithm D for further testing. None of the algorithms outperformed PD control, but bear in mind that the control gains were finely tuned to Deliverance s airframe. To test the robustness of adaptive control, we attached a flap near the center of Deliverance s right wing and deployed it mid flight. By doing this, we simulated what a change in the airframe may have on the controllability of the MAV. With the flap deployed, PD control can still fly the MAV, but with substantial steady state error on roll. However, algorithm D can be seen to adapt to the change in aerodynamics within a few seconds. When the flap goes back down, the parameters once again adapt back to their previous states. Flap attached to Deliverance To further test adaptive control s robustness, we flew two very different airframes; Ray, a heavier Zagi style airframe, and Phidipides, a fixed wing glider design. The plots presented in this paper are flight data from flying Algorithm D with the adaptive gains tuned while flying Deliverance. Both Ray and Phidipides fly well with these same gains. Phidipides in particular was very hard to tune PD gains for, but algorithm D had no problems. This shows the main strength of adaptive control: the ability to fly the same code on different airframes without having to hand tune gains. 4.7 Conclusion In this paper we have shown the derivation for four adaptive control schemes for pitch and roll. We also present flight test data for each of the schemes. Through comparing performance, stability, and simplicity, Algorithm D was chosen as the best algorithm to pursue for further testing and development on the autopilot. We show that the algorithm performs well on several different airframes without having to tune any gains. 43 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm A ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 24: Algorithm A pitch on Deliverance flying an hourglass 44 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm A ( Pitch) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 25: Algorithm A pitch parameters on Deliverance flying an hourglass 45 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm A ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 26: Algorithm A roll on Deliverance flying an hourglass 46 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm A ( Roll) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Figure 27: Algorithm A roll parameters on Deliverance flying an hourglass 47 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm B ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 28: Algorithm B pitch on Deliverance flying an hourglass 48 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm B ( Pitch) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 29: Algorithm B pitch parameters on Deliverance flying an hourglass 49 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm B ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 30: Algorithm B roll on Deliverance flying an hourglass 50 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm B ( Roll) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 31: Algorithm B roll parameters on Deliverance flying an hourglass 51 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm C ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 32: Algorithm C pitch on Deliverance flying an hourglass 52 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm C ( Pitch) Controlling Deliverance 50 100 150 200 250 300 350 50 100 150 200 250 300 350 Time ( seconds) Figure 33: Algorithm C pitch parameters on Deliverance flying an hourglass 53 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm C ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 34: Algorithm C roll on Deliverance flying an hourglass 54 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm C ( Roll) Controlling Deliverance 50 100 150 200 250 300 350 50 100 150 200 250 300 350 Time ( seconds) Figure 35: Algorithm C roll parameters on Deliverance flying an hourglass 55 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm D ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 36: Algorithm D pitch on Deliverance flying an hourglass 56 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm D ( Pitch) Controlling Deliverance 50 100 150 200 250 300 350 Time ( seconds) Figure 37: Algorithm D pitch parameters on Deliverance flying an hourglass 57 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm D ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 38: Algorithm D roll on Deliverance flying an hourglass 58 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm D ( Roll) Controlling Deliverance 50 100 150 200 250 300 350 Time ( seconds) Figure 39: Algorithm D roll parameters on Deliverance flying an hourglass 59 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) PD ( Pitch) Controlling Deliverance Desired Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 40: PD pitch on Deliverance flying an hourglass 60 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) PD ( Roll) Controlling Deliverance Desired Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 41: PD roll on Deliverance flying an hourglass 61 4.7 Conclusion 50 100 150 200 250 Algorithm D ( Pitch) Controlling Deliverance as Flap Deploys Desired Model Actual 50 100 150 200 250 Algorithm D ( Roll) Controlling Deliverance as Flap Deploys Desired Model Actual 50 100 150 200 250 Flap Time ( seconds) Figure 42: Algorithm D Pitch and Roll with Flap Deployment flying an hourglass. One transit around the hourglass is completed with the flap deployed. 62 4.7 Conclusion 50 100 150 200 250 Estimated Parameters of Algorithm D Controlling Deliverance as Flap Deploys 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 Flap Time ( seconds) Figure 43: Algorithm D Pitch and Roll parameters with Flap Deployment flying an hourglass. Notice k 2 for roll adapting to a new value when the flap is deployed and adapting back to the previous value when the flap is stowed. 63 4.7 Conclusion 50 100 150 200 250 PD ( Pitch) Controlling Deliverance as Flap Deploys Desired Actual 50 100 150 200 250 PD ( Roll) Controlling Deliverance as Flap Deploys Desired Actual 50 100 150 200 250 Flap Time ( seconds) Figure 44: PD Control with Flap Deployment flying an hourglass. Notice the steady state error of about 10 degrees in roll. 64 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Time ( seconds) Algorithm D ( Pitch) Controlling Ray 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Desired Model Actual Figure 45: Algorithm D pitch on Ray flying between two waypoints 65 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Estimated Parameters of Algorithm D ( Pitch) Controlling Ray 20 40 60 80 100 120 140 160 180 Time ( seconds) Figure 46: Algorithm D pitch parameters on Ray flying between two waypoints 66 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Time ( seconds) Algorithm D ( Roll) Controlling Ray Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 47: Algorithm D roll on Ray flying between two waypoints 67 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Estimated Parameters of Algorithm D ( Roll) Controlling Ray 20 40 60 80 100 120 140 160 180 Time ( seconds) Figure 48: Algorithm D roll parameters on Ray flying between two waypoints 68 4.7 Conclusion 50 100 150 200 250 300 350 400 Time ( seconds) Algorithm D ( Pitch) Controlling Phidipidies Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 49: Algorithm D pitch on Phidipides flying an hourglass 69 4.7 Conclusion 50 100 150 200 250 300 350 400 Estimated Parameters of Algorithm D ( Pitch) Controlling Phidipides 50 100 150 200 250 300 350 400 Time ( seconds) Figure 50: Algorithm D pitch parameters on Phidipides flying an hourglass 70 4.7 Conclusion 50 100 150 200 250 300 350 400 Time ( seconds) Algorithm D ( Roll) Controlling Phidipidies Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 51: Algorithm D roll on Phidipides flying an hourglass 71 4.7 Conclusion 50 100 150 200 250 300 350 400 Estimated Parameters of Algorithm D ( Roll) Controlling Phidipides 50 100 150 200 250 300 350 400 Time ( seconds) Figure 52: Algorithm D roll parameters on Phidipides flying an hourglass 72 5 Students Funded Under Project Student Topic Steve Griffiths Magnetometer integration Blake Barber Height above ground sensor Andrew Eldridge Magnetometer integration Josh Matthews Adaptive control Nathan Knoebel Adaptive control Steve Osborne Adaptive control References [ 1] J. Chahl, M. V. Srinivasan, and S. W. Zhang, Landing strategies in honeybees and applications to uninhabited airborne vehicles, The International Journal of Robotics Research, vol. 23, no. 2, pp. 101 110, 2004. [ 2] D. Shin and Y. Kim, Reconfigurable flight control system design using adaptive neural networks, 12, no. 1, Jan. 2004. [ 3] R. L. Broderick, Statistical and adaptive approach for verification of a neural based flight control system, in Digital Avionics Systems Conference, vol. 2, pp. 6. E. 1 1 6. E. 1 10, 24 28 Oct. 2004. [ 4] P. Melin and O. Castillo, A new neuro fuzzy fractal approach for adaptive model based control of non linear dynamic systems: The case of controlling aircraft dynamics, in 1999 IEEE International FUzzy Systems Conference Proceedings, August 1999. [ 5] M. Steinberg and A. Page, High fidelity simulation testing of intelligent and adaptive aircraft control laws, in Proceedings of the American Control Conference, pp. 3264 3268, 8 10 May 2002. [ 6] B. Kim and A. Calise, Nonlinear flight control using neural networks, IEEE Transactions on Control Systems Technology, vol. 20, no. 1, Jan. Feb. 1997. [ 7] M. Bodson and J. Groszkiewicz, Multivariable adaptive algorithm for reconfigurable flight control, IEEE Transactions on Control Systems Technology, vol. 5, no. 2, Mar. 1997. [ 8] D. Shore and M. Bodson, Flight testing of a reconfigurable control system on an unmanned aircraft, in Proceedings of the 2004 American Control Conference, pp. 3747 3752, June July 2004. [ 9] B. Porter and C. Boddy, Design of adaptive digital controllers incorporating dynamic pole assignment compensators for high performance aircraft, in Proceedings of the IEEE 1989 Aerospace and Electronics Conference, pp. 372 379, 22 26 May 1989. [ 10] G. Tao, S. Chen, J. Fei, and S. Joshi, An adaptive actuator failure compensation scheme for controlling a morphing aircraft model, in Proceedings of the 42nd IEEE Conference on Decision and Control, pp. 4926 4931, Dec. 2003. [ 11] S. Chen, G. Tao, and S. Joshi, An adaptive actuator failure compensation controller for mimo systems, in Proceedings of the 41st IEEE Conference on Decision and Control, pp. 4778 4783, December 2002. [ 12] J. Farrell, M. Polycarpou, and M. Sharma, Adaptive backstepping with magnitude, rate, and bandwidth constraints: Aircraft lonitude control, in Proceedings of the American Control Conference, pp. 3898 3904, June 2003. [ 13] K. Astr om and B. Wittenmark, Adaptive Control. Addison Wesley Publishing Company, 1989. 73 A Appendix: PNI Micromag Specifications 74 75 76
Click tabs to swap between content that is broken into logical sections.
Rating  
Permanent URL  http://hdl.lib.byu.edu/1877/65 
Title  Navigation and Control Technologies for Autonomous Micro Vehicles 
BYU Creator 
Beard, Randal McLain, Timothy W. 
Keywords  UAV's; Micro air vehicles; Adaptive control; Magnetometers; Attitude estimation; Optic flow sensor 
Description/Abstract  Sponsorship: Air Force Research Laboratory / Munitions Directorate. During this project we focused on four primary objectives which are listed below. 1. Magnetometer Integration. Integrate and flight test magnetometers with the current version of the autopilot (Spiral 1). 2. HeightAboveGround Sensor Integration. Integrate and flight test heightaboveground sensors with the current version of the autopilot (Spiral 1). 3. Automatic Gain Adjustment. Develop and flight test automatic gain adjustment algorithms that automatically tune the servo loops of the autopilot. 4. Automatic Trim Seeking. Develop and flight test automatic trim seeking algorithms to recursively estimate the trim values of the UAV. 
Source  BYU Magic Lab; 
Date Original  20050825 
Language 
English eng en 
Publication Type 
Technical Reports 
Type 
Text 
Format  application/pdf 
Owning Institution 
Brigham Young University 
BYU College 
Ira A. Fulton College of Engineering and Technology 
BYU Department 
Electrical and Computer Engineering 
ScholarsArchive Collection 
ElectricalandComputerEngineering EngineeringandTechnology 
Copyright Status  (c) 2005 Randal Beard and Timothy McLain; 
Copyright Use Information  http://lib.byu.edu/about/copyright/generic.php 
Full Text  Navigation and Control Technologies for Autonomous Micro Vehicles Sponsorship: Air Force Research Laboratory / Munitions Directorate. During this project we focused on four primary objectives which are listed below. 1. Magnetometer Integration. Integrate and flight test magnetometers with the current version of the autopilot (Spiral 1). 2. HeightAboveGround Sensor Integration. Integrate and flight test heightaboveground sensors with the current version of the autopilot (Spiral 1). 3. Automatic Gain Adjustment. Develop and flight test automatic gain adjustment algorithms that automatically tune the servo loops of the autopilot. 4. Automatic Trim Seeking. Develop and flight test automatic trim seeking algorithms to recursively estimate the trim values of the UAV. Final Report Award No. F08630 03 1 0017 Navigation and Control Technologies for Autonomous Micro Vehicles August 25, 2005 Principal Investigator/ Technical Point of Contact Randal W. Beard, Associate Professor Department of Electrical and Computer Engineering Brigham Young University Provo, Utah 84604 USA voice: ( 801) 422 8392 fax: ( 801) 422 0201 beard@ ee. byu. edu Co Principal Investigator Timothy W. McLain, Associate Professor Department of Mechanical Engineering Brigham Young University Provo, Utah 84604 USA mclain@ byu. edu AFRL/ MN Point of Contact Carrie Fowler AFRL/ MNAV 101 West Eglin Blvd, Suite 346B Eglin AFB, FL 32542 ( 850) 882 8876 x3383 carrie. fowler@ eglin. af. mil CONTENTS Contents Cover Page i Table of Contents iii 1 Objectives 1 2 Magnetometer Integration 1 2.1 Magnetometer Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2.2 Motor Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3.1 Ground Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3.2 Airborne Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4 Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4.1 2 Axis Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4.2 3 Axis Heading Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5 2 Axis Tilt Compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.6 Magnetometer Wind Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Height Above Ground Sensor 9 3.1 Landing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Orbit Routine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.1 Heading Field Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.2 Correcting Steady State Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2.3 Glide Slope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 Barometric Altimeter Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4 Ultrasonic Range Finder Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Optic Flow Based Landing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Automatic Trim and Gain Tuning 23 4.1 UAV Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 Scheme A: nonlinear MRAC using backstepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3 Scheme B: nonlinear MRAC without backstepping . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Scheme C: nonlinear MRAC grouping terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.4.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.4.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.5 Scheme D: nonlinear MRAC simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5.1 Pitch Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5.2 Roll Attitude Hold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 ii CONTENTS 4.6 Flight Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.6.1 Platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.6.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5 Students Funded Under Project 73 A Appendix: PNI Micromag Specifications 74 iii 1 Objectives During this project we focused on four primary objectives which are listed below. 1. Magnetometer Integration. Integrate and flight test magnetometers with the current version of the autopilot ( Spiral 1). 2. Height Above Ground Sensor Integration. Integrate and flight test height above ground sensors with the current version of the autopilot ( Spiral 1). 3. Automatic Gain Adjustment. Develop and flight test automatic gain adjustment algorithms that automatically tune the servo loops of the autopilot. 4. Automatic Trim Seeking. Develop and flight test automatic trim seeking algorithms to recursively estimate the trim values of the UAV. The last two objectives were handled within one uniform framework and will therefore be described together. 2 Magnetometer Integration Adding a magnetic compassing device to the sensing suite of the Kestrel Autopilot is very desirable for several reasons. The BYU autopilot uses GPS for an estimate of ground track. With the use of a magnetometer, true heading ( ) can also be estimated. The use of true heading instead of ground track could significantly improve the accuracy of our gimbaled camera tracking system. Without true heading the gimbaled tracker will be in error by at least the value of the crab angle and the crab angle can be significant in high winds. In addition, magnetic heading allows for the computation of the wind vector which can easily be estimated in real time by using magnetometer heading, GPS velocity, airspeed, and ground track heading. 2.1 Magnetometer Selection When first selecting a magnetometer to interface with the autopilot several models were considered including the Honeywell HMC 1022 2 axis analog magnetometer and the PNI Micromag 2 axis digital magnetometer. Both sensors have a footprint of less than 0.6 x 0.5 inches. However the PNI Micromag magnetometer was selected because it incorporates built in temperature compensation and uses a SPI bus protocol that was straight forward to implement on the autopilot. The PNI Micromag also has selectable resolution and sample rate and the technical support at PNI is excellent. PNI developed a 3 axis version of the Micromag within several months which is currently being used on the autopilot. A table of Specifications and pin outs is included in Appendix A. Figure 1: PNI Micromag two axis digital magnetometer 1 2.2 Motor Noise 2.2 Motor Noise It was theorized that the magnetometer heading would be adversely influenced by induced magnetic fields from large brushless motors that are used for propulsion onboard Unmanned Air Vehicles ( UAVs). Fortunately, it was found that motor noise could be mitigated by twisting and shielding the high current carrying wiring. The alternating nature of brushless motors causes high frequency low amplitude noise and can also introduce a bias on the magnetometer readings. The high frequency component of the noise can be filtered out digitally. The biasing component of the motor noise is correlated with motor current. Work is currently being done to remove the effects of the electric motor on the magnetometer readings. The magnetometer is mounted on the autopilot which is located about 4 to 6 inches from the motor. If possible this distance should be increased to further reduce the effects of noise. All current carrying wiring should also be placed as far from the magnetometer as possible in addition to being twisted and/ or shielded. 2.3 Calibration In order to use magnetometers effectively a good calibration must be employed. This is easier said than done. First the magnetometer must be placed in the UAV in its completed form. Every time that something metallic is added to the UAV or every time the magnetometer is relocated on the UAV a new calibration should be performed. The magnetometer should be mounted as far as feasibly possible from current carrying wires, ferrous metals, and magnetic components. Then a magnetically clean environment should be located for calibration. A magnetically clean environment is one that is distant from steel and iron structures, power lines, and magnetic items. As a rule of thumb the calibration should be performed at least 5 times the diameter away from the nearest ferrous and magnetic objects. For example, if a 10 foot automobile is nearby, move at least 50 feet away from the automobile before calibrating the magnetometers. Beware that underground power lines and nearby power transformers could significantly affect the calibration. 2.3.1 Ground Calibration If it is intended to calibrate on the ground, the person performing the calibration should be sure to remove any metal objects and powered devices, such as keys and cell phones, from the area. Next, the ground station operator should set Calibrate Magnetometers to a value of 1. This starts logging data for the X and Y axes and stores the maximum and minimum values ( X m a x , X m i n , Y m a x , Y m i n ) for each axis. Then the person calibrating should hold the plane as level as possible while slowly turning in a circle in the horizontal plane through two complete rotations. Next the ground station operator should set Calibrate Magnetometers to a value of 2. This stops logging magnetometer data and calculates the scale factors ( X scale, Y scale) and offset values ( X offset, Y offset) for each axis. To ensure good calibration it is recommended that the X and Y magnetometer data be plotted on the x y plane and displayed for the operator. The X and Y magnetometer data should appear to be an ellipsoid. Applying the scale factors and offsets should circularize the ellipsoid and shift its center toward zero on the plot ( see Figure 2). For a 3 axis magnetometer the Z  axis should be calibrated separately. To do this the ground station operator should set Calibrate Magnetometers to 3. This starts searching for maximum and minimum Z axis values ( Z m a x , Z m i n ). The airplane must then be rotated about either the X or Y  axis for a full rotation. It is also valid to only turn the aircraft completely upside down, gathering the data points in both orientations. Once the rotation is completed Calibrate Magnetometers should be set to 4 to calculate the Z scale factor and offset values ( Z scale and Z offset respectively). The scale factors and offsets are calculated using the following equations: X offset = ( X m a x + X m i n ) = 2 ; Y offset = ( Y m a x + Y m i n ) = 2 ; Z offset = ( Z m a x + Z m i n ) = 2 ; ( 1) X scale = ( X m a x X m i n ) = 2 ; Y scale = ( Y m a x Y m i n ) = 2 ; Z scale = ( Z m a x Z m i n ) = 2 : ( 2) Once the calibration is complete be sure to save the calibration parameters in the autopilot flash memory and to a file. 2 2.4 Heading Estimation  2000  1500  1000  500 0 500 1000 1500 2000  2000  1500  1000  500 0 500 1000 1500 2000 X axis counts Y axis counts raw data calibrated data Figure 2: Y axis data plotted against X axis data for a typical two axis calibration. The uncalibrated data is in red while the data with the applied scaling and offset is blue. 2.3.2 Airborne Calibration The most effect method for calibrating the magnetometer was found to be an airborne calibration. An airborne calibration places the UAV in an environment that is far from magnetic disturbances on the ground. To calibrate the magnetometers in the air the first step is to set up an orbit point with a large orbit radius ( between 200 and 300 meters). Next, get the UAV up in the air and flying in its commanded orbit. Then start the magnetometer X Y calibration by setting Calibrate Magnetometers to 1. This starts logging data for the X and Y axes and stores the maximum and minimum values ( X m a x , X m i n , Y m a x , Y m i n ) for each axis. Once the plane has completed at least one complete orbit, set Calibrate Magnetometer to 2 which stops logging magnetometer data and calculates the scale factors ( X scale, Y scale) and offset values ( X offset, Y offset) for each axis. To calibrate the Z  axis, fly on an east west flight path. Then set Calibrate Magnetometers to 3. Next, invert the UAV momentarily before returning to normal flight. Barrel rolls could also be performed to achieve the same goal. Finally, set Calibrate Magnetometer to 4. This completes the Z  axis calibration and calculates the Z scale factor and offset values ( Z scale and Z offset respectively). As an alternative method for calibrating the Z  axis, set up a north south flight path. Then set Calibrate Magnetometers to 3. Next perform a complete loop or split S maneuver. Then set Calibrate Magnetometers 4. In either case the goal of the calibration routine is to capture the maximum and minimum values of all three axes. Be sure to save the calibration parameters in the autopilot flash memory and to a file. Once the rotation is completed the scale factors and offsets are calculated using Equations ( 1). 2.4 Heading Estimation 2.4.1 2 Axis Heading Estimation To calculate heading, raw magnetometer readings ( B x and B y ) are adjusted using the scale factors and offsets calculated during calibration using the following equations: B 0 y = B y Y 0 Y g B 0 x = B x X 0 X g : 3 2.4 Heading Estimation Figure 3: 3 axis heading ( red) compared to 2 axis heading ( green) along with corresponding roll ( ) and pitch ( ) angles taken during a clockwise circular orbit. To estimate heading the following formula is used: = atan2 ( B 0 y ; B 0 x ) + declination The heading estimate on a 2 axis magnetometer is only valid when the magnetometer is oriented level. A 2 axis magnetometer will be in error by as much as 20 degrees with roll and pitch angles of up to 15 degrees ( see Figure 3). 2.4.2 3 Axis Heading Estimation To calculate heading, raw magnetometer readings ( B x , B y , and B z ) are adjusted using the scale factors and offsets calculated during calibration, where B 0 x , B 0 y , and B 0 z are the now calibrated x , y , and z components of the magnetic field, i. e., B 0 x = B x X 0 X g B 0 y = B y Y 0 Y g B 0 z = B z Z 0 Z g : Then the adjusted magnetometer readings are un rolled and un pitched using standard rotation matrices as follows: 0 @ X h Y h Z h 1 A = 0 @ c o s 0 s i n 0 1 0 s i n 0 c o s 1 A 0 @ 1 0 0 0 c o s s i n 0 s i n c o s 1 A 0 @ B 0 x B 0 y B 0 z 1 A : ( 3) The un pitched and un rolled horizontal components of the magnetic field then become: X h = B 0 x c o s B 0 y s i n s i n B 0 z c o s s i n Y h = B 0 y c o s B 0 z s i n : From these equations heading is computed as: = atan2 ( Y h ; X h ) + declination : 4 2.5 2 Axis Tilt Compensation Figure 4: 3 axis estimation of heading as the UAV flew in a clockwise orbital flight path. Corresponding roll and pitch angles are given in the lower plot This method been used successfully to estimated the actual heading of the UAV even while maneuvering at pitch and roll angles of over 35 degrees ( see figure 5). 2.5 2 Axis Tilt Compensation A 2 axis magnetometer can be pitch and roll compensated to achieve results comparable to those of a 3 axis magnetometer by estimated the z component of the magnetometer. This is done by using additional information that is known about the Earth s magnetic field such as inclination and magnitude of the Earth s field. Inclination is the angle between the Earth s surface and the total magnetic field vector. This angle, often referred to as the dip angle, mainly varies depending on latitude. Inclination is not to be confused with declination which is the angular difference between true north and magnetic north. The magnitude of the Earth s magnetic field is relatively constant, however, aperiodic variation of earth s magnetic field over time and place does occur. The Z  axis component can be estimated using the following expression: B 0 z = s i n B Total + s i n B 0 x c o s s i n B 0 y c o s c o s ; where is the inclination and B Total is is the magnitude of the Earth s magnetic field. Both of these values are only valid in the region where the magnetometer is to be used. In the case of the PNI magnetometer which outputs readings in terms of ADC counts rather than in Tesla or Gauss, The value for B Total must also be input in terms of ADC counts. After estimating the z component, the adjusted magnetometer readings are un rolled and un pitched using Equation ( 3) with B 0 z replaced by B 0 z . This method has been used successfully to estimate the actual heading of the UAV even while maneuvering at pitch and roll angles of over 35 degrees. In simulation this method of estimating the Z  axis is equally as good as using a 3 axis magnetometer ( see Figure 5). 2.6 MagnetometerWind Estimation Knowledge of true heading and ground track heading as well as airspeed and ground speed allows for the computation of the wind vector which can easily be estimated in real time. Using what is called the wind triangle, the wind speed and wind direction can both be found using simple trigonometry ( see figure 6). Using the wind triangle to subtract the air velocity vector from the ground velocity vector yields the wind vector. This method has accurately measured the wind speed and wind direction in flight on a UAV ( see figure 7). 5 2.6 Magnetometer Wind Estimation 0 10 20 30 40 50 60 70 80 90  200  100 0 100 200 Comparison Between Three magnetic heading estimation methods degrees 0 10 20 30 40 50 60 70 80 90  5 0 5 10 15 Time ( sec) degrees psi 3D Mag 2D Mag 2D Mag Est 3D Mag f q Constant roll of 10 Degrees Figure 5: A simulation comparison of all three methods discussed for estimating psi including using a 3 axis mag ( 3D Mag), 2 axis mag ( 2D Mag), and 2 axis tilt compensation by estimating the Z axis ( 2D Mag Est 3D Mag). Corresponding roll and pitch angles are given in the lower plot. Wind speed / Direction Airspeed / Heading ( y ) Ground speed / Ground track ( c ) Figure 6: Wind triangle showing that the addition of the air velocity vector and wind vectors should equal the ground velocity vector. Additional magnetometer data is shown in Figure 8 and 9. 6 2.6 Magnetometer Wind Estimation 0 20 40 60 80 100 120 300 320 340 Wind Estimation ( Actual Wind: from North West at 1 to 2 m/ s) Wind Direction ( degrees) 0 20 40 60 80 100 120 0 2 4 6 Wind Speed ( m/ s) 0 20 40 60 80 100 120 240 260 280 300 320 degrees sample number Estimated Reported Estimated Reported GPS Course Mag Heading Figure 7: Top: Estimated wind direction. Middle: Estimated wind speed. Bottom: GPS Course and Magnetometer heading. The actual wind speed and direction as reported on day of flight were 2 m/ s from the north west ( 315 degrees). 0 200 400 600 800 1000 1200 1400 1600 1800  4  2 0 2 4 x 10 4 x y z 0 200 400 600 800 1000 1200 1400 1600 1800  4  2 0 2 4 x 10 4 total field magnitude 0 200 400 600 800 1000 1200 1400 1600 1800  50 0 50 sample num degrees phi theta Figure 8: Top: Magnetometer flight data of X , Y , and Z axes. Middle: Total magnetic field vector calculated from X , Y and Z data. It is relatively constant which indicates that the magnetometer did not enter the presence of any significant magnetic disturbances. Bottom: Roll ( ) and Pitch ( ) data corresponding to the collected magnetometer data in flight. 7 2.6 Magnetometer Wind Estimation Figure 9: Heading estimates in a circular orbit using both 2 and 3 axis magnetometers. 8 3 Height Above Ground Sensor The usefulness of mini and micro unmanned air vehicles ( UAVs) depends largely on the ability to hand launch and recover the UAV in uncontrolled environments. Efforts to implement effective autonomous landing of mini UAVs in uncontrolled environments have been most hindered by lack of small, lightweight, effective sensors to determine height above ground ( HAG). Because of the lack of a suitable sensor for determining HAG, small UAVs have traditionally relied solely on the barometric altimeters with which they are already equipped for HAG estimation. This approach works well for controlled environments in which the altitude of the landing environment relative to the altitude of the takeoff environment is known, and the terrain is relatively smooth. However, in uncontrolled environments where this offset is not known, or where the terrain is not acceptably smooth, other methods must be utilized to effectively determine HAG. Two methods with acceptably small sensor payloads have been proposed as effective means of HAG estimation. These methods are: 1. Estimation of HAG via Ultrasonic Ranging 2. Estimation of HAG via Optic Flow It is proposed that these sensors can be combined with information from the barometric altimeter to reliably and accurately determine HAG. This accurate estimate of HAG could then be combined with robust control algorithms to effectively and safely land small UAVs in uncontrolled environments. This report presents the experimental results of efforts to estimate HAG and land small fixed wing UAVs in uncontrolled environments using each of the sensors outlined above. Because each of the sensors tested will couple with the barometric altimeter in determining HAG, the experimental results will be presented as follows: 1. Overview of the landing algorithm, 2. Landing results using only the barometric altimeter, 3. Landing results using an ultrasonic sensor combined with the barometric altimeter, 4. Landing results using an optic flow sensor combined with the barometric altimeter, 5. Comparison of results using each of the sensors. Each of these sections will address the specific sensor tested and the algorithms used to determine HAG from the data provided by the sensor. 3.1 Landing Algorithm 3.1.1 Overview The landing algorithm developed and tested for this project consists of two user defined waypoints: the approach point and the landing point ( Figure 10). LANDING POINT APPROACH POINT Figure 10: Landing Algorithm The approach point serves as a point about which the airplane may orbit until it reaches a specified altitude, at which point it will break out of its orbit and follow a glide slope to the landing point. The approach point is fully defined by the following parameters: 9 3.2 Orbit Routine 1. Airspeed The desired airspeed which will be maintained throughout the circle down portion of the landing. 2. Radius The radius of the orbit which will be maintained throughout the circle down portion of the landing 3. Start Altitude The altitude to which the plane will attempt to climb or descend while en route to the approach point. The true start altitude for the descent will be the altitude to which the plane has successfully climbed or descended by the time it is within 1.5 orbit radii of the landing point. 4. End Altitude  The altitude at which the airplane will break out of its orbit and begin its glide slope in to the landing point 5. Relative East Coordinate  The distance east of home ( in meters) around which the plane will orbit while descending from its start altitude to its end altitude. 6. Relative North Coordinate  The distance north of home ( in meters) around which the plane will orbit while descending from it start altitude to its end altitude. 7. Descent Rate  The rate ( in meters/ second) at which the plane will descend while orbiting the approach point. The landing point is fully defined by the following parameters: 1. Airspeed  The airspeed which will be maintained throughout the glide slope from approach point to landing point. 2. Flare Height  The height above ground at which the airplane will cut throttle and attempt to hold zero roll while maintaining its glide slope by commanding pitch to control around desired altitude. 3. Relative East Coordinate  The distance east of home ( in meters) at which the plane will touch down. 4. Relative North Coordinate  The distance north of home ( in meters) at which the plane will touch down. Altitude control throughout the landing routine is accomplished through activation of various PID loops. For the portion of the landing sequence in which the plane is en route from its last waypoint to the approach point altitude is maintained as follows: 1. Actual altitude is within a preset distance ( usually around 15 m) of desired altitude  Pitch from Altitude, Throttle from Airspeed. 2. Actual altitude is greater than the same preset distance from desired altitude and is less than desired altitude  Pitch from Airspeed, Full Throttle. 3. Actual altitude is greater than the same preset distance from desired altitude and is greater than desired altitude  Pitch from Airspeed, Zero Throttle. For the portion of the landing sequence in which the plane is circling down to its end altitude and then following its glide slope to the landing point altitude is maintained using the pitch from altitude and throttle from airspeed loops. Throughout the entire landing sequence, until the flare at which desired roll is held at zero, heading is controlled using the roll from heading PID loop. Roll rate and pitch rate are also controlled to provide damping. Desired headings are commanded based on the generation of a heading field which forms a limit cycle around the desired orbit radius. A full description of this algorithm is provided in the following section. 3.2 Orbit Routine 3.2.1 Heading Field Generation The ability to implement the landing routine described in the previous section required the development and implementation of a control algorithm which would allow the airplane to fly accurate orbits, at minimal orbit radii, even in high wind conditions. The orbit algorithm implemented generates a desired GPS heading based on the position of the airplane such that the corresponding heading field forms a limit cycle around the desired orbit radius ( Figure 11). Such a heading field is generated using Equations ( 4). 10 3.2 Orbit Routine Figure 11: Typical Heading Field for a 50 Meter Radius Orbit ( Desired Orbit Radius Indicated by Solid Circle). d = 8 > < > : center t a n 1 r c d ; d 2 r c center r ; r c d 2 r c center + + + r ; d < r c ( 4) In the above equations d is the distance from the center of the orbit, d is the desired course, center is the compass heading to the center of the orbit, r c is the desired orbit radius, is the included angle between d and center for d = 2 r c , is the complimentary angle to , and is a gain corresponding to the field strength of the orbit. Increasing makes the field more reactive around the desired orbit radius and decreasing makes the field less reactive as shown in Figure 12. r is a parameterization of the distance from the center of the orbit based on r c . The formulation of r is such that 0 r 1 and r = 1 when d = r c . The value of r is calculated according to Equations ( 5). r = ( 1 d r c r c ; r c d 2 r c d r c ; d < r c ( 5) Figure 12: High Gain Heading Field ( Left) and Low Gain Heading Field ( Right) for a 50 Meter Radius Orbit ( Desired Orbit Radius Indicated by Solid Circle). The desired heading from the heading field described by equations ( 4) is used to generate a roll command for the airplane according to d = k ( d a ) ( 6) 11 3.2 Orbit Routine 3.2.2 Correcting Steady State Error Because the heading control represents a type one system and flying an orbit based on desired heading is essentially tracking a ramp in heading, the actual course of the airplane tracks the desired course with a steady state error. This steady state error in course represents a lag between desired and actual course which translates to actual orbits that are larger than the commanded orbit. The magnitude of the steady state error in orbit radius is a function of the desired radius of the orbit, the groundspeed of the airplane, and the heading field gain . Two different approaches were used to overcome this effect. The first approach involved adapting the radius command to remove the steady state error. The second approach, which was more effective in high winds, involved calculating a nominal roll angle based on the desired heading rate and the coordinated turn equation. Adaptive Radius Commands Because the actual orbit radius of the airplane is always exterior to the commanded radius, for every dynamically feasible orbit there exists some smaller radius command that will result in the actual desired orbit. The adaptive radius command progressively increases or decreases an offset between the desired orbit radius entered by the operator and the desired orbit radius passed to the orbit algorithm. This offset is continuously incremented or decremented until the plane is orbiting at the radius requested by the operator. Orbit results with and without the adaptive radius command are shown below in Figure 13. Telemetry from actual flight tests of the adaptive radius method are shown in Figure 14. Figure 13: Circular Orbits ( From Simulation) with and without an adaptive radius command. Nominal Roll Angle Based on Desired Heading Rate The adaptive radius command approach works extremely well for situations where wind speed does not exceed 50% of airspeed. However, high winds cause the groundspeed of the airplane to vary significantly based on the course of the airplane relative to the wind direction. This variation in groundspeed causes a corresponding variation in the radius offset necessary to maintain the desired orbit. A better approach for high wind conditions was developed which involved calculating a nominal roll angle that resulted in the heading rate necessary to fly the desired orbit. The heading field was then used to give roll commands about this nominal roll value. Calculating the nominal roll angle first required taking the derivative of the heading field with respect to time. This resulted in the formation of a heading rate field over the same space as the heading field. The equations for the heading rate field are given by _ d = V r c r 2 c + 2 d 2 d ( r 2 c + d 2 ) 3 2 ( 7) = V d s i n ( + r ) + r ( 1 ) V c o s ( + r ) r c ( 8) = V d ( s i n ( + r ) r c o s ( + r ) ) ( 9) An example of a heading rate field for a given heading field strength and constant groundspeed is shown in Figure 15. Equations ( 4) combined with Equations ( 7) ( 9) allow the calculation of both a desired heading and a desired heading rate for any location relative to the center of the orbit. The desired heading rate can then be converted to a 12 3.2 Orbit Routine Figure 14: Actual Flight Results for 50, 70, 100, and 150 m Radii Orbits Flown in 2 3 m/ s Wind . nominal roll angle using the standard coordinated turn equation nom = _ g V g ( 10) and the control law described in Equation ( 6) can be adapted as nom = _ g V g ( 11) d = nom + k ( d a ) ( 12) The control law shown in Equation ( 12) is more effective at tracking orbits in high winds, but somewhat less steady at tracking orbits in low wind conditions than is the adaptive radius command method. A comparison of the two methods flying a 100 meter radius orbit in simulation with a constant wind speed of 8 m/ s is shown in Figure 16. 3.2.3 Glide Slope The glide slope implemented in this landing routine was a simple linear glide slope calculated by decreasing desired altitude as a function of the distance from which the glide slope began. This formulation allows the desired altitude to be driven negative for distances from the initiation of the glide slope which are farther than the desired landing point. This makes the glide slope robust with respect to negative offsets between measured HAG and actual HAG. Expected overshoot or undershoot of the landing point can approximated as Expected Over / Undershoot 1 t a n a + b ; ( 13) where is the glide slope, a is the altitude control error ( difference between desired and actual HAG), and b is the measurement error. The relationship shown in Equation ( 13) indicates that larger glide slopes will tend to reduce the overall overshoot or undershoot error for the landing provided the plane is capable of matching the commanded glide 13 3.3 Barometric Altimeter Based Landing Figure 15: Typical Heading Rate Field. slope. However, larger glide slopes also cause increased impact velocities at touch down. Testing showed a glide slope of anywhere from 7 to 12 degrees to provide a good balance between accuracy and impact at touch down for the MAGICC Lab airplanes. 3.3 Barometric Altimeter Based Landing 3.3.1 Theory The Kestrel autopilot comes equipped with a barometric pressure sensor as part of its standard sensor suite. Before a given flight, the barometric pressure sensor is zeroed such that it references the barometric pressure at the takeoff location as the pressure corresponding to zero altitude. Pressure decreases with altitude according to P = P ref Z h h ref g h d h : ( 14) If the density of air ( ) is assumed to be sufficiently constant over the range of altitudes being measured, then altitude above the takeoff point can be determined by h = P ref P actual g = ( P ref P actual ) ; ( 15) where the constant of proportionality ( ) can be determined by taking readings of the pressure sensor at a number of different known reference altitudes and fitting a straight line to the data set. Once the constant of proportionality has been determined, readings of the pressure sensor from the analog to digital converter can be converted directly into relative altitudes above the takeoff point. This sort of calibration allows the barometric pressure sensor to provide valid and accurate estimates of HAG where the terrain is sufficiently flat. Where this is not true, barometric pressure can still provide a rough estimate of HAG until other sensors are able to collect valid readings of true HAG. After another sensor has reported a valid reading of true HAG variations in barometric pressure from the pressure at which true HAG was determined can be used to further estimate the true HAG. 3.3.2 Application Temperature Compensation The pressures reported by the pressure sensor on the Kestrel autopilot proved to be somewhat sensitive to the temperature of the sensor. Particularly this would become a problem when the autopilot was allowed to heat up as it sat static during calibration and programming, and then cooled as airflow over the airfoil dissipated heat during flight. This calibration at high temperature and subsequent cooling would cause the autopilot to consistently read approximately 20 meters high. The increase in pressure readings with decreased sensor temperature was empirically determined to be linear over the range of operating temperatures. The constant of proportionality was likewise empirically determined and implemented as P actual = P measured + ( T measured T ref ) ( 16) 14 3.3 Barometric Altimeter Based Landing Figure 16: A Comparison of Nominal Roll Angle and Adaptive Radius Based Methods ( Simulation). This temperature compensation makes the pressure readings from the autopilot robust to the variations in temperature that are encountered in standard autopilot operation. Compensation for Airflow Over the Airfoil Another source of error in the readings from the autopilot s pressure sensor stems from the fact that the MACICC Lab airplanes are not normally equipped with static pressure tubes. While in flight a considerable offset in the absolute pressure measured by the pressure sensor is introduced by airflow over the airfoil in which the autopilot is housed. The increased speed of the airflow over the airfoil results in a Bernoulli effect that causes the pressure measured along the airfoil to be less than the true absolute pressure. This low pressure region in which the actual pressure measurements are taken causes the altitude measurement computed by the autopilot to consistently read high. The relationship between the velocity of the airflow over the airfoil and the induced offset in the pressure reading was determined empirically through wind tunnel testing to be a second order function of airspeed. This functional relationship can be combined with the temperature correction term in Equation ( 16) to calculate pressure as P actual = P measured + C 1 v 2 diffpres + C 2 v diffpres + ( T measured T ref ) ( 17) where v diffpres is the airspeed past the airplane as measured by the airplane s pitot probe attached to a differential pressure port. Once the constants ( C 1 , C 2 , ) have been empirically determined, the offset can be computed each time through the loop, and the barometric altitude reading can be adjusted accordingly. Testing has shown the constants determined for one airplane to be valid for other airplanes of a similar type regardless of slight variations in the airfoil and in autopilot location. This allows the constants for a fleet of similar airplanes to be determined through a rigorous wind tunnel calibration of a single airplane. Further testing is required to determine how robust these constants are with respect to major changes in airframe or autopilot location. This method of compensating for airflow over the airfoil allows a large number of similar airplanes to be flown, using the exact same constants, without each being equipped with static tubes which can add unnecessary weight and complication. 3.3.3 Results Using the calibration described above and testing over relatively short flights, the MAGICC Lab airplanes consistently read within 1  1.5 meters of the actual altitude relative to home. This allows relatively accurate landings where the altitude of the landing point can be safely assumed to be the same as the altitude of home. This works especially well when trying to land the airplane at or near its home location which is most often the case. 15 3.4 Ultrasonic Range Finder Based Landing If the actual altitude of the airplane matched the desired altitude perfectly, then overshoot or undershoot would be solely a function of error in HAG estimation. For glide slopes of 7 to 12 degrees and a measurement error of 1  1.5 meters, a corresponding overshoot or undershoot error of around 4  12 meters can be expected ( see Equation ( 13)). Actual results for 26 consecutive landing runs are shown below in Figure 17. The first 22 autoland runs were each from the main approach point. These were followed by one run apiece from each of the four alternate approach points. Figure 17: Results of Twenty six Consecutive Autolands from Five Different Approach Points. An average overshoot or undershoot error of 4  12 meters was estimated. Actual results from these 26 landing show an average over/ undershoot error of 7.6 meters with a standard deviation of 5.4 meters and a median value of 6 meters. 3.4 Ultrasonic Range Finder Based Landing 3.4.1 Theory Ultrasonic sonar transducers determine distance from an object by measuring the time of flight of an ultrasonic signal which is emitted by the transducer and then reflected by the object. The speed of sound in an ideal gas is given by c = p R T k ; ( 18) where is the adiabatic constant ( 1.402 for air), R is the universal gas constant ( 297.05 kJ/( kg K) for air), and T is the absolute temperature in Kelvin. Linearizing Equation ( 18) around 0 degrees C gives a convenient approximation for the speed of sound as c air = 3 3 1 : 5 + 0 : 6 T c : ( 19) 3.4.2 Application Sensor The MUST01 Ultrasonic Sonar Transducer manufactured by Mekatronix was selected and used as the ultrasonic ranging module and transducer for these tests. This sensor provided both the ranging board and transducer in one integrated package as well as allowing a sensing range of up to 12.2 meters with a resolution of 4.8 cm. The entire unit has an overall diameter of 4.3 centimeters and a total mass of 17 grams. The sensor emits a 40  50 kHz sound wave and has a nominal beam angle of 15 degrees. The MUST01 unit emits a signal 10 times a second and returns an analog voltage proportional to each signal s time of flight. This value is alpha filtered on the ranging board according to V out = 0 : 9 V new + 0 : 1 V t e x t l a s t : ( 20) 16 3.4 Ultrasonic Range Finder Based Landing Figure 18: MUST01 Ultrasonic Sonar Transducer. The voltages from the ultrasonic sensor are read into the autopilot via the autopilot s analog to digital converter. The ultrasonic sensor was calibrated by holding the sensor at a number of fixed distances from a surface with good acoustical properties and then recording the analog to digital converter values corresponding to that distance. The A/ D values were then related directly to their corresponding distances through a constant of proportionality. This constant was determined through a least squares linear fit of the data. Compensation for Airplane Attitude Because the ultrasonic sensor has a nominal beam angle of 15 degrees and reports the closest item from which it receives a valid return signal, the reading from the ultrasonic sensor may not always reflect HAG depending on the attitude of the airplane. If the airplane is either rolled or pitched sharply the signal emitted by the ultrasonic sensor will reflect off the ground at some distance significantly displaced from the point located directly beneath the airplane. The distance from the airplane to this displaced point will be significantly greater than the distance from the airplane to the point directly beneath it. Thus this reading will be significantly greater than the airplane s actual HAG. However, if the attitude of the airplane is known, and the earth can be assumed to be sufficiently flat over the span from the point directly below the airplane to the point at which the signal is reflected, then the actual HAG can still be determined by h actual = h measured c o s ( j j ) c o s ( j j ) + d s i n ; ( 21) where is one half the nominal beam angle, and refer to pitch angle and roll angle respectively, and d is the lateral distance from the center of the airplane at which the sensor is mounted. Equation ( 21) is only valid, and only necessary, when either , , or both are greater than . It is important to note that when this equation must be used to determine HAG, the maximum distance at which true HAG can be accurately determined becomes a dynamic function of airplane attitude rather than a constant ( 12.2 meters for the sensor tested). The airplane attitude dependant value of maximum distance at which true HAG can be accurately determined can be found setting h measured equal to the maximum distance at which true HAG can be determined for level flight, and then solving Equation ( 21) for h actual. Another important consideration when using an ultrasonic sensor to determine HAG is that there will be an angle of incidence associated with every surface at which the ultrasonic signal will not be reflected back toward the sensor. Equation ( 21) is not valid if the angle of incidence is greater than this angle. Instead there is a discontinuity at which the maximum value of true HAG that can be determined goes immediately to zero because no return signal is received. Flight tests have shown that for all outdoor environments tested this angle of incidence is acceptably small such that this discontinuity is not encountered under normal flight conditions. Limitations Due to Propeller Noise A major limitation of the ultrasonic sensor is its tendency to pick up false readings due to noise from the propeller. For the MAGICC Lab airplanes with the motor running at full throttle to near full throttle any reading over 8 meters would be rejected in favor of noise from the propeller. This makes it impossible to trust readings from the ultrasonic sensor at a range greater than 8 meters. For our testing we set the cutoff at which we would trust readings from the ultrasonic sensor at 7.9 meters. This severely limits the usefulness 17 3.5 Optic Flow Based Landing of the ultrasonic sensor by severely limiting its maximum attainable range. This also introduces uncertainty as to whether the sensor will be more or less susceptible to noise from a different combination of motor and propeller. The sensor s maximum range then becomes very platform dependant rather than strictly sensor dependant. The magnitude of the noise which the sensor picks up has been minimized by placing the sensor on the far wingtip of the airplane, as far from the motor and propeller as possible. This, however, introduces a significant payload to the extremity of the airplane causing balancing issues. Also, if the sensor were to be used on an airplane where it could not be mounted a significant distance from the propeller, a significant decrease in maximum range should be expected. Environmental Factors In order for the ultrasonic sensor to determine distances it must transmit and then receive back a signal of acceptable magnitude. For the case of the airplane this signal must be of acceptable magnitude such that it will be accepted by the sensor as the significant signal over noise from the propeller. Testing has shown that this is extremely reliable for flights over surfaces that are sufficiently acoustically reflective such as streets or well trimmed yards. However, for flights over more dense foliage ( fields with tall grass for example) the signal can be sufficiently muffled such that the sensor does not receive an acceptably large return signal until it is 1  2 meters off the ground. Thus maximum attainable range for the ultrasonic sensor is not only dependant on the particular airframe it is flown in ( especially sensor placement on that airframe) and that airplane s attitude, but it is also highly dependant on the acoustical properties of its flight environment. 3.4.3 Results In order for the ultrasonic sensor to facilitate autonomous landing in truly uncontrolled environments where landing altitude relative to takeoff altitude is not known, it must allow the airplane to circle down at a steady rate until true HAG is detected, and then break out of its orbit and follow a glide slope to the landing point. The ultrasonic sensor was not able to reliably perform this task. The sensor s extremely limited range coupled with loss of ultrasonic readings over dense foliage and loss of ultrasonic readings due to airplane attitude were the main contributing factors to the sensor s unreliability. 3.5 Optic Flow Based Landing 3.5.1 Theory Optic flow sensors can be used to determine HAG by relating the flow of features across an imaging array to the speed the imaging array is moving past the surface it is imaging and its distance from that surface. The number of pixels a given object moves in the imaging plane can be combined with data from the airplane s IMU and GPS to determine HAG according to [ 1] h = x 2 t a n fov 2 p n _ t s 2 c o s c o s ; ( 22) where h represents HAG, x represents the displacement of the sensor parallel to the plane being imaged over the sample time t s , represents the average number of pixels of displacement of features across the imaging plane over the same sample time, fov represents the field of view of the sensor, p n is the number of pixels in the imaging array in the direction of motion of the sensor, _ is the average pitch rate of the sensor over the sample period, is the average pitch of the sensor over the sample period, and is the average roll of the sensor over the sample period. 3.5.2 Application Sensor Testing of autonomous landing using optic flow to determine HAG was performed using Agilent s ADNS 2610 optic flow sensor. The ADNS 2610 has a small form factor, measuring only 10 x 12.5 mm and runs at 1500 frames per second. It requires a light intensity of at least 80 mW/ m2 at a wavelength of 639 nm or 100 mW/ m2 at a wavelength of 875 nm. The ADNS 2610 measures the flow of features across an 18 x 18 CMOS imager. It outputs two values, d x and d y , representing the total optic flow across the sensor s field of view in both the x and y directions respectively. These values are stored in a buffer which can store accumulated optic flow up to 128 pixels. Once the buffer is full no more optic flow readings are accumulated until the buffer has been read. Reading the buffer clears it and allows it to begin reaccumulating optic flow measurements. The flow data in the camera y direction corresponds to lateral motion of the airplane and for purposes of establishing HAG can be ignored. The flow data in the camera 18 3.5 Optic Flow Based Landing x direction can be combined with data from the IMU, GPS, and the field of view of the sensor to determine HAG according to Equation ( 22). Optics A critical consideration in outfitting an airplane with an optic flow sensor is the setup of the optics. Narrow angle lenses increase the distance at which the airplane will be able to pick up appreciable optic flow. This increases the distance at which the sensor will be able to accurately measure HAG for a given groundspeed. This comes at the expense of increasing the longitudinal size of the sensor. Narrow angle lenses can also decrease the low end range of the sensor by causing overflow in the optic flow sensor s dx register or by causing the angular flow rate seen by the sensor to become too fast for the sensor to register. Wide angle lenses allow for smaller longitudinal sensor size, but require that larger features be available in the environment. They also decrease the distance at which appreciable optic flow will be recorded for a given groundspeed. A particular advantage of wider angle lenses is that they are less susceptible to noise introduced by angular oscillations of the airplane in pitch and roll. In addition to focal length, lens diameter is an important consideration. This is especially important with regard to the corresponding f stop value of the lens. The f stop is defined as the ratio of a lens focal length to its diameter and is a measure of the amount of light that the optics allow to pass through to the imaging array. Lenses with larger f stops allow less light to pass through. Keeping in mind that the imaging array requires a light intensity of 80  100 W/ m2 in the correct spectral range it is important to select a lens with a low f stop in order to increase the sensor s operational range. Because the sensor automatically adjusts shutter speed to keep the light intensity at an ideal level, there is no disadvantage to selecting a lens with a lower f stop other than increased size. Another important implication of the f stop is that selecting a lens with an increased focal length ( a narrower angle lens) will not only increase the longitudinal size of the sensor, but also increase the lateral size of the sensor if the f stop is to be maintained. For the purposes of these tests lenses were selected which yielded fields of view of 6.5, 2.5, and 1.2 degrees when mounted on the sensor. These lenses had accompanying f stops of 2.0, 2.0, and 2.5 respectively. Each of the configurations is shown in Figure 19. Left FOV: 6.5o Long. Size: 25 mm Lat. Size: 30 mm Mass: 15 g F stop: 2.0 Center FOV: 2.6o Long. Size: 35 mm Lat. Size: 30 mm Mass: 23 g F stop: 2.0 Right FOV: 1.2o Long. Size: 50 mm Lat. Size: 30 mm Mass: 23 g F stop: 2.5 Figure 19: 1.2, 2.5, and 6.5 degree Field of View Optics Configurations. Comparable results were achieved using both the 2.5 degree and 1.2 degree field of view setups. However, the 1.2 degree field of view lens offered the slightly greater range while performing better over feature poor surfaces such as asphalt. The one disadvantage of the 1.2 degree field of view setup up was its increased susceptibility to noise caused by pitch and roll oscillations. Sample Rate Another critical consideration when equipping an airplane with an optic flow sensor is the selection of sample rate. The resolution of the sensor ( measured in units of distance per sensor count) is given by r = l i m h ! 0 0 @ h fov 2 p n t a n 1 V t s 2 ( h h ) t a n 1 V t s 2 h 1 A : ( 23) 19 3.5 Optic Flow Based Landing For a given sensor with a fixed field of view the sensor s resolution is a function of sample rate, groundspeed, and HAG. Plots of sensor resolution for a fixed field of view and varying sample rates and groundspeeds are shown in Figure 20. Figure 20: Plots of Resolution as a function of HAG and Groundspeed and HAG and Sample Rate. For high values of r each sensor counts corresponds to a relatively large distance measure. Correspondingly the noise on the sensor measured in sensor counts is amplified by r into higher magnitude noise in the HAG measurement. This makes large values r undesirable because they can significantly decrease the sensor s signal to noise ratio. It is possible to decrease r by decreasing the sample rate as shown in Figure 20. This, however, can be undesirable because it decreases the sensor s response time and can lead to overflow in the sensor s d x register. Furthermore, for a fixed sample rate, maintaining a low value of r for relatively large values of HAG results in unnecessarily low values of r for lower values of HAG. This means unnecessarily low sample rates are being maintained where significantly faster sample rates could be used without significant loss in signal to noise ratio. Such an approach also results in premature sensor overflow and corresponding loss of low end range. To solve these problems two approaches were developed using a dynamic rather than static sample rate. The first approach maintains an r value along a specified resolution curve regardless of groundspeed of the plane. The second approach uses the current velocity estimate and HAG estimate to maintain a constant value of r regardless of groundspeed and HAG. Constant Curve Resolution In general, flight control systems for UAVs utilize a feedback loop to control around a commanded airspeed rather than groundspeed. Even when the feedback loop is closed around groundspeed, the commanded groundspeed may or may not be realizable due to wind conditions. In either case in the presence of wind it is reasonable to assume that groundspeed may vary significantly depending on the course of the airplane relative to the wind direction. As shown in Figure 20 variations in groundspeed can alter the resolution curve significantly. However, if the sample rate is set according to t s = V ; ( 24) where represents the sample rate divisor, the value of r will follow a fixed resolution curve as a univariate function of HAG. This is useful because it eliminates the variations in resolution due to velocity and provides reliable repeatable resolutions regardless of wind conditions or commanded airspeeds. Constant Value Resolution When an estimate of the current HAG is available, the correct sample rate can be calculated such that the sensor resolution can be maintained constant regardless of HAG or groundspeed. This sample rate is found by solving Equation ( 23) for t s . The result is t s = l i m h ! 0 0 B B @ 1 V 2 ( h h ) V 2 h 1 t a n h fov 2 p n r V 2 4 h ( h h ) 1 C C A ; ( 25) where r now represents a user defined parameter that defines the resolution of the sensor for all combinations of HAG and groundspeed. In practice, the actual value of t s is capped at both the high and low end, but allowed to change in order to maintain a fixed resolution inside of some saturation block. In flight tests using a dynamic sample rate to fix sensor resolution cut noise drastically for high values of HAG while significantly increasing the sample rate and thus 20 3.5 Optic Flow Based Landing improving the response time for low values of HAG. In practice Equation ( 25) is implemented by dropping the limit expression and fixing h to some reasonably small number. Flight test results using both constant curve resolution and constant value resolution are shown in Figure 21. The spike in the optic flow measurement at sample 400 in the constant value resolution plot and 275 in the constant curve resolution plot corresponds to flight over the asphalt road. In practice such spikes are ignored by ignoring values for which the optic flow sensors reports poor surface quality. Figure 21: Optic HAG Results from Glideslope Descents Using Fixed Value Resolution ( Left) and Fixed Curve Resolution ( Right). 3.5.3 Results Ranging The accuracy of the distance measurements calculated from optic flow was tested by comparison with a laser rangefinder. Testing was performed by mounting the optic flow sensor and laser rangefinder perpendicular to the motion of a vehicle driving along a freeway at varying distances from the freeway sound barrier wall. Range values were recorded for both the optic flow calculations and the laser rangefinder. The results of this test are shown in Figure 22. The results show very close agreement and validate the readings obtained from the optic flow calculations. Figure 22: Range Data from Optic Flow Calculations Compared to Laser Rangefinder Data. Range data from the optic flow sensor in flight tests showed reliable, repeatable readings for any height above ground less than 40 meters. For the purposes of landing based on HAG calculations from optic flow, values above 40 meters were discarded as noise. 21 3.6 Conclusion Landing Using the fixed value resolution sample rate and the 2.65 degree field of view optics, HAG measurements based on optic flow were used to perform twenty seven consecutive autonomous landings from a single approach point. The results are shown in Figure 19. For each of these autonomous landings HAG was calculated according to h = h bar + h offset ( 26) h offset = ( h bar h optic ) + ( 1 ) h offset : ( 27) The value of h was updated every time through the loop and the value of h offset was updated every t s seconds where t s varied with the current HAG estimate and groundspeed. The readings from the barometric pressure sensor were allowed to drift over the twenty seven landings in order to simulate an unknown landing point altitude. The barometric altimeter drifted to 15 meters of inaccuracy over the first 22 landings and was then rezeroed to introduce negative offsets between true HAG and barometric altitude. The value of h offset ( corresponding to the drift in the barometric altimeter) for each of the landings is shown in the third subplot of Figure 19. The second subplot of Figure 23 shows the distance from the landing point for each of the twenty seven landings. Figure 23: Results of 27 Consecutive Optic Flow Based Autonomous Landings. The average distance between the desired and actual touchdown points for autonomous landings using optic flow data was 4.3 meters. This represents a 56% improvement over the average distance from the desired touchdown point using only the barometric altimeter. Landings based on HAG from optic flow also had a much narrower spread than those based on barometric altitude alone with a standard deviation of 2.2 for optic flow based landings compared to 5.4 for barometric. Because the circle down portion of the landing routine is based on barometric altitude until HAG readings fall under a threshold at which they are trusted, it is important to note that the distance from the desired touchdown point does not increase with increasing error in barometric altitude. This demonstrates the ability of the optic flow sensor to land the airplane accurately even when the altitude of the landing point relative to the takeoff point is not known. 3.6 Conclusion The first step necessary to land a vehicle autonomously is the design of control algorithms that will allow the vehicle to accurately perform a landing maneuver. A method for controlling the landing maneuver based on heading fields has been presented. This method has been tested both in simulation and in hardware and has proven effective, repeatable, and robust to wind. The method has been used effectively to perform hundreds of autonomous landings in a wide variety of wind conditions. 22 In addition to effectively tracking a landing maneuver, in order to accomplish true autonomous landing in uncontrolled environments, the UAV must be able to detect HAG reliably. In environments where the terrain is sufficiently flat and the altitude of the landing point relative to the takeoff point is known, a well calibrated barometric altimeter ( with both airflow and temperature compensation) is sufficient to accurately determine HAG and land the UAV. In environments where this is not the case fusion of data from the barometric altimeter and the optic flow sensor provides a feasible and robust method for estimating HAG. This estimate of HAG can be combined with robust landing algorithms to successfully land mini and micro fixed wing UAVs within meters of their desired landing point. No other sensor combination to date with a similar payload and computational cost has demonstrated the same capability. 4 Automatic Trim and Gain Tuning Traditional autopilot design involves taking the equations of flight and linearizing them about an operating point. Controllers are then designed with these linearized models. A popular choice for such controllers are PID controllers. Although such designs are easy to implement and can work well, this method requires that the plane remain within the design flight conditions and are typically not very failure tolerant. To allow for a larger flight envelope, gain scheduling can be used to bridge different PID controllers designed around different operating points. However such designs still require tuning at the different operating points, do not address fault tolerance, and must be tuned for different platforms. Adaptive control was developed to address some of these problems. Examples of different adaptive control approaches are neural network, least squares estimation, and Lyapunov based methods. One adaptive approach is to use neural networks. This often entails training a network offline which is used in the autopilot. However, to compensate for training errors and to adapt to changes in the aircraft ( i. e. battle damage, system failures, etc.), a separate online neural network is also used. This online network is then constantly adapting to changes on the plane and training error. [ 2], [ 3], [ 4], [ 5], and [ 6] discuss this approach. Another approach identifies the plane parameters using recursive least squares techniques. These parameters are then used to adjust the controller. Such controllers have the ability to quickly converge on the plane parameters [ 7], and therefore adapt quickly to failures. [ 8], [ 9], and [ 5] use least squares estimation. Also, a Lyapunov based approach can be used. In this approach parameters are adjusted according to a Lyapunov function to assure stability. These algorithms typically have a small code footprint and require little computing power. [ 10], [ 11], and [ 12] are examples of Lyapunov based algorithms. In our research, we deal with Miniature Aerial Vehicles ( MAVs). We have found that PID values are different for each MAV, even if they have the same geometry. This is mostly due to the manufacturing process. The desire to have one autopilot that works on any of our MAVs with only one set of gains has lead us into adaptive control. Because of the size of our MAVs, we need to use a small autopilot. This means we have code space, memory, and computational restrictions on our autopilot. Also, because of the differing platforms and an imperfect manufacturing process, we do not have an accurate model of our MAVs. This means we are not able to use neural networks because they require a good model of the MAV, require some computational power to implement, and are not well suited for being placed on MAVs of differing geometries. We also cannot use the least squares algorithms which use inverse matrices that require memory and computer power to handle. Such resources come at a premium. This has lead us to use the Lyapunov based approach. Using model reference design and Lyapunov theory, we will show that a stable, adaptive algorithm can be derived that is computationally simple. This new autopilot design will automatically adapt itself to whatever MAV it is flying. Its inputs will be desired pitch and roll angles. We assume that both heading and altitude looped are closed outside of the current control design and feeding the autopilot with the desired pitch and roll angles. It should be noted that adaptive control algorithms are primarily used to compensate for actuator and control surface failure. To the authors knowledge it is not being used to make a universal autopilot that can be used on a variety of different platforms. This makes our approach somewhat unique. Further research should be done to see how many different platforms our autopilot can fly on. 23 4.1 UAV Model In this paper we will show four different Model Reference Adaptive Control ( MRAC) schemes. All four schemes will use the same reference model, but will differ on how much is lumped into the regressor s bias term. The first scheme will employ backstepping to handle derivatives on angular rates. The second does not use backstepping; instead we assume that the angular rate derivatives change slowly enough that we can assume they are slow varying parameters. The third and forth schemes further simplify by lumping more slowly changing signals into the bias term. Flight data from three different MAV designs will also be shown. 4.1 UAV Model If we assume that the inertia matrix has the form J = 0 @ J x 0 J x z 0 J y 0 J x z 0 J z 1 A ; and that the aerodynamic moments are linear, then the rotational dynamics of the UAV are given by _ = p + q s i n t a n + r c o s t a n ( 28) _ = q c o s r s i n ( 29) _ = q s i n s e c + r c o s s e c ( 30) p _ = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V + C p a a + C p r r ( 31) q _ = J x z J y ( r 2 p 2 ) + J z J x J y p r + 1 2 J y V 2 c S C m 0 + C m + C m q c q V + C m e e ( 32) r _ = 3 p q 4 q r + 1 2 V 2 S b 2 C r 0 + C r + C r p b p 2 V + C r r b r 2 V + C r a a + C r r r ; ( 33) where , , roll, pitch, yaw angles d , d desired roll and pitch angles m , m reference model roll and pitch angles k m reference model gain p , q , r roll, pitch, yaw rates , angle of attack, sideslip angle , flight path heading angle, flight path angle e , a , r elevator, aileron, rudder commands V airspeed S , b , c wing area, wing span, average cord length J i j inertial moments about i and j axies Table 1: Definition of terms 24 4.1 UAV Model = ( J x J z J 2 x z ) 1 = J x z ( J x J y + J z ) = ; 2 = ( J z ( J z J y ) + J 2 x z ) = 3 = ( ( J x J y ) J x + J 2 x z ) = 4 = ( J x z ( J x J y + J z ) ) = C p 0 = J z C 0 + J x z C n 0 C p = J z C + J x z C n C p p = J z C p + J x z C n p C p r = J z C r + J x z C n r C p a = J z C a + J x z C n a C p r = J z C r + J x z C n r C r 0 = J x z C 0 + J x C n 0 C r = J x z C + J x C n C r p = J x z C p + J x C n p C r r = J x z C r + J x C n r C r a = J x z C a + J x C n a C r r = J x z C r + J x C n r : For all schemes we use the following first order models, _ m = k m d m ( 34) for pitch and _ m = k m d m ( 35) for roll. We define tacking error as = m and ( 36) = m ( 37) for pitch and roll, respectively. The derivatives are _ = _ _ m = q c o s r s i n _ m and ( 38) _ = _ _ m = p + q s i n t a n + r c o s t a n _ m : ( 39) We also make the following assumptions: 25 4.2 Scheme A: nonlinear MRAC using backstepping A1: The following states can be measured: , , p , r , q . A2: The pitch angle is limited to where < = 2 . A3: The roll angle is limited to where < = 2 . A4: The pitch angle = + where is the flight path climb angle. A5: The pitch angle = where is the flight path heading angle. A6: For the pitch controllers is constant. A7: For the roll controllers is constant. 4.2 Scheme A: nonlinear MRAC using backstepping 4.2.1 Pitch Attitude Hold For pitch attitude hold, we will work with Equations ( 29) and ( 32). Our objective is to design a controller so that the pitch angle follows the reference model ( 34). Our desire is to pick e such that Equation ( 38) becomes _ = 1 ; where 1 > 0 . If we define q d e s 4 = _ m 1 + r s i n c o s ; and let q = q d e s in Equation ( 38), we get the desired result. Using the actual pitch rate but adding and subtracting q d e s we obtain _ = q c o s r s i n _ m = q + q d e s q d e s c o s r s i n _ m = q d e s c o s r s i n _ m + q q d e s c o s = 1 + q c o s ; ( 40) where q 4 = q q d e s . Differentiating q we get _ q = _ q _ q d e s = J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y C m 0 + C m ( ) + C m q c q V + C m e e q _ d e s : Define the Lyapunov function candidate V = 1 2 2 + 1 2 q 2 : 26 4.2 Scheme A: nonlinear MRAC using backstepping Differentiating we get _ V = _ + q _ q = 1 2 + q c o s + q _ q = 1 2 + q q _ d e s + J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y h C m 0 + C m ( ) + C m q c q V + C m e e i : ( 41) If we knew the coefficients, we would like to pick e to satisfy 2 q = q _ d e s + J x z J y ( p 2 r 2 ) + J z J x J y p r + 1 2 V 2 S c J y h C m 0 + C m ( ) + C m q c q V + C m e d e s e i : Solving for d e s e gives d e s e = J y 1 2 V 2 S c C m e 2 q c o s + q _ d e s + 1 1 2 V 2 S c C m e J x z ( p 2 r 2 ) ( J z J x ) p r 1 C m e C m 0 + C m ( ) + C m q c q V : ( 42) If we define K 4 = 0 B B B B B B B B B B B @ J y 1 2 S c C m e J x z 1 2 S c C m e J z J x 1 2 S c C m e C m 0 C m C m e C m C m e C m q c C m e 1 C C C C C C C C C C C A and 4 = 0 B B B B B B @ 2 q c o s + q _ d e s V 2 p 2 r 2 V 2 p r V 2 1 q V 1 C C C C C C A ; ( 43) then d e s e = T K : Let ^ K be the current estimate of K , then we will use the control e = T ^ K : Therefore, plugging into Equation ( 41) we get _ V = 1 2 + q q _ d e s + J x z J y ( p 2 r 2 ) + J z I x J y p r + 1 2 V 2 S c J y h C m 0 + C m + C m q c q V + C m e e + d e s e d e s e i = 1 2 2 q 2 + q 1 2 V 2 S c J y C m e ( e d e s e ) = 1 2 2 q 2 + 1 2 S c C m e J y ! V 2 q T ( ^ K K ) = 1 2 2 q 2 1 2 S c C m e J y ! V 2 q T K ; 27 4.2 Scheme A: nonlinear MRAC using backstepping where K 4 = K ^ K . To cancel the last term, let the modified Lyapunov function candidate be V 2 = 1 2 2 + 1 2 q 2 + 1 2 V S c C m e I y ! 1 2 K T 1 K : ( 44) where is a diagonal matrix. Then _ V 2 = 1 2 2 q 2 1 2 S c C m e I y ! V 2 q T K + 1 2 V S c C m e I y ! _ K T 1 K : Observing that _ K = _ K _ ^ K = _ ^ K since K is constant, gives _ V 2 = 1 2 2 q 2 1 2 V S c C m e I y ! h 1 _ ^ K + V q i T K : ( 45) Therefore letting _ ^ K = V q ; ( 46) gives _ V 2 = 1 2 2 q 2 : The algorithm can be summarized as follows: Algorithm 1 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: q _ d e s 4: Compute: = 2 q c o s + q _ d e s V 2 p 2 r 2 V 2 p r V 2 1 q V T , 5: Update the gain estimate according to: _ ^ K = V q 6: Compute the elevator command: e = T K . 28 4.2 Scheme A: nonlinear MRAC using backstepping 4.2.2 Roll Attitude Hold Using Equation ( 35) as our roll model reference, we want Equation ( 39) to be _ = , so we define p d e s to be p d e s 4 = q s i n t a n r c o s t a n + _ m : We define p as p 4 = p p d e s ( 47) and take its derivative _ p = _ p _ p d e s : ( 48) Substituting equation ( 31) into equation ( 48) we get the following. _ p = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V + C p a a + C p r r p _ d e s ( 49) We define our Lyapunov candidate function V = 1 2 2 + 1 2 p 2 ( 50) and take it derivative _ V = _ + p _ p = p + q s i n t a n + r c o s t a n _ m + p _ p : ( 51) We add and subtract p d e s to get _ V = p + p d e s + q s i n t a n + r c o s t a n _ m + p _ p = 1 + p + p _ p = 1 2 + p + p _ p = 1 2 + p + _ p : ( 52) To arrive at the Lyapunov equation we want, we set + _ p = 2 p . If we now plug this into equation ( 49) and solve for a , calling it d e s a , gives us d e s a = 4 V 2 S b C p a 2 p + p _ d e s 4 1 V 2 S b C p a p q 4 2 V 2 S b C p a q r C p 0 C p a C p C p a C p p C p a b p 2 V C p r C p a b r 2 V C p r C p a r : ( 53) Using the assumption that = where is the flight path heading angle gives us d e s a = 4 V 2 S b C p a 2 p + p _ d e s 4 1 V 2 S b C p a p q 4 2 V 2 S b C p a q r C p 0 C p a C p C p a + C p C p a C p p C p a b p 2 V C p r C p a b r 2 V C p r C p a r : ( 54) which can be represented as d e s a = T K ( 55) 29 4.2 Scheme A: nonlinear MRAC using backstepping where = 0 B B B B B B B B B B @ 2 p + p _ d e s V 2 p q V 2 q r V 2 1 p V r V r 1 C C C C C C C C C C A K = 0 B B B B B B B B B B B B B B B B B @ 4 b V 2 S C a 4 1 b V 2 S C p a 4 2 b V 2 S C p a C p 0 C p a + C p C p a C p C p a b C p p 2 C p a b C p r 2 C p a C p r C p a 1 C C C C C C C C C C C C C C C C C A : ( 56) However, we do not know the terms in K and instead only have an estimate ^ K so really we have a = T ^ K . Adding and subtracting d e s a to equation ( 49) we get _ p = 1 p q 2 q r + 1 2 V 2 S b 2 C p 0 + C p + C p p b p 2 V + C p r b r 2 V ( 57) + C p a a d e s a + d e s a + C p r r p _ d e s _ p = V 2 S b C p a 4 a d e s a 2 p _ p = V 2 S b C p a 4 T ^ K K 2 p _ p = V 2 S b C p a 4 T K 2 p + _ p = V 2 S b C p a 4 T K 2 p ( 58) Defining a new Lyapunov function candidate V 2 = 1 2 2 + 1 2 p 2 + V S b C p a 8 K T 1 K ( 59) where is a diagonal matrix. We differentiate the Lyapunov function to get V _ 2 = _ + p _ p + V S b C p a 4 _ K T 1 K = 1 2 + p + p _ p + V S b C p a 4 _ K T 1 K ( 60) Substituting the results from ( 57) we get V _ 2 = 1 2 + p + _ p + V S b C p a 4 _ K T 1 K = 1 2 + p V 2 S b C p a 4 T K 2 p + V S b C p a 4 _ K T 1 K = 1 2 2 p 2 + p V 2 S b C p a 4 T K + V S b C p a 4 _ K T 1 K = 1 2 2 p 2 + V S b C p a 4 p V T + _ K T 1 K : ( 61) Setting p V T + _ K T 1 = 0 we arrive at the update law _ ^ K = p V : ( 62) The algorithm can be summarized as follows. 30 4.3 Scheme B: nonlinear MRAC without backstepping Algorithm 2 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: p _ d e s 4: Compute: = 2 p + p _ d e s V 2 r _ + p q V 2 q r V 2 1 p V r V r T , 5: Update the gain estimate according to: _ ^ K = p V 6: Compute the aileron command: a = T ^ K . 4.3 Scheme B: nonlinear MRAC without backstepping 4.3.1 Pitch Attitude Hold Using the same definition of tracking error for pitch as above, we design controller without using backstepping. Again, our desire is to pick e such that Equation ( 38) becomes _ = 1 ; ( 63) where 1 > 0 . Solving for q in equation ( 32), and assuming q _ is slow varying parameters, we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 64) If we knew the parameters we could we could command d e s e = T K where T = 0 B B B B @ 1 k m ( d m ) r s i n V c o s 1 p 2 r 2 p r 1 C C C C A K = 0 B B B B B B B B @ c C m p C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e J x z 1 2 V 2 S c C m e J z J x 1 2 V 2 S c C m e C m C m e 1 C C C C C C C C A ( 65) Adding and subtracting this in equation ( 64) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 66) 31 4.3 Scheme B: nonlinear MRAC without backstepping Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 67) We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 68) We summarize the algorithm. Algorithm 3 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: T = 1 k m ( d m ) r s i n V c o s 1 p 2 r 2 p r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the elevator command: e = T ^ K . 32 4.3 Scheme B: nonlinear MRAC without backstepping 4.3.2 Roll Attitude Hold Assuming p _ and r _ are slow varying parameters and solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 69) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 70) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 71) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 72) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us = 0 B B B B B B B B B B @ 1 + q s i n t a n + r c o s t a n k m ( d m ) V 1 V 2 p q V 2 q r V 2 1 r V r 1 C C C C C C C C C C A and K = 0 B B B B B B B B B B B B B B B B @ b C p p 2 C p a 4 S b p _ 4 1 S b 4 2 S b C p 0 C p a C p C p a C p C p a C p r 2 C p a C p r C p a 1 C C C C C C C C C C C C C C C C A : ( 73) 33 4.4 Scheme C: nonlinear MRAC grouping terms Adding and subtracting d e s a from Equation ( 70) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 74) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 75) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 76) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 4 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: = 1 + q s i n t a n + r c o s t a n k m ( d m ) V 1 V 2 p q V 2 q r V 2 1 r V r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 4.4 Scheme C: nonlinear MRAC grouping terms 4.4.1 Pitch Attitude Hold We assume that q _ , 1 c o s , p 2 r 2 , and p r terms are slow varying parameters. Our desire is to pick e such that Equation ( 38) is _ = 1 ; ( 77) 34 4.4 Scheme C: nonlinear MRAC grouping terms where 1 > 0 . Solving for q in equation ( 32) we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 78) If we knew the parameters we could we could command d e s e = T K where T = 0 @ 1 k m ( d m ) r s i n V 1 1 A and K = 0 B B @ c C m p c o s C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e + J x z 1 2 V 2 S c C m e p 2 r 2 J z J x 1 2 V 2 S c C m e p r C m C m e 1 C C A : ( 79) Adding and subtracting this in equation ( 78) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 80) Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 81) We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 82) We summarize the algorithm. 35 4.4 Scheme C: nonlinear MRAC grouping terms Algorithm 5 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d m ) V c o s 1 T , 4: Update the gain estimate according to: _ ^ K = V c o s 5: Compute the elevator command: e = T ^ K . 4.4.2 Roll Attitude Hold Let the derivative of the tracking error be _ = p + q s i n t a n + r c o s t a n _ m : ( 83) We make the assumption here that the terms p _ , r _ , q s i n t a n + r c o s t a n , p q , q r , and r V are slow varying parameters Solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 84) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 85) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 86) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 87) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us = 0 B B @ 1 k m ( d m ) V 1 r 1 C C A and K = 0 B B B B B @ b C p p 2 C p a C p b i a s C p C p a C p r C p a 1 C C C C C A ; ( 88) 36 4.4 Scheme C: nonlinear MRAC grouping terms where C p b i a s = 4 V 2 S b p _ + C p 0 C p a + b C p p 2 C p a ( q s i n t a n + r c o s t a n ) 4 1 S b p q V 2 + 4 2 S b q r V 2 C p C p a C p r 2 C p a r V Adding and subtracting d e s a from Equation ( 89) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 89) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 90) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 91) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 6 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m d m . 3: Compute: = 1 k m ( d ) 2 V 1 r T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 37 4.5 Scheme D: nonlinear MRAC simple 4.5 Scheme D: nonlinear MRAC simple 4.5.1 Pitch Attitude Hold We assume that q _ , 1 c o s , p 2 r 2 , and p r terms are slow varying parameters. We want to pick e such that _ = 1 ; ( 92) where 1 > 0 . Solving for q in equation ( 32) we get q = V C m q c 2 V 2 S c [ I y q _ J x z r 2 p 2 ( I z I x ) p r ] C m 0 + C m C m C e e : If we now plug this into equation ( 38) we get _ = V c o s C m q c 2 V 2 S c [ J y q _ J x z r 2 p 2 ( J z J x ) p r ] C m 0 + C m C m C e e r s i n _ m ( 93) If we knew the parameters we could we could command d e s e = T K where T = 1 k m ( d m ) V 1 ! T and K = 0 B B @ c C m p c o s C m e J y q _ 1 2 S c C m e C m o C m e + C m C m e ( ) + J x z 1 2 V 2 S c C m e p 2 r 2 J z J x 1 2 V 2 S c C m e p r C m C m e 1 C C A : ( 94) Adding and subtracting this in equation ( 95) we get _ = V c o s C m q c 2 V 2 S c J y q _ J x z r 2 p 2 ( J z J x ) p r C m 0 + C m C m C e d e s e + e d e s e r s i n _ m = 1 + V c o s C m q C m e c T K ( 95) Assuming is constant, we define our Lyapunov function candidate as V 4 = 1 2 2 + c o s 2 C m q C m e c K T 1 K and take its derivative we get _ V = _ + c o s C m q C m e c _ K T 1 K = 1 2 + V c o s C m q C m e c T K + c o s C m q C m e c _ ^ K T 1 K = 1 2 + c o s C m q C m e c V T + _ ^ K T 1 K ( 96) 38 4.5 Scheme D: nonlinear MRAC simple We now let V T + _ ^ K T 1 = 0 and we get our update law _ ^ K = V ( 97) We summarize the algorithm. Algorithm 7 Nonlinear MRAC for pitch attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d m ) V 1 T , 4: Update the gain estimate according to: _ ^ K = V c o s 5: Compute the elevator command: e = T ^ K . 4.5.2 Roll Attitude Hold We make the assumption here that the terms p _ , r _ , q s i n t a n + r c o s t a n , p q , q r , and r V are slow varying parameters Solving for p in equation ( 31) we find that p = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r : ( 98) Substituting this equation in equation ( 39) give us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m : ( 99) But since we want _ = 1 to guarantee that the tracking error goes to 0, we define d e s a to be d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 100) Using the assumption that = , then d e s a = b C p p 2 V C p a 1 _ + q s i n t a n + r c o s t a n _ m C p 0 C p a C p C p a + C p C p a C p r C p a r C p r 2 V C p a r + 4 V 2 S b p _ 4 1 V 2 S b p q + 4 2 V 2 S b q r ( 101) Grouping the knowns and unknowns together so we get the form d e s a = T K gives us T = 0 @ 1 k m ( d ) V 1 r 1 A T and K = 0 B B @ b C p p 2 C p a C p b i a s C p r C p a 1 C C A ; ( 102) 39 4.5 Scheme D: nonlinear MRAC simple where C p b i a s = 4 V 2 S b p _ + C p 0 C p a + b C p p 2 C p a ( q s i n t a n + r c o s t a n ) 4 1 S b p q V 2 + 4 2 S b q r V 2 C p C p a ( ) C p r 2 C p a r V Adding and subtracting d e s a from Equation ( 99) gives us _ = 8 V S b 2 C p p p _ 8 1 V S b 2 C p p p q + 8 2 V S b 2 C p p q r 2 V C p 0 b C p p 2 V C p b C p p C p r C p p r 2 V C p a b C p p a + d e s a d e s a 2 V C p r b C p p r + q s i n t a n + r c o s t a n _ m = 1 2 V C p a b C p p a d e s a = 1 2 V C p a b C p p T ^ K T K = 1 2 V C p a b C p p T K : ( 103) If we define our Lyapunov candidate function to be V = 1 2 2 + V C p a b C p p K T 1 K ( 104) and take its derivative to get _ V = _ + 2 V C p a b C p p K T 1 K = 1 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 V C p a b C p p T K + 2 V C p a b C p p K T 1 K = 1 2 2 C p a b C p p V T + K T 1 K : ( 105) Setting V T + K T 1 = 0 we arrive at our update law _ ^ K = V : The algorithm can be summarized as follows. Algorithm 8 Nonlinear MRAC for roll attitude hold. 1: Obtain , p , r , , q , from sensors. 2: Update the reference model according to _ m = k m ( d m ) . 3: Compute: = 1 k m ( d ) V 1 T , 4: Update the gain estimate according to: _ ^ K = V 5: Compute the aileron command: a = T ^ K . 40 4.6 Flight Data 4.6 Flight Data 4.6.1 Platform The algorithms were flown on MAVs equipped with the Kestrel Autopilot running a 29 MHz Rabbit microcontroller with 512K Flash and 512K RAM. The autopilot is equipped with a variety of sensors, including rate gyros, accelerometers, an absolute pressure sensor for measuring altitude, a differential pressure sensor for measuring airspeed, and a GPS receiver. The autopilot is programmed in C. The adaptive control algorithms run at about 80 90 hertz with this platform. The airframes flown were a Zagi style flying wing design ( Deliverance), a larger and heavier flying wing design ( Ray), and a fixed wing glider design ( Phidipides). In the initial phase of testing, we used Deliverance to tune adaptive control gains and to obtain a comparison between the four algorithms presented in this paper. Afterward, we narrowed our testing to one algorithm ( D) and flew with the other airframes. Deliverance Weight: 0.734 k g Wingspan: 1.207 m Mean geometric chord: 0.254 m Ray Weight: 1.321 k g Wingspan: 1.511 m Mean geometric chord: 0.305 m Phidipides Weight: 0.969 k g Wingspan: 1.524 m Mean geometric chord: 0.168 m Nose to tail: 0.914 m 41 4.6 Flight Data Pitch k m 3.0 Pitch 1 45.0 Pitch 1 0.06 Pitch 1 0.01 Pitch Leakage gain 1 0.001 Pitch Leakage gain 2 0.0005 Roll k m 4.0 Roll 1 140.0 Roll 1 0.005 Roll 1 0.001 Roll Leakage gain 1 0.001 Roll Leakage gain 2 0.0001 Table 2: Adaptive control gains used in Algorithm D 4.6.2 Implementation While Lyapunov theory provides for stability, adaptive control engineers have several techniques available to ensure good stability in practice. Many times, flight conditions do not provide sufficient excitation for all of the adaptive parameters. Sensor noise and other disturbances can also cause parameter drift. [ 13] Some of the parameters may diverge such that when flight conditions change quickly ( i. e., we switch from straight, level flight into a climbing turn), the plane can go unstable. We use three techniques in our implementation to improve stability: dead zone, leakage, and normalization. With dead zone, adaption is stopped once the actual value gets within a certain region around the desired value. Once we get into a flight pattern with little excitation, the ^ K parameters will adapt until the plane is close to the desired value. Then adaption freezes until flight conditions change such that we are no longer within the no adaption ( dead ) zone. Leakage subtracts off a small percentage of each ^ K parameter each time the code is executed. This prevents parameters from ever becoming too large. So as not to affect performance, the leakage gain is kept small, about 0.0001 for most parameters. Normalization is similar to leakage. Each time through the control loop, all of the ^ K parameters are normalized with respect to themselves. This prevents one parameter from become much larger than the others, but preserves the relative weighting between each parameter. Because all of our MAVs do not have a rudder we have implemented the controller with r = 1 . 4.6.3 Results We wanted to compare the four algorithms against each other and against the usual PD controller used in our lab to control pitch and roll on the MAVs. To obtain this comparison, we flew an hourglass shaped pattern with legs between 200 and 300 meters long on the same airframe ( Deliverance) with each of the four algorithms. The flights incorporated auto takeoff and auto land functions, giving the adaptive algorithms full control for the entire flight. We also obtained data for the flight plan flown with more conventional PD control, with gains tuned by hand. For all of the adaptive algorithms, transients die out within about 3 seconds for pitch and about 10 seconds for roll. Along with the adaptive control plots, we show the values of the ^ K parameters adapting with time. In most cases, these parameters settle to a steady value after adaption has been running for a while. Note that while the parameters represent certain combinations of aerodynamic factors and other estimated quantities, they do not, in most cases, converge to these true values. The Lyapunov stability proof used to derive the algorithms does not guarantee that these parameters will converge to the true values. Rather, the proof shows that the parameters stabilize to values that will cause the tracking error to go to zero. The results show that each adaptive control algorithm s performance is comparable to PD control. Before each algorithm could be flown, several adaptive control gains ( namely k m , 1 , and ) had to be tuned. Algorithm D, with the fewest number of gains, was easiest to tune. Algorithm D also is the simplest algorithm to implement and takes up the smallest amount of autopilot code space. 42 4.7 Conclusion Algorithm Pitch Roll A 2.0 2.5 B 2.0 3.5 C 2.0 5.2 D 2.0 3.9 PD control 1.5 3.5 Table 3: Average error in degrees for each algorithm  from last third of hourglass flight tests on Deliverance We used average error as a metric to compare between the algorithms. We define average error as the average of p ( d e s i r e d a c t u a l ) 2 over the sample. Each figure included in this paper displays average error over the total sample as well as the average error over each third of the flight time. This shows the gradual improvement in performance characteristic of adaptive control as the parameters become more finely tuned. In most of the figures, the last third of the flight has less average error than the first third. The longer the adaptive algorithms run, the better the performance will be. The adaptive gains k m , 1 , 2 , and were very difficult to tune for most of the algorithms, especially algorithm C. Algorithms A and B were tuned to give good performance, but did not have very good stability. Based on low average error, good stability, and simplicity of implementation, we chose to focus on algorithm D for further testing. None of the algorithms outperformed PD control, but bear in mind that the control gains were finely tuned to Deliverance s airframe. To test the robustness of adaptive control, we attached a flap near the center of Deliverance s right wing and deployed it mid flight. By doing this, we simulated what a change in the airframe may have on the controllability of the MAV. With the flap deployed, PD control can still fly the MAV, but with substantial steady state error on roll. However, algorithm D can be seen to adapt to the change in aerodynamics within a few seconds. When the flap goes back down, the parameters once again adapt back to their previous states. Flap attached to Deliverance To further test adaptive control s robustness, we flew two very different airframes; Ray, a heavier Zagi style airframe, and Phidipides, a fixed wing glider design. The plots presented in this paper are flight data from flying Algorithm D with the adaptive gains tuned while flying Deliverance. Both Ray and Phidipides fly well with these same gains. Phidipides in particular was very hard to tune PD gains for, but algorithm D had no problems. This shows the main strength of adaptive control: the ability to fly the same code on different airframes without having to hand tune gains. 4.7 Conclusion In this paper we have shown the derivation for four adaptive control schemes for pitch and roll. We also present flight test data for each of the schemes. Through comparing performance, stability, and simplicity, Algorithm D was chosen as the best algorithm to pursue for further testing and development on the autopilot. We show that the algorithm performs well on several different airframes without having to tune any gains. 43 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm A ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 24: Algorithm A pitch on Deliverance flying an hourglass 44 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm A ( Pitch) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 25: Algorithm A pitch parameters on Deliverance flying an hourglass 45 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm A ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 26: Algorithm A roll on Deliverance flying an hourglass 46 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm A ( Roll) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Figure 27: Algorithm A roll parameters on Deliverance flying an hourglass 47 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm B ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 28: Algorithm B pitch on Deliverance flying an hourglass 48 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm B ( Pitch) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 29: Algorithm B pitch parameters on Deliverance flying an hourglass 49 4.7 Conclusion 50 100 150 200 250 300 Time ( seconds) Algorithm B ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 30: Algorithm B roll on Deliverance flying an hourglass 50 4.7 Conclusion 50 100 150 200 250 300 Estimated Parameters of Algorithm B ( Roll) Controlling Deliverance 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 50 100 150 200 250 300 Time ( seconds) Figure 31: Algorithm B roll parameters on Deliverance flying an hourglass 51 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm C ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 32: Algorithm C pitch on Deliverance flying an hourglass 52 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm C ( Pitch) Controlling Deliverance 50 100 150 200 250 300 350 50 100 150 200 250 300 350 Time ( seconds) Figure 33: Algorithm C pitch parameters on Deliverance flying an hourglass 53 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm C ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 34: Algorithm C roll on Deliverance flying an hourglass 54 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm C ( Roll) Controlling Deliverance 50 100 150 200 250 300 350 50 100 150 200 250 300 350 Time ( seconds) Figure 35: Algorithm C roll parameters on Deliverance flying an hourglass 55 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm D ( Pitch) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 36: Algorithm D pitch on Deliverance flying an hourglass 56 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm D ( Pitch) Controlling Deliverance 50 100 150 200 250 300 350 Time ( seconds) Figure 37: Algorithm D pitch parameters on Deliverance flying an hourglass 57 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) Algorithm D ( Roll) Controlling Deliverance Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 38: Algorithm D roll on Deliverance flying an hourglass 58 4.7 Conclusion 50 100 150 200 250 300 350 Estimated Parameters of Algorithm D ( Roll) Controlling Deliverance 50 100 150 200 250 300 350 Time ( seconds) Figure 39: Algorithm D roll parameters on Deliverance flying an hourglass 59 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) PD ( Pitch) Controlling Deliverance Desired Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 40: PD pitch on Deliverance flying an hourglass 60 4.7 Conclusion 50 100 150 200 250 300 350 Time ( seconds) PD ( Roll) Controlling Deliverance Desired Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 41: PD roll on Deliverance flying an hourglass 61 4.7 Conclusion 50 100 150 200 250 Algorithm D ( Pitch) Controlling Deliverance as Flap Deploys Desired Model Actual 50 100 150 200 250 Algorithm D ( Roll) Controlling Deliverance as Flap Deploys Desired Model Actual 50 100 150 200 250 Flap Time ( seconds) Figure 42: Algorithm D Pitch and Roll with Flap Deployment flying an hourglass. One transit around the hourglass is completed with the flap deployed. 62 4.7 Conclusion 50 100 150 200 250 Estimated Parameters of Algorithm D Controlling Deliverance as Flap Deploys 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 Flap Time ( seconds) Figure 43: Algorithm D Pitch and Roll parameters with Flap Deployment flying an hourglass. Notice k 2 for roll adapting to a new value when the flap is deployed and adapting back to the previous value when the flap is stowed. 63 4.7 Conclusion 50 100 150 200 250 PD ( Pitch) Controlling Deliverance as Flap Deploys Desired Actual 50 100 150 200 250 PD ( Roll) Controlling Deliverance as Flap Deploys Desired Actual 50 100 150 200 250 Flap Time ( seconds) Figure 44: PD Control with Flap Deployment flying an hourglass. Notice the steady state error of about 10 degrees in roll. 64 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Time ( seconds) Algorithm D ( Pitch) Controlling Ray 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Desired Model Actual Figure 45: Algorithm D pitch on Ray flying between two waypoints 65 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Estimated Parameters of Algorithm D ( Pitch) Controlling Ray 20 40 60 80 100 120 140 160 180 Time ( seconds) Figure 46: Algorithm D pitch parameters on Ray flying between two waypoints 66 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Time ( seconds) Algorithm D ( Roll) Controlling Ray Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 47: Algorithm D roll on Ray flying between two waypoints 67 4.7 Conclusion 20 40 60 80 100 120 140 160 180 Estimated Parameters of Algorithm D ( Roll) Controlling Ray 20 40 60 80 100 120 140 160 180 Time ( seconds) Figure 48: Algorithm D roll parameters on Ray flying between two waypoints 68 4.7 Conclusion 50 100 150 200 250 300 350 400 Time ( seconds) Algorithm D ( Pitch) Controlling Phidipidies Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 49: Algorithm D pitch on Phidipides flying an hourglass 69 4.7 Conclusion 50 100 150 200 250 300 350 400 Estimated Parameters of Algorithm D ( Pitch) Controlling Phidipides 50 100 150 200 250 300 350 400 Time ( seconds) Figure 50: Algorithm D pitch parameters on Phidipides flying an hourglass 70 4.7 Conclusion 50 100 150 200 250 300 350 400 Time ( seconds) Algorithm D ( Roll) Controlling Phidipidies Desired Model Actual 1 2 3 4 1: Total Flight 2: First Third of Flight 3: Second Third of Flight 4: Last Third of Flight Average Error Figure 51: Algorithm D roll on Phidipides flying an hourglass 71 4.7 Conclusion 50 100 150 200 250 300 350 400 Estimated Parameters of Algorithm D ( Roll) Controlling Phidipides 50 100 150 200 250 300 350 400 Time ( seconds) Figure 52: Algorithm D roll parameters on Phidipides flying an hourglass 72 5 Students Funded Under Project Student Topic Steve Griffiths Magnetometer integration Blake Barber Height above ground sensor Andrew Eldridge Magnetometer integration Josh Matthews Adaptive control Nathan Knoebel Adaptive control Steve Osborne Adaptive control References [ 1] J. Chahl, M. V. Srinivasan, and S. W. Zhang, Landing strategies in honeybees and applications to uninhabited airborne vehicles, The International Journal of Robotics Research, vol. 23, no. 2, pp. 101 110, 2004. [ 2] D. Shin and Y. Kim, Reconfigurable flight control system design using adaptive neural networks, 12, no. 1, Jan. 2004. [ 3] R. L. Broderick, Statistical and adaptive approach for verification of a neural based flight control system, in Digital Avionics Systems Conference, vol. 2, pp. 6. E. 1 1 6. E. 1 10, 24 28 Oct. 2004. [ 4] P. Melin and O. Castillo, A new neuro fuzzy fractal approach for adaptive model based control of non linear dynamic systems: The case of controlling aircraft dynamics, in 1999 IEEE International FUzzy Systems Conference Proceedings, August 1999. [ 5] M. Steinberg and A. Page, High fidelity simulation testing of intelligent and adaptive aircraft control laws, in Proceedings of the American Control Conference, pp. 3264 3268, 8 10 May 2002. [ 6] B. Kim and A. Calise, Nonlinear flight control using neural networks, IEEE Transactions on Control Systems Technology, vol. 20, no. 1, Jan. Feb. 1997. [ 7] M. Bodson and J. Groszkiewicz, Multivariable adaptive algorithm for reconfigurable flight control, IEEE Transactions on Control Systems Technology, vol. 5, no. 2, Mar. 1997. [ 8] D. Shore and M. Bodson, Flight testing of a reconfigurable control system on an unmanned aircraft, in Proceedings of the 2004 American Control Conference, pp. 3747 3752, June July 2004. [ 9] B. Porter and C. Boddy, Design of adaptive digital controllers incorporating dynamic pole assignment compensators for high performance aircraft, in Proceedings of the IEEE 1989 Aerospace and Electronics Conference, pp. 372 379, 22 26 May 1989. [ 10] G. Tao, S. Chen, J. Fei, and S. Joshi, An adaptive actuator failure compensation scheme for controlling a morphing aircraft model, in Proceedings of the 42nd IEEE Conference on Decision and Control, pp. 4926 4931, Dec. 2003. [ 11] S. Chen, G. Tao, and S. Joshi, An adaptive actuator failure compensation controller for mimo systems, in Proceedings of the 41st IEEE Conference on Decision and Control, pp. 4778 4783, December 2002. [ 12] J. Farrell, M. Polycarpou, and M. Sharma, Adaptive backstepping with magnitude, rate, and bandwidth constraints: Aircraft lonitude control, in Proceedings of the American Control Conference, pp. 3898 3904, June 2003. [ 13] K. Astr om and B. Wittenmark, Adaptive Control. Addison Wesley Publishing Company, 1989. 73 A Appendix: PNI Micromag Specifications 74 75 76 



1 

A 

B 

C 

D 

E 

F 

G 

H 

I 

J 

K 

L 

M 

N 

O 

P 

R 

S 

T 

U 

V 

W 

Y 

Z 


