Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Robust adaptive control for unmanned aerial vehicles
(USC Thesis Other)
Robust adaptive control for unmanned aerial vehicles
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ROBUSTADAPTIVECONTROLFORUNMANNED AERIALVEHICLES by NazlıE.Kahveci ADissertationPresentedtothe FACULTYOFTHEGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (ELECTRICALENGINEERING) December2007 Copyright 2007 NazlıE.Kahveci DEDICATION Tomybelovedparents Nes . e&AliRızaKahveci N.E.K. ii ACKNOWLEDGMENTS I am grateful to my advisor, Dr. Petros A. Ioannou for providing me with the chal- lenging opportunity to be his research assistant during the course of my graduate stud- ies. It has been an exciting and enriching experience for it has allowed me to benefit greatly from his expertise and guidance. I am also indebted to my dissertation commit- tee members, Dr. Michael G. Safonov, Dr. Edmond A. Jonckheere and Dr. Larry G. Redekopp. Theirfeedbackanddevotionoftimearehighlyappreciated. It has been a pleasure to be a member of the Center for Advanced Transportation Technologies at the University of Southern California. I would like to thank Dr. Maj D. Mirmirani for his comments on my research during our collaborations in aerospace controlapplications,Dr. LeonardM.Silvermanforourmathematicaldiscussionsinthe earlierstagesofmystudies,Dr. KonstantinosPsounisforhissuggestionswhileserving as a member of my guidance committee, and Dr. James E. Moore, II for memorable meetingsontransportationengineering. TheControlsandDynamicsBranchinNationalAeronauticsandSpaceAdministra- tionDrydenFlightResearchCentersuppliedtheflighttestdataduringourinteractions. I wouldliketothankthemfortheirhelpinkeepingmystudiesconsistentwiththeirexper- imentalresults. IalsowouldliketoacknowledgeDr. MatthewC.Turnerfordiscussions ontheiranti-windupcompensatordesignincollaborationwithDr. IanPostlethwaite. iii I consider myself fortunate to have attended numerous control systems lectures in Electrical Engineering Department. They have been nurturing and stimulating due to theprofoundknowledgeofthefacultyandhavesignificantlycontributedtomyresearch efforts. I also wish to thank all my colleagues for being inspiring friends with many of whomIhaveenjoyedinterestingtechnicalconversationsovertheyears. Last but not least I would like to thank my family. My parents, Nese and Ali Riza Kahveci have always been an immense influence on my perspectives and an invaluable source of encouragement throughout my education along with my sister, Harika. I owe my sincerest gratitude to my grandfather, Osman Kahveci who is sorely missed but whosememorywillcontinuetokeepmemotivatedtoexpandmyhorizonsfurther. iv TABLEOFCONTENTS Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii ListofTables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii ListofFigures . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 OptimumTrajectoryGenerationforStaticSoaringofUAVs . . . . . . . . 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 AircraftModelandControlObjective . . . . . . . . . . . . . . . . . . 11 2.3 OptimumTrajectoryGeneration . . . . . . . . . . . . . . . . . . . . . 17 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3 AdaptiveLinearQuadraticControlDesignwithAnti-windup Augmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 AdaptiveControlDesignforAircraftSubjecttoActuatorSaturation . . 29 3.2.1 LQControlDesignfortheNominalAircraftModel . . . . . . . 29 3.2.2 RobustAdaptiveLaw . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.3 AdaptiveAnti-windupCompensatorDesign . . . . . . . . . . . 39 3.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.1 VerificationoftheDecide-To-ThermalAlgorithm . . . . . . . . 47 3.3.2 Adaptive LQ Controller for the Unknown Aircraft Model with NoInputConstraints . . . . . . . . . . . . . . . . . . . . . . . 49 3.3.3 EffectsofSaturationLimitsontheOptimalAutonomous SoaringFlightProblem . . . . . . . . . . . . . . . . . . . . . . 50 3.3.4 OptimalSoaringFlightPerformancewithAdaptiveAnti-windup CompensatorDesign . . . . . . . . . . . . . . . . . . . . . . . 50 v 3.3.5 Flight Performance Using an Alternative Method: Including a FiltertoAvoidSaturation . . . . . . . . . . . . . . . . . . . . . 51 3.3.6 ActuatorFailureAccommodation . . . . . . . . . . . . . . . . 52 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4 A Stochastic Approach to Optimal Soaring Problem and Robust Adaptive LQGControl. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 PerformanceLossDuetoIncorrectThermalData . . . . . . . . . . . . 63 4.3 Performance Loss Due to Vertical Gust Disturbances in Inter-thermal Glide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 PerformanceLossDuetoSensorInaccuracies . . . . . . . . . . . . . . 71 4.5 PerformanceLossRecovery: AdaptiveLQGDesign . . . . . . . . . . . 74 4.6 AdaptiveLQRandAdaptiveLQGPerformanceComparison . . . . . . 78 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5 OptimalStaticSoaringofUAVsUsingVehicleRoutingwithTime Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 AnAnalogy: StaticSoaringProblemasaVRP . . . . . . . . . . . . . 84 5.3 StaticSoaringona2DMap . . . . . . . . . . . . . . . . . . . . . . . . 85 5.4 MathematicalFormulationofOptimalStaticSoaringasaVRPTW . . . 88 5.5 FlightArcCostandOptimalFlightVelocity . . . . . . . . . . . . . . . 90 5.6 AnAlgorithmicFrameworkBasedonShortestPathSearch . . . . . . . 92 5.7 ASimulationScenario . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.8 OtherVRPInterpretationsforStaticSoaring . . . . . . . . . . . . . . . 103 5.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6 AHeuristicSearchAlgorithmforManeuveringofUAVsacrossDense ThermalAreas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Soaring Assigned Set of Thermals and Mathematical Formulation as a CVRP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 AParallelSavingsBasedHeuristicasaSolutionMethodology . . . . . 110 6.3.1 Sectorization . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.3.2 Kruskal’sAlgorithm . . . . . . . . . . . . . . . . . . . . . . . 111 6.3.3 RouteConstruction . . . . . . . . . . . . . . . . . . . . . . . . 112 6.3.4 RelaxationandParallelSavingsBasedHeuristics . . . . . . . . 113 6.3.5 Near-OptimalRouteSelection . . . . . . . . . . . . . . . . . . 115 6.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 vi 7 GeneticAlgorithmsforShortestPathRoutingofAutonomousGliders . . 122 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 7.2 OptimizationinTermsofBiologicalEvolution . . . . . . . . . . . . . 124 7.3 MathematicalFormulationofThermalSoaringProblem: TheVRP Revisited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.4 AGA-BasedSPAlgorithmAdaptedtoThermalSoaringProblem . . . 127 7.4.1 PopulationInitialization . . . . . . . . . . . . . . . . . . . . . 128 7.4.2 Selection(Reproduction) . . . . . . . . . . . . . . . . . . . . . 128 7.4.3 Crossover(Recombination) . . . . . . . . . . . . . . . . . . . 129 7.4.4 Mutation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.4.5 RepairFunction . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.5 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 8 StabilityAnalysisforIndirectAdaptiveControlDesignwith Anti-windupCompensation . . . . . . . . . . . . . . . . . . . . . 139 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 8.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 8.3 NominalControlDesign . . . . . . . . . . . . . . . . . . . . . . . . . 142 8.4 AnAdaptiveLQControlDesignwithAnti-windupAugmentation . . . 148 8.5 StabilityAnalysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 8.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 9 FurtherResults: IndirectAdaptiveControlforSystemswithInput RateSaturation . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 9.2 InputRateSaturation . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 9.3 AnIndirectAdaptiveControlDesign . . . . . . . . . . . . . . . . . . . 166 9.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 10 SummaryandConcludingRemarks . . . . . . . . . . . . . . . . . . 170 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 vii LISTOFTABLES 5.1 FlightTimeTable(I):SimulationScenarioinSection5.7 . . . . . . . . 100 5.2 FlightTimeTable(II):SimulationScenarioinSection5.7 . . . . . . . . 100 5.3 FlightTimeTable(III):SimulationScenarioinSection5.7 . . . . . . . 101 6.1 FlightTimeTable: SimulationScenarioinSection6.4 . . . . . . . . . . 118 7.1 FlightTimeTable(I):SimulationScenarioinSection7.5 . . . . . . . . 133 7.2 FlightTimeTable(II):SimulationScenarioinSection7.5 . . . . . . . . 134 viii LISTOFFIGURES 2.1 Polarcurveinterpolation(Section2.2) . . . . . . . . . . . . . . . . . . 16 2.2 Thermalprofile:W=15ft/sec, R=1ft(Section2.3) . . . . . . . . 18 2.3 Decide-To-ThermalAlgorithm(Section2.3) . . . . . . . . . . . . . . . 20 3.1 Adaptive LQ Control Structure with Anti-windup Augmentation (Sec- tion3.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2 Distancecoveredviathermalsoaringandpureglidefortheknownmodel caseintheabsenceofsaturationlimits(Section3.3.1) . . . . . . . . . . 54 3.3 Altitudechangeviathermalsoaringandpureglidefortheknownmodel caseintheabsenceofsaturationlimits(Section3.3.1) . . . . . . . . . . 54 3.4 Optimal velocity defining the reference, horizontal velocity of the air- craft,andthermalstrengthasfunctionsoftime(Section3.3.2) . . . . . 55 3.5 Controlinputdirectlyappliedtotheactuator(Section3.3.2) . . . . . . 55 3.6 Instabilityofthesystemresponseforthesaturatedactuatorcase(Section 3.3.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.7 Control input sent to the actuator for the unknown model where satura- tionlimitsaretobeapplied(Section3.3.3) . . . . . . . . . . . . . . . . 56 3.8 Horizontal velocity, V x of the unknown model using the adaptive anti- windupaugmentationinthepresenceofsaturation(Section3.3.4) . . . 57 3.9 Controlinputmodifiedfortheunknownmodeltorecoveroptimalsoar- ingperformance(Section3.3.4) . . . . . . . . . . . . . . . . . . . . . 57 3.10 Theanti-windupcompensatormodification: y aw 1 (Section3.3.4) . . . 58 3.11 Theanti-windupcompensatormodification: y aw 2 (Section3.3.4) . . . 58 ix 3.12 Apropersecondorderprefilterchosentoavoidsaturationbysmoothing the transitions in the reference signal for the control of the unknown model(Section3.3.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.13 Controlinputusingaproperprefiltersoastoavoidthesaturationlimits (Section3.3.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.14 Controlefficiencylossby20%at t=15sec withnosignificantperfor- mancedegradationobserved(Section3.3.6) . . . . . . . . . . . . . . . 60 3.15 For the unknown model in the presence of saturation and control effi- ciencyloss,controlinputwithanti-windupmodification(Section3.3.6) 60 4.1 The flight distance covered for different thermal data provided (Section 4.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Thecorrespondingaltitudelossfordifferentthermaldataprovided(Sec- tion4.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3 The flight distance covered using the adaptive LQR in the presence of verticalgust(Section4.3) . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.4 ThecorrespondingaltitudelossusingtheadaptiveLQRinthepresence ofverticalgustdisturbances(Section4.3) . . . . . . . . . . . . . . . . 69 4.5 The flight velocity using the adaptive LQR in the presence of vertical gustdisturbances(Section4.3) . . . . . . . . . . . . . . . . . . . . . . 70 4.6 ThecontrolinputusingtheadaptiveLQRinthepresenceofverticalgust disturbances(Section4.3) . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.7 The flight distance covered using the adaptive LQR in the presence of verticalgustdisturbancesandmeasurementnoise(Section4.4) . . . . . 72 4.8 ThecorrespondingaltitudelossusingtheadaptiveLQRinthepresence ofverticalgustdisturbancesandmeasurementnoise(Section4.4) . . . 72 4.9 The flight velocity using the adaptive LQR in the presence of vertical gustdisturbancesandmeasurementnoise(Section4.4) . . . . . . . . . 73 4.10 ThecontrolinputusingtheadaptiveLQRinthepresenceofverticalgust disturbancesandmeasurementnoise(Section4.4) . . . . . . . . . . . . 73 4.11 AdaptiveLQGControlStructure(Section4.5) . . . . . . . . . . . . . . 77 x 4.12 The flight distance comparison in the presence of vertical gust distur- bancesandmeasurementnoise(Section4.6) . . . . . . . . . . . . . . . 79 4.13 The corresponding altitude loss comparison in the presence of vertical gustdisturbancesandmeasurementnoise(Section4.6) . . . . . . . . . 79 4.14 The flight velocity comparison in the presence of vertical gust distur- bancesandmeasurementnoise(Section4.6) . . . . . . . . . . . . . . . 80 4.15 The control input comparison in the presence of vertical gust distur- bancesandmeasurementnoise(Section4.6) . . . . . . . . . . . . . . . 80 5.1 The UAV assigned to autonomously navigate from start to destination and the set of thermals detected to be available for soaring in the flight region(Section5.7) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 StaticSoaringProblemona2Dgraph(Section5.7) . . . . . . . . . . . 102 5.3 The thermals assigned by the solution method in order to achieve opti- malsoaringflightperformance(Section5.7) . . . . . . . . . . . . . . . 103 6.1 The thermal map with the direct flight distances between the thermal locationsdepictedsothatthetimecostsofthecorrespondingflightarcs canbelisted(Section6.4) . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.2 Thesolutionalgorithmiteratedonthethermalmap(Section6.4) . . . . 119 6.3 Anear-optimalflightpathreachedusingtheproposedsolutionmethod- ology(Section6.4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.1 Thefeasibledirectflightdistancesbetweenthethermallocationsshown onthethermalmapsothatthetimecostsoftherespectiveflightarcscan becalculated(Section7.5) . . . . . . . . . . . . . . . . . . . . . . . . 131 7.2 The time costs of the feasible flight arcs displayed on the thermal map (Section7.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7.3 TheoptimalsolutionobtainedbyiteratingtheGA(Section7.5) . . . . . 136 8.1 Thenominalcontrolstructurewithanti-windupcompensator(Section8.3)147 8.2 The adaptive control structure with anti-windup compensator (Section 8.4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 9.1 Arate-limiterwithfirst-orderlag(Section9.2) . . . . . . . . . . . . . . 161 9.2 Aninfinitelyfastactuatormodel(Section9.2) . . . . . . . . . . . . . . 161 xi 9.3 Therate-limitedsystem(Section9.2) . . . . . . . . . . . . . . . . . . . 162 9.4 Thecontrolstructurewithanti-windupcompensatorforrate-limitedsys- tems(Section9.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 9.5 The adaptive control structure for rate-limited systems with unknown plantparameters(Section9.3) . . . . . . . . . . . . . . . . . . . . . . 168 xii ABSTRACT Theobjectiveofmeetinghigherendurancerequirementsremainsachallengingtask foranytypeandsizeofUnmannedAerialVehicles(UAVs). Accordingtorecentresearch studiessignificantenergysavingscanberealizedthroughutilizationofthermalcurrents. Thenavigationstrategiesfollowedacrossthermalregions,however,arebasedonrather intuitiveassessmentsofremotepilotsandlackanysystematicpathplanningapproaches. Various methods to enhance the autonomy of UAVs in soaring applications are investi- gatedwhileseekingguaranteesforflightperformanceimprovements. The dynamics of the aircraft, small UAVs in particular, are affected by the environ- mental conditions, whereas unmodeled dynamics possibly become significant during aggressiveflightmaneuvers. Besides,thedemandedcontrolinputsmighthaveamagni- tuderangebeyondthelimitsdictatedbythecontrolsurfaceactuators. Theconsequences ofignoringtheseissuescanbecatastrophic. SupportingthisclaimNASADrydenFlight Research Center reports considerable performance degradation and even loss of stabil- ity in autonomous soaring flight tests with the subsequent risk of an aircraft crash. The existingcontrolschemesareconcludedtosufferfromlimitedperformance. Consideringtheaircraftdynamicsandthethermalcharacteristicswedefineavehicle- specific trajectory optimization problem to achieve increased cross-country speed and extended range of flight. In an environment with geographically dispersed set of ther- mals of possibly limited lifespan, we identify the similarities to the Vehicle Routing Problem (VRP) and provide both exact and approximate guidance algorithms for the xiii navigation of automated UAVs. An additional stochastic approach is used to quantify the performance losses due to incorrect thermal data while dealing with random gust disturbancesandonboardsensormeasurementinaccuracies. Oneofthemaincontributionsofthisresearchisanoveladaptivecontroldesignwith anti-windupcompensation. Ouranalysisontheindirectadaptiveschemerevealsthatthe perturbation terms due to parameter errors do not cause any unbounded signals in the closed-loop. Thestabilityoftheadaptivesystemisestablished,andthepropertiesofthe proposed control scheme are demonstrated through simulations on a UAV model with input magnitude saturation constraints. The robust adaptive control design is further developedtoextendourresultstorate-saturatedsystems. xiv 1. INTRODUCTION AvarietyofUAVsarecurrentlyemployedforabroadrangeofmissionsincludingrecon- naissance, surveillance, search and rescue, law enforcement, border patrol, forecasting, target acquisition, disaster management, and media. These vehicles all have substan- tial limitations in terms of flight endurance. In recent years soaring flight has sparked interest in UAV applications where significant efforts have been focused on improving energydensitytomaximizeflightenduranceandrange. Soaring is generally defined as the art of keeping a motorless aircraft airborne, or sustainedinflightbeyondsimplysinkinginstillair,usingonlytheaircurrentsthatoccur in the atmosphere. [Short, 2005] Glider pilots employ soaring as it allows the glider to gain altitude with no expenditure of an onboard power supply. Upon investigation of UAVs in soaring applications, flight tests have verified the potential for considerable endurance,rangeandcross-countryspeedimprovements. Takingadvantageofanalogiesbetweentheaerodynamicdesignofbirdsandaircraft isa popular approach in the literature. A summaryofwork priortowind tunnelstudies onlivebirdsispresentedin[Parrott,1970]wherethemainconcernistheaerodynamics of nonflapping flight without acceleration called equilibrium gliding. A wind tunnel studyforglidingflightcanbefoundin[RosenandHedenstrom,2001]. UAVs differ in crucial ways from transport aircraft, and these differences might enable them to better exploit the wind energy. In general, UAVs are lighter in weight which makes this energy source potentially more effective. Because there is no human 1 on board, UAVs can perform more drastic maneuvers than manned aircraft in flying wind-optimaltrajectories.[QiandZhao,2005] The overall performance of gliders in thermal soaring applications is based on two separate missions, namely the gliding phase that utilizes energy and the thermal climb- ing phase that extracts energy using rising air currents. While soaring thermals the aircraftstoresenergyinthepotentialformwhichcanbeefficientlyusedduringitsglide. Inordertoensurethattheaircrafthasthepotentialtomeettheperformancerequire- ments, the particular aerodynamic demands involved must be taken into account for the aircraft design. The relevant compromises are reported in [Cone, 1964]. Once the designiscompletedandthepartsoftheaircraftareassembled,carefullychosensoaring strategies prove an additional source of flight performance enhancements which in turn provide further feedback for designers. As such, there has been considerable interest in moderngliderdesignandsoaringflightduringthelastfewdecades. ThesoaringtechniquesfortheproperlydesignedUAVscurrentlyrequiresubstantial human intervention. The lack of efficient thermal utilization strategies decreases the autonomy of these vehicles. We consider the effects of the rising air currents on the dynamics of the aircraft and define mission specific trajectory generation problems in ordertodeveloppathplanningalgorithmsforautonomousUAVs. The soaring performance of a UAV depends a great deal on the thermal character- istics including the thermal strength and location. Onboard sensor measurements are typically corrupted by noise, whereas gusts and turbulence effects can cause consider- ableperformancedegradationduringinter-thermalglide. We consider the performance losses due to incorrect thermal data and the effect of stochastic gust disturbances. In the presence of additional sensor inaccuracies the deteriorationinthesystemresponseissimulated. Therecoveryoflossesduetoincorrect thermaldataispossibleonlyifdataresourcesbeingusedcouldbeenhanced. However, 2 an adaptive tracking control implied by a Linear Quadratic Regulator (LQR) design inherently provides the optimal control when the aircraft is subject to gust represented byaGaussianwhitenoise. The LQR design needs to be improved if in addition stochastic noise is posed at each sensor measurement channel. By including an adaptive Kalman-Bucy filter and modifying the adaptive law accordingly, we implement an adaptive Linear Quadratic Gaussian(LQG)controldesign. Thesystemresponseisshowntomeettheperformance requirementsinthepresenceofstochasticprocessandmeasurementnoise. Variationsintheclimbrateofatmosphericaircurrentsandthelimiteddurationperi- odsforwhichthesecurrentsareavailableforsoaringflightposeextraconstraintsonthe thermalsoaringproblem. Iftheseconstraintsarenoteffectivelyaddressed,theresultant performancelossofthevehiclemightoffsettheadvantagesofemployingsoaringflight. The thermals are most commonly detected at locations outside the original flight path. Itisthereforenotstraightforwardhowtoformulatethesoaringproblemsothatthe overallflightperformancecanbeimproved. Inanefforttoresolvethisissuewepropose that the optimal static soaring problem can be efficiently adapted to a VRP, a generic combinatorial optimization problem well-known in academic literature. This problem formulationprovescapableofbeingconvenientlyextendedtotheVRPwithTimeWin- dows (VRPTW) considering duration constraints such as limited thermal lifespan. The variations in the air currents can also be readily investigated using the formulation for theDynamicVehicleRoutingProblem(DVRP). Based on our VRPTW formulation of thermal soaring flight, we develop an exact solution method including Preprocessing, Route Optimization and Route Validation. The optimal flight route along with the optimal set of horizontal flight velocities for corresponding inter-thermal flight arcs are determined, and the total flight time mini- mizationisachievedwhilecertaintemporalcriteriaareessentiallysatisfied. Theability 3 toeffectivelyplantheflightpathisthusimproved,andconsiderableincreaseinthelevel of autonomy of a soaring UAV is attained. We present a simple illustrative example where the simulation results approve that the VRP approach motivates further research inthediscussionofoptimalsoaringflightforUAVs. For the discussion of optimal flight across dense thermal regions the use of exact soaring algorithms might appear to be computationally cumbersome. We present a rel- evantproblemscenariowheretheobjectiveistoclimbtheassignedthermalsinthearea of interest and complete the flight mission in minimum time. An approximate solution methodologyisdevelopedbasedondividingthemaneuveringareaintomainsectorsand applyingaminimalspanningtreealgorithmtocoverthesetofthermalsdetectedineach sector. AParallelSavingsBasedHeuristicisincludedimprovingthepathdecisionpro- cesswhileincorporatingamaximumdistanceconstraint. Thesimulationsarepresented inordertoverifytheperformanceoftheproposednear-optimalsoaringstrategy. Wefurtherdeviseanear-optimalsoaringalgorithmfortheShortestPath(SP)routing ofUAVsusingGeneticAlgorithms(GAs)basedonthenaturalevolutionparadigm. The thermal locations are defined as the genes on the chromosomes describing the feasible flightroutesforthevehicle. TheflightpathproposedbytheGAiscomparedtothesolu- tion defined by the exact method developed for thermal soaring flight. The simulation results confirm the proposed route to be coinciding with the optimal one obtained by iteratingtheFloyd-WarshallAlgorithmonthethermalmap. ThesoaringperformanceofaUAVisafunctionofboththestrategyusedtogenerate thetrajectoryofthevehicleandthecontroldesignemployedtotracktherespectiveflight path. Recentadvancesintheareaofstaticsoaringassumeknownlineargliderdynamics and no actuator saturation limits. In practice the dynamics of the glider change with flightandenvironmentalconditions,andtheactuatorsmovingthecontrolsurfaceshave 4 mechanical limits. We consider the optimal flight problem in the presence of actuator saturationandlargeparametricuncertaintiesinthedynamicsofthevehicle. There are essentially two remedies for saturating actuators. The input constraints areconsideredinthecontroldesignphasewithseverelimitationsontheachievableper- formance, or the system is augmented with linear conditioning. If the latter approach is chosen, the system performance is not affected when the actuators are operating linearly, but it is modified during and following a saturation event to ensure stability and eventual escape from saturation restoring the intended linear behavior reasonably quickly.[WestonandPostlethwaite,2000] A problem closely linked to the anti-windup compensation is the bumpless transfer problemwhereasupervisingsystemoverseesmultiplecontrollersdesignedforthesame linear system and switches among the controllers avoiding large transients due to the incompatible initial conditions of the controller connected to the plant and the plant itself.[ZaccarianandTeel,2002] Traditional robust control theory furnishes techniques to deal with modeling uncer- tainties provided they are somehow small, but for larger uncertainties it is frequently advisable to use some sort of adaptation to close the control loop. [Pait and Kassab, 1997] Along with the advances in control theory, the content of adaptive control is also growing. Stochasticadaptivesystems[ChenandGuo,1991,DuncanandPasik-Duncan, 1992,Ren and Kumar, 1994], adaptive neural network control [Ge et al., 1998,Ge et al., 2001], the reduced order adaptive control issues [Ioannou and Kokotovic, 1983], adaptive control of time-varying plants [Tsakalis and Ioannou, 1993], adaptive control design for asymptotic tracking and disturbance attenuation for nonlinear systems [Pan and Basar, 1998], adaptive control for input-output models [Khalil, 1996], and as a dif- ferent approach adaptive unfalsified controllers [Safonov and Tsao, 1995,Wang et al., 2006] are discussed. Several effective adaptive control schemes that can accommodate 5 actuatorfailureswithbothsystemparameterandactuatorfailureuncertaintiesaredevel- oped.[Taoetal.,2004] Output regulation of linear systems with input amplitude and rate saturations are discussed using nonadaptive controllers in [Saberi et al., 1999]. The research in robust and constrained control is summarized in [Tarbouriech and Garcia, 1997]. Adaptive controlapplicationsforsystemswithactuatorandsensornonlinearitiesareinvestigated andanadaptiveinversecontrolapproachisdevelopedin[TaoandKokotovic,1996]for commonnonlinearitiessuchasdead-zone,backlashandhysteresis. Althoughsaturationisoneothertypicalnonlinearity,relevantdiscussionintheadap- tive control context is rather immature. In this research we specify the main shortcom- ings of the current anti-windup design methodologies requiring accurate models and bring these methods due to their wide applicability to the attention of the adaptive con- trolcommunity. We develop an adaptive anti-windup compensator to complement a rather well- formulatedadaptivecontrollersothatthesaturationnonlinearitiescanalsobeaddressed for models with unknown plant parameters. We rely upon the adaptive control the- ory and the stability properties offered by the Linear Matrix Inequality (LMI) based anti-windup compensation and construct a combination of both. This approach raises questions about how the stability of the overall system can be established and thus how onecansurmounttheobstaclesthatimpedefurtherdevelopmentsintheadaptivecontrol theoryforinputconstrainedsystems. If the desired linear performance can be recovered for nominal systems subject to input saturation by implementing a Linear Quadratic (LQ) control augmented with an anti-windup compensator, our analysis on its adaptive counterpart reveals the fact that theperturbationtermsduetoplantparametererrorsintheadaptiveschemedonotcause 6 any unbounded signals in the closed-loop, and the system remains stable. Continuous- time lumped Hurwitz systems with input saturation nonlinearities and unknown plant parameterscaneffectivelybenefitfromtheresultingcontrolstructure. The indirect adaptive control design with anti-windup augmentation for stable sys- tems with input magnitude saturation nonlinearities is also extended to the control of plants with input rate constraints. The adaptive control design can be employed to recover the tracking performance of systems with unknown plant parameters by sup- pressingtheadverseeffectsofdominantinputratesaturation. We expect the proposed trajectory generation methodology to be a useful tool in improving flight range and endurance of both manned and unmanned aircraft as it can also be used to alleviate the piloting tasks for the former. The proposed control design, on the other hand, addresses several issues relating not only to UAV applications but to any practical system involving uncertainties and control input constraints which are either inherent because of the nature of the physical devices employed in the system or intentionallyposedduetoefficiencyconcerns. The terminology and the mathematical notations used are standard and presumably clear from the context. Each chapter includes an introduction for specific background material such as basic concepts and related work, and a brief summary appears in the conclusionsectionofeachchapterhighlightingtheimportanttopicscovered. InChapter 2thetrajectorygenerationmethodologyisdescribed. Chapter3introducestheadaptive control design with anti-windup augmentation. A stochastic approach to the optimal soaringproblemisprovidedinChapter4. Chapter5isdevotedtotheoptimalstaticsoar- ing discussion using the VRPTW, followed by a heuristic search algorithm in Chapter 6 for the navigation of UAVs across dense thermal regions with reduced computational efforts. Chapter7employsGAsforSProutingofautonomousgliders. InChapter8the stabilityanalysisfortheindirectadaptivecontroldesignwithanti-windupcompensation 7 is provided in detail. An adaptive control design for systems with input rate-saturation constraints is developed in Chapter 9. Summary and concluding remarks appear in Chapter 10. Appendices document the application of Pontryagin’s Minimum Principle, Bellman-GronwallLemmaandBarbalat’sLemma. 8 2. OPTIMUMTRAJECTORYGENERATION FORSTATICSOARINGOFUAVS 2.1 Introduction The performance of UAVs is restricted by the energy requirements, mainly on-board storagelimitationsinsizeandweightofbattery/electricsources. Solarorfuelcellpower sources, too, are limited by the size and weight restrictions. An effective method for improvingrange,enduranceandcross-countryspeedisextractingenergybyinteracting withthewind,thesignificanceofwhichisevidentforunpoweredaircraftsuchasgliders. Soaringisabiologicallyinspiredtechnologywhichisaverybasicskillemployedby birds. Today, known as a worldwide sport, soaring is referred to flight in an unpowered heavier-than-air aircraft. A glider is powered by gravity as it always sinks through the still air. It is appropriate for soaring applications since it has a streamlined body and long,narrowwingsgivingitacombinationoflowsinkrateandflatglide. Gliderwings aredesignedsothatperformanceimprovementisrealizedthroughenergysavings. Sev- eralfixed-wingconfigurationstylesarepresentedin[Lazos,2005]. In [Chen, 1981] the basic principle of soaring flight is said to be addressed by the 1883notesofLordRayleigh: “Ifthepropercombinationofatmosphericconditionsand topographicalfeaturesproduceaircurrentswhichriseasfastasorfasterthantheglider sinks (in still air), then the machine will remain aloft or climb.” Two major sources of 9 energyaregivenforsoaringflightin[Short,2005]. Aircurrentshavinganupwardtrend called thermals are useful for static soaring while irregularities in the natural wind are efficientfordynamicsoaring. Examplesofstaticsoaringincludeslopesoaringabovethewindwardslopesofhills, mountain ranges and coasts, as well as thermal soaring, which is widely used by pilots for cross-country flights. It is also mentioned in [Short, 2005] that dynamic soaring is stillnotcompletelyunderstood,andthatthesoaringpilotandhiscraftutilizethevarying velocities and directions of the wind, being able to remain airborne in horizontal or evendescendingaircurrents. Onemayreferto[Hendriks,1974,Chen,1981,Boslough, 2002,Sachs and da Costa, 2003,Zhao and Qi, 2004,Zhao, 2004,Wharington, 2004, Parle, 2004,Yang and Zhao, 2005,Lissaman, 2005,Short, 2005] for details on dynamic soaring. Details on static soaring can be found in [Metzger and Hedrick, 1974,Chen, 1981,Reichmann, 1993,Wharington and Herszberg, 1998,Short, 2005]. It is reported in[Allen,2005]thattheenduranceofasmallgliderUAVcanbeextendedfrom2hours to8hoursinwinterandto14hoursinsummerbyusingsoaringflight. Glider aerotow control, which should consider tug, glider and the tow rope dynam- ical models, is presented in [de Matteis and Menichetti, 2000,Milne, 2001], dynamic soaring control problem is considered in [Sachs and da Costa, 2003,Zhao, 2004,Whar- ington, 2004,Yang and Zhao, 2005,Kyle et al., 2005], thermal centering in static soar- ing is discussed in [Reichmann, 1993,WharingtonandHerszberg, 1998], andalanding approachisthemainsubjectof[PiersonandChen,1979]. Ourfocusisonthecontroldesignforstaticsoaring,mainlyforthermalsoaringwhen thegliderisinflightafterithasbeentowed. Thenatureofslope-triggeredthermalsthat riseobliquelyisexplainedin[Reichmann,1993]. Itismentionedthatsuddenincreasein windspeedatthesummitcanpossiblycausethethermaltoleandownwindasitascends further,whereasonlyverticallyrisingthermalsareconsideredinourwork. 10 The basic method of thermal soaring is to use rising currents of warm air as they are detected. Thermals or updrafts, buoyant plumes of air in the lower atmosphere, can rise very rapidly, allowing the glider to attain substantial increase in altitude. Although a glider continuously descends as it glides in still air, it climbs if it is flown through air which rises faster than its rate of sink. Glide ratio, the analog of the lift-to-drag ratio(L/D)inpoweredflight,istheratioofforwardairspeed(V x )torelevantsinkspeed (−V z )instillair. Itcorrespondstotheratioofthedistanceflowntothelaunchheight. A glideratioof 20:1 forexamplemeansthattheglidertravelsforward 20 feetforevery foot it descends. The performance of a glider can be characterized by its polar curve, a graphthatshowsitssinkrate, V z atdifferentairspeeds, V x ,againanalogoustothedrag polar in powered flight. As discussed in [Maughmer, 2003], a speed polar in which the sinkratedoesnotchangetoorapidlywithchangesinairspeedisdesired. Theaircraftmodelusedinourworkisasmallelectric-poweredUAV.Themodelhas been appropriately modified so thattheUAV isconsideredto bean unpoweredaircraft. The aircraft model and the control objective are presented in Section 2.2. By making several assumptions based on gliders and the nature of thermals, an optimization prob- lem is defined in Section 2.3. In the same section the optimum trajectory is generated by a decision algorithm developed to be used for autonomy. Conclusions are drawn in Section2.4. 2.2 AircraftModelandControlObjective Aircraftdynamicsareusuallymodifiedtoaccountforatmosphericdisturbancessuchas winds,gustsandturbulence,[Nelson,1998]andarewritteninthestate-spaceformas: ˙ x =Ax+Bu+C d d (2.1) 11 where x,u,and d arethestatevector,thecontrolinput,andthegustdisturbanceterm respectively. The C d matrixincorporatesthedisturbance, d intotheaircraftdynamics. Inourworkthethermaleffectismodeledastheatmosphericdisturbanceinthelinear dynamics. Just like atmospheric turbulence, which is a random phenomenon generally discussed in a statistical way defining a Markov atmosphere model, encountered lift at a specific thermal location can also be regarded as random and be modeled using a stochasticapproach. Ontheotherhand,theUAVcanusemeasuredaircraftaccelerations as well as pressure-altitude changes to detect thermals once it arrives at the thermal locationsasmentionedin[Wharington,1998]. Therearealsoremotedetectionmethods that are discussed in [Derr and Little, 1970,Noonkester, 1976,Chadwick and Gossard, 1983,Ito et al., 1989,Abreu et al., 1992,Shannon et al., 2002]. We assume that the thermalsareremotelydetected,anddonotdiscussthespecificsofthedetectionmethods that can be found in these references. Our simulations are carried out with thermals of knownstrengths. Aircraft velocity (V), angle of attack (α), pitch rate (q), pitch angle (θ), and alti- tude(h)areusedtodefinethelongitudinaldynamics. Thedecoupledlateral/directional dynamics of the aircraft are defined through yaw angle (ψ), bank angle (φ), sideslip angle(β),rollrateandyawrate. The angle of attack (α), the sideslip angle (β) and the velocity (V) are defined in termsofroll,pitchandyawvelocitycomponentsasfollows: α =tan −1 ( ˜ ω ˜ u ) (2.2) β =sin −1 ( ˜ v V ) (2.3) V = (˜ u 2 + ˜ v 2 + ˜ ω 2 ) 1/2 (2.4) where(˜ u,˜ v,˜ ω)isthesetofvelocitycomponentsinabody-fixedcoordinate. 12 ThevelocityoftheaircraftcanbewrittenintermsofEuleranglesand(˜ u,˜ v,˜ ω)as: V x V y V z = cθcψ sφsθcψ−cφsψ cφsθcψ+sφsψ cθsψ sφsθsψ+cφcψ cφsθsψ−sφcψ −cθ sφcθ cφsθ × ˜ u ˜ v ˜ ω (2.5) where c and s standforcosineandsinefunctionsrespectively. Using(2.2,2.3,2.4,2.5)wederivethevelocitycomponents: V x V y V z = cθcψ sφsθcψ−cφsψ cφsθcψ+sφsψ cθsψ sφsθsψ+cφcψ cφsθsψ−sφcψ −cθ sφcθ cφsθ × cβcα sβ cβsα ×V (2.6) The lateral controls are unimportant in almost all flight conditions except when the glider is soaring thermals at a speed close to stall speed with limitations on sharp turn ratesduetosignificantsink. Thelateral/directionalstatesarezeroatthetrimconditions, andwefocusonthehorizontalvelocitycomponentdefinedas: V x =V x0 +δV x (2.7) whereV x0 isthehorizontaltrimmedvelocity,andδV x isthedeviation. Weapproximate δV x asalinearfunctionofthelongitudinalmodel’sincrementalstates: x s = [δV δα δq δθ δh] T (2.8) Thelinear δV x equationcanthenbeusedtomeasuretheglider’sperformance,and δV x canbereferredtoastheperformanceoutput. 13 Thestatesareavailableformeasurementandaredefinedastheoutputofthesystem throughtheidentityC matrix. Thelinearincrementalstatespacemodelforlongitudinal dynamicsincludingtheperformanceoutputinthesystemdescriptionisgivenas: ˙ x s =Ax s +Bu+Gω +Hδ T (2.9) y =Cx s z p =C p x s where A,B,C,G,H :system matrices C p :performance output matrix x s :incrementalstatevectorforlongitudinal aircraft model, [δV δα δq δθ δh] T u :elevator deflection ω :thermal strengthfunction δ T :negative throttle setting z p :performance output, δV x In (2.9) the control input, u is the elevator deflection, and the influence of the thermal’sverticalvelocity, ω isincorporatedintotheaircraftdynamicsthroughthe G matrix. The original UAV model which is a powered aircraft, is modified by applying a constant disturbance, Hδ T for the entire flight with a negative throttle setting of δ T =−0.3 tocanceltheeffectofthethrustinputforthetrimconditionsandthusobtain azero-thrust(gliding)aircraftmodel. 14 For the particular set of trim conditions given as [ V 0 α 0 q 0 θ 0 h 0 ] T = [40.00004 0.06308 0 0.00994 4000] T ,thenominalsystemmatricesare A = −0.08029 0.44721 0 −0.560010.00005 −2.30141 −6.87585 0.84355 0.04017 0.00141 0.00002 −32.58844−7.80618−0.00007 0 0 0 1 0 0 −0.04997 −0.69718 0 0.69718 0 B= h 0.01158 −0.69954 −37.47948 0 0 i T G= h 0.63434 −9.45270 −44.06030 0 0 i T H= h 4.0986 0 0 0 0 i T andtheperformanceoutputmatrixis C p = h 1.0000 −0.0006 0 0.0006 0 i Theelementsofthe G matrixwoulddeviatefromthenominalvaluesifthethermal isrisingobliquely,andthesystemmatriceswhichdefinethegliderdynamicsdependon changingflightandenvironmentalconditions. If we run a series of simulations in the absence of thermal effects using the aircraft model in glide flight at a set of different horizontal velocities, the respective vertical velocitiescanbetabulatedastheaircraft’ssinkrates. A quadratic curve-fitting application using the (V x ,V z ) pairs as data points renders theexpressionforthepolarcurveas: 15 20 40 60 80 100 120 140 160 180 200 −18 −16 −14 −12 −10 −8 −6 −4 V x (ft/sec) V z (ft/sec) ( V x , V z ) pairs obtained via simulations Quadratic approximation to the polar curve Figure2.1: Polarcurveinterpolation(Section2.2) V z =aV 2 x +bV x +c (2.10) a =−0.000391, b = 0.026, c =−6.3 whichisacloseapproximationtotheperformanceofthegliderintherangeofinterest. TheinterpolationisgraphicallypresentedinFigure2.1. The glider starts its flight from the current location, and is assigned to reach a spe- cific destination point. It utilizes the thermals on its flight path as they are detected. The objective is to reach the destination safely in minimum time which corresponds to optimizing the glider’s performance. We analyze the longitudinal motion aiming to synthesize the control input that guarantees tracking any suggested horizontal velocity reference, V opt intheallowedrange. Accordingtotheclassicaltheoryofcrosscountrysoaringflight,notonlyshouldthe glider climb in discrete thermals, but it should also fly at the optimal speed during the inter-thermal glide. An efficient method to generate the corresponding flight trajectory isdevelopedinthenextsection. 16 2.3 OptimumTrajectoryGeneration Thedetailsofenteringathermal,centering,flightinthethermal,andleavingthethermal are discussed in [Reichmann, 1993,Wharington and Herszberg, 1998]. In our work detailsrelevanttothelateralmotionareomittedasboththeavailabilityandfunctionality ofthelateralcontrolsareassumed. Theverticalvelocityofthethermal,whichisreferredtoasitsstrength,isingeneral a slowly varying function of time. We assume it to be time-invariant, at least while the gliderissoaringwithinthethermal. Thisapproach,asartificialasitmayappear,isalso supportedbythetime-invariantmodelofthethermalprofileconsideredin[Wharington andHerszberg,1998]: ω(D) =We (−D/R) 2 (2.11) where W is the vertical movement, R is the characteristic radius of the thermal core representingthesizeofthethermal,and D isthelateraldistancetothethermalcenter. A thermal profile (W =15ft/sec and R=1ft), and the glider’s expected optimum flightpathareillustratedinFigure2.2. It is possible to conduct a detailed statistical analysis using the Global Positioning System (GPS) traces to obtain more accurate thermal models for a site and season as presentedin[Cochrane,1999]. Theinterestedreadermayalsoreferto[Allen,2006]for thederivationofthermalmodelsusingballoonandsurfacemeasurementdata. Withadialtypevariometer,asdiscussedin[Piggott,2002],aMacCreadyringcanbe mountedroundthefaceoftheinstrumentandcalibratedinspeedstofly. Thefunctionof this movable speed ring developed by MacCready is explained in [Reichmann, 1993], and some applications are briefly illustrated. If the pilot sets the ring at the expected climbrate,itcommandsthecorrectspeed-to-fly. 17 Figure2.2: Thermalprofile:W=15ft/sec, R=1ft(Section2.3) We assume that the glider UAV has a map of the flight region which stores position and vertical velocity data for the next thermal on the flight path, fairly accurately using sensor measurements from onboard and ground-based sources. This map is referred to as the Cognitive Map here since the glider receives environmental information through the map to form a basis for decisions. A one-step-look-ahead policy is applied and the datais processedto evaluate‘glide-and-thermal’or‘glide-and-bypass’optionstomake a decision. Thermals that are too weak are bypassed since they are not worth climbing. On the other hand, if the vertical velocity threshold for thermals used in the decision logic is too high, the glider never switches into thermal soaring mode. During the first flight tests by NASA Dryden Flight Research Center due to the ineffectiveness of the logicusedtodetermineifanupdraftwasstrongenoughtosoar,theglidermissedalmost allthethermals. When the decision made is in favor of thermal soaring, the same theory applied in the design of the speed ring for soaring pilots can be employed to determine the speed-to-fly for the glider UAV. If the decision is to bypass the thermal, the Cognitive 18 Map updates the data for the next thermal in turn, and the decision process starts over again. Assoaringatathermallocationiscompleted,itischeckedwhethertheglidercan switch safely to the maximum allowed speed for the rest of the flight and reach its final destinationinpureglidemodewithouttheneedforanearlylandingorcrashlanding. When the glider climbs inside a thermal, altitude gain is measured. We require that thealtitudelossintheinter-thermalglideisentirelyregainedbeforeleavingthethermal. As soon as this local constraint is satisfied, the glider resumes gliding at the speed the decisionalgorithmassignsbasedonthenextthermaldata. In the soaring problem defined in [Metzger and Hedrick, 1974] a constraint is imposed such that the altitude at the destination is the same as the initial altitude. We impose the same constraint at each thermal being utilized. In [Metzger and Hedrick, 1974] only one thermal is considered which is located at the arrival point. Instead of a single thermal located at the destination appearing as a must to soar, in our approach a decision is first made for each thermal in the flight direction on whether or not to enter thethermal. Ateachthermalbeingutilized,theconstrainttorecoverthealtitudelossin the inter-thermal glide after leaving the previous thermal is then applied. The decision processissummarizedasaflowchartinFigure2.3. The global optimal solution can be obtained if the entire flight path is considered byremovingtheinter-thermalaltitudelossrecoveryconstraintsatthethermallocations. Intuitively, it might be more efficient to make longer use of a strong thermal to climb higher before leaving in search for the next one. It might also be preferable to leave a relatively weak thermal sooner in hope of finding a stronger one later on. However, it is with these local constraints that one can break the original time-optimal problem into a sequence of smaller ones making the optimization problem tractable. Notably, according to the Visual Flight Rules, pilots of gliders are not permitted to continue climbing at the cloud base. The UAV does not have a pilot, but there are cases where 19 Figure2.3: Decide-To-ThermalAlgorithm(Section2.3) a pilot is remotely flying the airplane using vision sensors fitted on the UAV. Even if the UAV is in autonomous flight, the glider needs data from onboard or ground-based measurementsforthespeed-to-flydecisioninapproachingthenextthermal. Atorabove the cloud base, sensor accuracy or communication capability with the ground might be limited, and such flight is not encouraged for that very reason. This implies that the constraint to leave a thermal at some point, even if it is highly efficient to climb it for longertime,hasapracticalinterpretation. Using our algorithm the next thermal under consideration is utilized if W≥W TH holds where W TH is a fixed parameter. If more than one step look-ahead is allowed in terms of the thermal data, it would be advantageous to adjust the vertical velocity threshold accordingly, allowing the glider to skip the thermals on the flight path which are strong but comparably weaker than the ones located possibly further away. This wouldresultinaflightpathwithanoverallreducedflighttime. 20 Theoptimalvelocitytoapproachadetectedthermalcanbefoundbysolvingacon- strainedoptimizationproblemfor t total =t glide +t thermal (2.12) over V 0 def = {V x ∈R | V min ≤V x ≤V max } asfollows: t f = min Vx∈V 0 t total (2.13) subjectto t glide = Z x(t b ) x(t 0 ) (1/V x )dx (2.14) t b =t 0 +t glide (2.15) h(t f ) =h(t 0 ) (2.16) ⇒ t thermal = h(t 0 )−h(t b ) W (2.17) ˙ x s =Ax s +Bu+Gω +Hδ T (2.18) ω(t) = 0, t<t b W, t≥t b (2.19) where x isthedistancecoveredinflightdirection, x(t 0 ) istheinitialposition, t b isthe time to reach the thermal, x(t b ) is the thermal location, h(t b ) is the altitude at which the UAV arrives at the thermal location, W is the constant value approximation to the thermalstrength. We define in (2.12) an allowed interval, V min ≤V x ≤V max for the velocity in the flightdirection. Thevalueof V max representsthespeedabovewhichtheflightbecomes unsafe due to unstable oscillations of the control surfaces and thus should be avoided, 21 whereastheminimumglidespeed, V min issettobeclosetothestallspeedatwhichthe flowsoverthewingsseparateanddonotprovidesufficientlift,andthatthegliderstalls. For a more complex situation where the glider must fly over high trees, hills or low mountains, extra flight safety constraints should be included in order to increase situational awareness. The polar curve approximation can then be used to predict the correspondingaltitudelossasthegliderfliesatanyhorizontalvelocitytowardsthedes- tination. Ifthelowestallowedflightaltitudeisdefinedapriori,achosen V x canappear tobesafeorunsafe. However,iftheallowedaltitudelossistoorestrictiveforaspecific destination,a V x ∈V 0 satisfyingthesafetyconstraintmaynotexist. Insuchasituation, thermalsoaringbeforereachingthedestinationisrequiredwithnootheralternative. SinceweuseaUAVmodelwhichisoriginallypowered,onecanconsiderthesitua- tionwherethepowerisbroughtintoprovidethrust. ThespecificUAVmodelusedhere can possibly be powered when the aircraft is flying close to the lower altitude thresh- old. The option of having thrust available when needed would change the character of the optimization problem and improve the flight performance significantly. With that optiontheaircraftcouldflyfasterandrisklosingmorealtitudeevenwhenastrongther- mal ahead is not expected since it could use power if the limit of altitude loss tolerance is reached before arriving at the next thermal. In that case the cost function needs to be modified to assign a penalty for fuel consumption or to limit the use of thrust under certain conditions. An alternative is discussed in [MacCready, 1998] where the regen- erativebattery-augmentedsoaringispresentedinthesoaringperformanceoptimization problem. Similarly,anewstrategycouldbedevisedtomakethemosteffectiveuseofthe battery, and battery recharge could become advantageous when the glider is precluded fromclimbingbecauseithasreachedthecloudbase. Thesecasesarenotdiscussedany furtherhere. 22 For a specific glider the vertical velocity, V z is dependent mainly on the horizontal velocity,V x intermsoftheglideratio. Usingthequadraticapproximationtothisrelation weobtainasimplifiedproblemstatementforeachglide-and-thermalinterval: V opt = arg min Vx∈V 0 t total (2.20) subjectto t total =t glide +t thermal (2.21) = x(t b )−x(t 0 ) V x + h loss (V x ,V z ) W where h loss (V x ,V z ) =−{x(t b )−x(t 0 )}× V z V x (2.22) V z ∼ =aV 2 x +bV x +c (2.23) In the above formulation the glider dynamics are used implicitly to obtain the approximatepolarcurverelationandrepresent V z intermsof V x . The total distance to be covered is only a factor to be considered in the Decide- To-Thermal algorithm at the stage where it is decided whether the glider can fly at the maximumallowedspeedwithnothermalsoaringahead. If the decision is made to soar the next thermal, the greater the distance to that thermal location, x(t b )−x(t 0 ), the longer the time spent in the inter-thermal glide. The optimal glide speed itself does not directly depend on the term x(t b )−x(t 0 ) if a feasiblesolutionisassumedtoexistinthefirstplace. Itispossiblethatthenextthermal location is too far away from the current position for the glider to cover the distance, x(t b )−x(t 0 ) for any V x , and therefore it has to land before completing the task. In 23 such a case there are obviously no feasible solutions to the static soaring problem and nosenseindiscussingtheoptimalityofthesolution. The application of Pontryagin’s Minimum Principle (see Appendix A) to the opti- mization problem in (2.20) proves that our problem definition is consistent with the approachusedin[MetzgerandHedrick,1974]. Nonlinear Programming (NLP) describes in general the optimization problems where either the cost function is nonlinear or the constraint set is specified by non- linear equations and inequalities. Although NLP lies within the continuous problem category,severaloptimizationproblemsofhybridcharacterarealsostronglyconnected withNLP.[Bertsekas,1999] The optimization problem we discuss in its simplified form in (2.20) is continuous and is considered as an NLP problem. The relevant optimality conditions are investi- gatednext. GivenaconvexsubsetofRsuchas: V 0 ={V x ∈R :V min ≤V x ≤V max }, (2.24) thefunctionf:V 0 →Rwhere t total =f(V x ) = x(t b )−x(t 0 ) V x + h loss (V x ,V z ) W (2.25) isnotaconvexfunctionof V x sincethecondition f(αx+(1−α)y)≤αf(x)+(1−α)f(y) (2.26) 24 is not true in general ∀x, y ∈ V 0 , ∀α∈ [0, 1]. The mapping becomes a concave functionforsomevaluesof W ,andthusthefirstordernecessaryoptimalitycondition: ∇ Vx f(V opt ) = dt total dV x | Vx=V opt = 0 (2.27) where dt total dV x =− x(t b )−x(t 0 ) V 2 x − x(t b )−x(t 0 ) W V x dVz dVx −V z V 2 x (2.28) isnotsufficientforoptimality. For dt total /dV x tobezero,however,onemusthave − x(t b )−x(t 0 ) V 2 x − x(t b )−x(t 0 ) W V x dVz dVx −V z V 2 x = 0 (2.29) implyingthat ⇒W +V x dV z dV x −V z = 0 (2.30) ⇒W +V x (2aV x +b)−(aV 2 x +bV x +c) = 0 (2.31) ⇒W +aV 2 x −c = 0 (2.32) andthus V opt mustsatisfy V opt = r c−W a > 0, c−W a > 0 (2.33) Thesecondordernecessaryoptimalityconditionrequires ∇ 2 Vx f(V opt ) tobepositive definite: d 2 t total dV 2 x | Vx=V opt > 0 (2.34) 25 where d 2 t total dV 2 x =−{x(t b )−x(t 0 )} (−2) V 3 x − x(t b )−x(t 0 ) W 1 V 4 x × (2.35) dV z dV x +V x dV 2 z dV 2 x − dV z dV x V 2 x − V x dV z dV x −V z (2V x ) =−{x(t b )−x(t 0 )}[(−2/V 3 x )+(1/W)(1/V 4 x )× {[2aV x +b+V x (2a)−(2aV x +b)]V 2 x −[V x (2aV x +b)−(aV 2 x +bV x +c)](2V x )}] =−{x(t b )−x(t 0 )} −2/V 3 x +(1/W)(2c/V 3 x ) =−(x(t b )−x(t 0 ) | {z } >0 (−2/V 3 x ) | {z } <0 [1−(c/W)] ∇ 2 Vx f(V opt ) > 0 is satisfied if and only if 1− c W > 0, which defines a default minimum threshold for the thermal strength that should in general be set at a higher valueinpracticalapplications. Notefurtherthatforasetof m thermalssoaredalongtheflightpath,ifglidingand thermalsoaringtimesareadded,totalflighttimecanbewrittenas: m X i=1 t i total = m X i=1 (t i glide +t i thermal ) (2.36) The total flight time in (2.36) is suboptimal since the optimization is performed locally,i.e.,segment-by-segmentseparatelyforeachcruiseinterval. This trajectory generation algorithm assumes thermals to be available for soaring at the same strength at all times once they are detected. Duration constraints such as limited thermal lifespan and the variations in these air currents are further discussed in Chapter4. 26 2.4 Conclusion In this chapter we give a literature review on soaring flight including the background and the recent advances in soaring methods. We introduce an aircraft model which is a proper modification of a small electric-powered UAV. The basic assumptions of the model are stated, and the corresponding details of the static soaring flight are specified. Formulatingtheflightrequirementsintheformofanoptimizationproblemwhereatmo- spheric thermals are utilized improving the flight range and endurance, we propose an optimum trajectory generation algorithm in order to define the flight path for soaring UAVsandthusincreasetheautonomyinUAVoperations. 27 3. ADAPTIVELINEARQUADRATIC CONTROLDESIGNWITH ANTI-WINDUP AUGMENTATION 3.1 Introduction Oneofthemostcommonmemorylessnonlinearitiesinphysicalsystemmodelsisinput saturation. Saturation characteristics are common in all practical amplifiers (electronic, magnetic,pneumatic,orhydraulic),motors,andotherdevices;theyarealsoused,inten- tionally,aslimiterstorestricttherangeofavariable.[Khalil,2002] Although the final aim is to design an adaptive control which takes the actuator saturationlimitsintoaccount,wesimplychoosetostartourdesignbyinitiallyignoring these limits. The motivation behind is to use the two-step design procedure where the controlfornosaturationlimitsisdesignedfirsttobeaugmentedafterwardswithananti- windup compensator in order to improve the system behavior when the control signal undergoes saturation as discussed in [Turner and Postlethwaite, 2004]. However, we further update the anti-windup compensator parameters along with the parameters of our LQ controller using the model parameter estimates generated by a robust adaptive law.[Kahvecietal.,2007d] 28 The details of the adaptive LQ control design with anti-windup augmentation are presentedinSection3.2. Theefficacyofthecontroldesignisdemonstratedviasimula- tionspresentedinSection3.3,followedbyourconclusionsinSection3.4. 3.2 AdaptiveControlDesignforAircraftSubjectto ActuatorSaturation An LQ control law with disturbance rejection is discussed for the nominal glider UAV model. The controller gains for the LQ design are updated using the model parameter estimates generated by a robust adaptive law. We address the problem of actuator sat- uration limits which is overcome by employing the adaptive version of an LMI based anti-windupcompensator. 3.2.1 LQControlDesignfortheNominalAircraftModel Allsystemstatesareavailableinourproblem,andalthoughwearemainlyinterestedin thehorizontalvelocitycomponentintheflightdirection,itisstillfavorabletoexploitthe availability of system states considering the case as a single-input multi-output system, asubclassofmulti-inputmulti-outputcase. Theclassicalcontroltheoryoflinearsystemsbasedonfrequencyresponseandroot- locustechniques,asdiscussedin[Mooij,1998],arenotsimpletoapplytomulti-variable plants, whereas state feedback systems are particularly suitable. For our design pur- poses,ratherthancomputingthecontrollergainsbyanypoleplacementmethod,anLQ control design based on the relevant mathematical definition of an LQ cost criterion is chosen due to its straightforwardness in systematic implementation. Moreover, if the statemeasurementsavailableforfeedbackarecorruptedbynoise,aKalman-Bucyfilter 29 for optimal state estimation can be employed using the stochastic-deterministic dual- ism,andtheSeparationPrinciplecanbeappliedinordertosolvetheoptimal-estimation problemandthedeterministiccertainty-equivalencecontrolproblemseparately. Having the flexibility to refine the design to include a separate optimal estimator is regarded as anadditionalmotivationforourchoiceofthecontrolscheme. IntheLQcontroldesignweassumethatthesystemmatricesaretime-invariant,and thattheoptimizationintervalfortheperformanceindexapproachesinfinity. Thesteady- stateregulator problemfortheaugmentedsystemisdiscussed, simplifyingthesolution to solving Algebraic Riccati Equation (ARE) instead of Riccati Equation in its more generalform. GivenanLQRproblemthatdoesnotincludeameasureofterminalcontrolaccuracy, if (A,C) is observable, and (A,B) is controllable, a solution to the steady-state LQR problemexists. Theuniquepositivedefinitesolution P totheARE: A T P +PA+Q lq −PBR −1 lq B T P = 0 (3.1) R lq > 0 Q lq =C T C definestheoptimalfeedbackgain, K =R −1 lq B T P (3.2) suchthattheclosed-loopsystemisasymptoticallystable. Note that under weaker conditions of stabilizability and detectability, existence and stabilizationareguaranteed,but P mayonlybepositivesemidefinite. 30 Ourgoalistodesignacontrolinputforthedisturbedsystem ˙ x s =Ax s +Bu+d (3.3) z p =C p x s and ensure that the performance output, z p tracks a given reference signal while the constantdisturbanceterm, d isbeingrejected. We first reduce the tracking problem in (3.3) to an equivalent disturbance rejection problem,andthenapplyaconstantdisturbancerejectionmethodasdiscussedin[Dorato etal.,1995]byaugmentingthesystemproperly. Given the optimal horizontal velocity, V opt , a reference signal, z r is defined to be thecorrespondingdeviationfromthehorizontaltrimmedvelocity, V x0 : z r def =V opt −V x0 (3.4) We can obtain a set of compatible states using the left pseudo inverse of the C p matrixasaCommandStateGenerator(CSG): z r =C p x r ⇒ x r =C T p (C p C T p ) −1 z r (3.5) Thismethodisindeedquitecommoninmodernflightcontrolsystemswhichinvolve multipleredundantsensorstomeasurethesystemstates,directlyorindirectly. Forreferenceinputshaping,asimplefirstorderlow-passprefilterisincludedinour simulationsbutisignoredinthecontroldesignphase. Wedefinethestateerror, e usingthestatevectorandthestatereference: e def = x s −x r (3.6) 31 Foranyconstantstatereference, x r wecanwrite ˙ x s − ˙ x r =Ax s −Ax r +Ax r +Bu+d (3.7) ¯ d def = Ax r +d (3.8) ⇒ ˙ e =Ae+Bu+ ¯ d (3.9) z p =C p x s −C p x r +C p x r (3.10) ⇒z p =C p e+C p x r =C p e+z r (3.11) z new def = z p −z r =C p e (3.12) wherethedisturbance,d actingontheoriginalsystemandthestatereferenceterm,Ax r are combined to form an artificial disturbance, ¯ d acting on the state error dynamics. Asymptotically rejecting ¯ d would also reject d and regulate the state error, implying inturnareferencetrackingresponsefortheoriginalsystem. Torejectartificialdisturbance, ¯ d weaugmentthestatevectorandsolvetheregulator problemfortheaugmentedsystem ˙ x aug =A aug x aug +B aug v (3.13) z new =C aug x aug where x aug = ˙ e z new , v = ˙ u, A aug = A 0 C p 0 , B aug = B 0 , C aug = h 0 0 0 0 0 1 i . 32 IftheAREfortheaugmentedsystemiswrittenas: A T aug P +PA aug +Q aug −PB aug R −1 aug B T aug P = 0 (3.14) Q aug = 0 0 0 I ×Q lq , Q lq > 0 R aug =R lq > 0 with the assumptions that (A aug ,C aug ) is observable and (A aug ,B aug ) is controllable, thecontrollerparameterscanbeobtainedusing K =R −1 aug B T aug P (3.15) Thecontrolinput, v fortheaugmentedsystemis v =−Kx aug =− h K 1 K 2 i x aug =−K 1 ˙ e−K 2 z new (3.16) where K 1 ∈R 1×5 and K 2 ∈R. The control input for the originalsystem in(3.3) isthen foundby integrating(3.16) whichresultsinaPIcontrolleroftheform: u =−K 1 e−K 2 Z t 0 C p e(τ)dτ (3.17) Theaugmentationin(3.13)isnotintendedtominimizetheeffectofanydisturbance ingeneralonthesystemin(3.3),butitratherhelpstoasymptoticallyrejectthenegative throttle setting term acting as a constant additive disturbance which enters in the initial state of the system, and also to track the set of constant reference values that form the desiredresponse. 33 At this point of the design, if one concludes that the system response denies a good behaviorthroughexaminingthesimulationresultsforthetrackingproblem,thetracking propertiescanbeimprovedusinganasymptoticallystablereferencemodel,(A r ,B r ,C r ) with asymptotically stable matrix, A r . The LQ design problem can then be solved for theextendedsystem( ¯ A, ¯ B, ¯ C)asgivenin[Sima,1996]where ¯ A = A 0 0 A r , ¯ B = B 0 , ¯ C = C 0 0 C r . 3.2.2 RobustAdaptiveLaw As the flight conditions change, the system matrices also deviate from their nominal values, and thus we aim to have a controller that can compensate for possible changes intheglidermodelparametersduringflight. Assuch,adaptationinthetrackingcontrol designbecomesessentialtoguaranteerobustnesstomodelingerrorsandtohandlenoise and unexpected changes in the system that result in time varying A and B matrices. Sincemostmultiplicativedisturbances(eitherstatedependentorcontroldependentsuch as control efficiency loss) canbemodeledasperturbationsin A and B matrices, they can also be compensated by the adaptive control. This is the strength adaptive control methodologyaddstoourdesign. Anonlineestimatorisoftenimplementedusingadigitalcomputer. Wethereforeuse the discrete approximation of the continuous-time online parameter estimator based on thecontinuous-timemodelofthesystem. For parameter identification purposes the continuous-time model with unknown A and B matrices ˙ x s =Ax s +Bu+Gω +Hδ T (3.18) 34 isfirstrewrittenwithmeasurabletermscollectedontheleft-handside: ˙ x s −Gω−Hδ T =Ax s +Bu (3.19) To avoid high frequency sensor noise amplification by the derivative term, we filter the differentiated signal, and thus all other terms, by a first order filter 1/(s +λ), introducingthetuningparameter λ> 0 : s s+λ x s − 1 s+λ Gω− 1 s+λ Hδ T | {z } def = z = 1 s+λ Ax s + 1 s+λ Bu (3.20) Foranyfixedflightconditionswethusobtain z = h A B i | {z } def = θ ∗T 1 s+λ x s 1 s+λ u | {z } def = φ (3.21) andtheestimationmodelis ˆ z i =θ T i φ, i = 1,2,...,5 where θ i (t) = h θ i1 θ i2 θ i3 θ i4 θ i5 θ i6 i T (3.22) = h ˆ a i1 ˆ a i2 ˆ a i3 ˆ a i4 ˆ a i5 ˆ b i i T istheestimateof θ ∗ i attime t,and θ ∗ i ∈R 6 isthe i th columnofthematrix θ ∗ ∈R 6×5 . Note also that ˆ a ij is the estimate of a ij where a ij is the ij th element of A∈R 5×5 , and ˆ b i istheestimateof b i where b i isthe i th elementof B∈R 5 . Theregressorvectoris φ = 1 s+λ h x s 1 x s 2 x s 3 x s 4 x s 5 u i T (3.23) 35 where x s i isthe i th elementofthestatevector, x s ∈R 5 . Let k =1,2,... be the discrete-time index, and T be the sampling period where estimationisdoneat t=kT . Thediscrete-timeestimationmodelcanthenbewrittenas: ˆ z i (kT) =θ T i (kT)φ(kT), i = 1,2,...,5 (3.24) FortheadaptivelawweemploythediscreteversionoftheLeastSquaresAlgorithm (LSA)givenin[IoannouandSun,1996]. WeapplyamodifiedLSAwithrobustweight- ingseparatelyforeach i=1,2,...,5 asfollows: P i (kT) =P i ([k−1]T)− c i (kT)P i ([k−1]T)φ(kT)φ T (kT)P i ([k−1]T) m 2 (kT)+c i (kT)φ T (kT)P i ([k−1]T)φ(kT) (3.25) i (kT) = z i (kT)−θ T i ([k−1]T)φ(kT) m 2 (kT) (3.26) θ i (kT) =θ i ([k−1]T)+ p c i (kT)P i (kT)φ(kT) i (kT) (3.27) whichcanalsobewrittenas: θ i (kT) = θ i ([k−1]T)+ c i (kT)P i ([k−1]T)φ(kT)[z i (kT)−θ T i ([k−1]T)φ(kT)] m 2 (kT)+c i (kT)φ T (kT)P i ([k−1]T)φ(kT) (3.28) where c i (kT) = c i1 , φ T (kT)P i ([k−1]T)φ(kT)≥δ i c i2 , otherwise (3.29) 36 istherobustweightingfor c i1 > >c i2 > 0,and δ i > 0. Notethattheestimationerror, i (kT) isnormalizedwhere m 2 (kT) = 1+φ T (kT)φ(kT) (3.30) is a proper normalization signal, and P i (0), c i1 , c i2 , δ i , i = 1,2,...,5 are design parameters. To eliminate possible parameter drift due to modeling errors, estimates of param- eters are forced to lie inside a compact set by using parameter projection so that their boundednessisguaranteed. Itisassumedthatthelowerandupperlimitsofeachelementθ ∗ ij ofθ ∗ i ,j=1,2,...,6, areknownsothat L ij ≤θ ∗ ij ≤U ij forsome L ij , U ij ∈R, i=1,...,5, j=1,2,...,6. Theprojectionsetisacrossproductofthefollowingfiveclosedsetsontherealline: S i ={θ i ∈R 6 | L ij ≤θ ij ≤U ij , j = 1,2,...,6} (3.31) foreach i=1,2,...,5.Theorthogonalterm-by-termprojectionisthengivenby ¯ θ ij (kT)=Pr{θ ij (kT)}= θ ij (kT),L ij ≤θ ij ≤U ij L ij , θ ij (k)<L ij U ij , θ ij (k)>U ij (3.32) where θ ij (kT) istheestimateofthecorrespondingelementofthe θ ∗ matrixat t=kT, and is calculated through our modified LSA with robust weighting. The LQ controller parameters are updated accordingly using the projection ¯ θ ij (kT) in (3.32) instead of θ ij (kT) given in (3.28). Note that since the control design is based on the explicit model such that the model parameters are estimated online and are used to calculate the controller parameters, we implement an indirect adaptive control design in the end. 37 Interestedreadermightreferto[Kahvecietal.,2007a]foranapproximateimplementa- tionoftheadaptiveLQcontrolstructure. The nominal glider UAV model is controllable and observable. However, due to changes in the system parameters, the estimated model might become unobservable or uncontrollable. In addition, even when the estimated model is both controllable and observable,itmightappeartobejusttheoppositebecauseofround-offerrorsinnumer- ical calculations. To prevent the system error in such cases, prior to control updates, lowerlimitconditionsforthesingularvaluesofestimatedlinearsystem’scontrollability and observability matrices are defined. If the minimum singular values for control- lability and observability matrices go below these limits at any discrete-time instant, t=kT,k=1,2,...,andthereforetheestimatedmodelisconcludedtobeuncontrollable or unobservable, the LQ controller gains are chosen relying on the system parameter estimatesat t=(k−1)T . Our adaptation criterion is conservative in the sense that we are making update-or- notdecisionsdependingonthecompletecontrollabilityandobservabilityofthesystem insteadofdirectstabilizabilityanddetectabilitymeasures. Stabilizabilityanddetectabil- ity are indeed sufficient for a solution to the ARE corresponding to the augmented sys- tem to exist. However, stabilizability and detectability check methods, as can be found in [Anderson and Moore, 1990] are rather complex and are not preferred to be imple- mentedforthatreason. As polar curves obtained from corresponding linear models differ slightly for dif- ferent flight conditions, (V x ,V z ) relationship could also be updated using the parameter estimatesoftheglidermodel. Itisclearthatamoreaccurate(V x ,V z )relationshipwould result in a more proper cruise speed decision. The realization, however, would be com- plicated. As the linear model parameters change, the estimated system parameters may not converge to their real values. Moreover, at least three points are needed to obtain 38 anyroughsecondorderpolarcurveapproximationtotheperformanceofalinearmodel, which implies that for the estimated model at least three sets of (V x ,V z ) pairs from the simulationsoftheinitialresponseareneeded. Apossiblewaytoapproximateapolarcurveonlineistohavepolarcurvesforseveral different sets of models tabulated beforehand and to apply an online scheduling. This polarcurve,ontheotherhand,wouldaddanextraloadofcalculationsandmightinvolve largererrors. We use the same simulation-based polar curve in all our applications including the cases for changing flight and environmental conditions. Yet, the results show no sub- stantialneedforfurtherupdates. 3.2.3 AdaptiveAnti-windupCompensatorDesign The actuator dynamics, included as 20/(s + 20) in the simulations, can simply be ignoredduringthecontroldesignduetothenegligibleeffectsobservedandinordernot to increase the order of the controlled system. However, it is not reliable to exclude the upperandloweractuatorsaturationlimitsinthecontroldesignphase. An important practical problem we observe with the aircraft model is the control inputsaturation. Theadmissiblecontrolinputsformasubsetof R definedas: U ={u∈R |− ¯ u≤u≤ ¯ u, ¯ u = 46deg} (3.33) Thisadmissiblecontrolinputsetcanbeincludedintheaircraftmodeldescription: ˙ x s =Ax s +Bsat(u)+Gω +Hδ T (3.34) y =Cx s 39 z p =C p x s sat(u) =sign(u)min(|u|, ¯ u) where the previously proposed solution of employing an adaptive LQ controller does notworkingeneral. In the cheap control problem the control effort is viewed as inexpensive, and thus the weighting matrix, R lq in the LQ design is chosen to approach zero. The control inputinourcaseislimited,butchoosing R lq largeandthusweightingthecontrolinput heavily in the quadratic cost for the LQ problem would only describe the trade-off in thedesignin thesensethatcontroleffortisexpensive. However,itwouldnotexplicitly guaranteethatthecontrolinputwouldremaininbetweenthetwolevels, ± ¯ u. Under challenging circumstances where the system demands unreasonable control usage, allowing R lq to vary exponentially with respect to control inputs that exceed physical limitations is remarkedin[FerrariandStengel, 2004]. Therelevantdiscussion isbeyondourscope. Itisstatedin[Grimmetal.,2003]thatPIorPIDcontrolsaremostseriouslyaffected byinputsaturation,andthestateoftheintegratormightwinduptoexcessivelylargeval- ues. As discussed in [Cao et al., 2002], one approach is to take constraints into account atthecontroldesignstep,whereasasecondoneistofirstignoreactuatorsaturationand then add an anti-windup compensator to the implemented control design to weaken the adverse influence of saturation. Likewise, we have a PI controller in the loop, and we choose to design a proper anti-windup compensator to minimize the performance loss duetotheinputconstraints. Wealsoobservethatifthesaturationavoidanceisaimedincontrolorinputshaping filter design phase, there would be a great loss in the performance due to the inevitably conservativeapproach,primarilybecauseofunknownsystemparametersinourcase. 40 Notethattheanti-windupcompensatorbecomesactiveandleadstoimprovedbehav- ior when the control signal undergoes saturation as discussed in [Turner and Postleth- waite,2004],whereasthesmallsignalresponseofthesystemisleftasintended,andthe local performance of the controller is not restricted. This in turn implies that the track- ing and disturbance rejection features of our LQ control design, which totally ignores theactuatorsaturationlimits,aretobevalidiftheanti-windupcompensatorisdesigned properlytorecoverthelinearperformance. Since we have an uncertain model and an adaptive controller, the anti-windup com- pensator itself should also have an adaptive nature. We propose that the parameters of theintegratedanti-windupcompensatorarealsoupdatedinafashionconsistentwiththe model parameter estimates and the relevant LQ controller parameters when the aircraft is in flight. We choose to update the anti-windup compensator parameters at discrete- timeinstantst=kT,k=1,2,... whereT isthesamplingperiodweusefortheadaptive estimationoftheaircraftparameters. Asaresultwealsocombinetheanti-windupaug- mentationwiththeadaptivelaw. It is shown in [Boyd et al., 1994] that a wide variety of problems arising in system andcontroltheorycanbereducedtoafewstandardconvexorquasiconvexoptimization problems involving LMIs. A quadratically stabilizing anti-windup bumpless transfer compensator synthesis method is proposed in [Mulder et al., 2001] where the relevant synthesis problem is reduced to a convex optimization problem over LMIs. It is dis- cussed in [Grimm et al., 2003] that such a construction method might possibly yield unfeasibleLMIconstraintsinseveralsituations. The anti-windup augmentation synthesis problem with zero-order and plant-order anti-windup compensators are reduced to the feasibility of a set of convex LMIs in [Grimm et al., 2003] where zero-order design is shown to be feasible provided a quasi- common Lyapunov function between the open-loop and unconstrained closed-loop 41 exists, and thus it still might not be feasible at all times. Plant-order linear anti-windup compensation, however, is shown to be always feasible for large enough closed-loop L 2 gain if the plant is asymptotically stable. Due to the extra restriction on the zero- order anti-windup feasibility, we prefer to design and update a full-order anti-windup compensator. The nominal linear aircraft model used is asymptotically stable. However, there is noguaranteethattheestimatedmodelstaysstableasthesystemparameterschangewith time. Ontheotherhand,asymptoticstabilityofthemodelisnotonlyasufficientbutalso the necessary condition for the feasibility of the full-order anti-windup compensation. Weproposethatattimes t=kT whenthefull-orderdesignisfoundtobeinfeasible,the estimationat t=(k−1)T iskeptalongforsimplicity. Analternativemethodwouldbe to use an approximation to the LMI solution to relax the stability condition. The latter, however,isbeyondthescopeofourdiscussion. We construct the proper anti-windup compensator based on the method given in [Turneretal.,2004]. Weexcludethedisturbancetermsintheaircraftmodelthroughthe assumptionthattheyarerejectedbythefeatureoftheLQcontrol. Intheabsenceofthesaturationlimitswecanthensummarizethemodelas: ˙ x s =Ax s +Bu (3.35) y =Cx s wheretheperformanceoutput, z p isalsoomittedfromthesystemdescription. OurLQ controldesignin(3.17)formsatransferfunctionmatrix K 1 +K 2 C p /s from −e = (x r −x s ) (3.36) 42 to u,thecontrolinputtotheactuator. Apossiblerealizationofthecontrolis ˙ x c =A c x c +B c (x r −x s ) (3.37) u =C c x c +D c (x r −x s ) where A c = 0, B c =C p , C c =K 2 , D c =K 1 . Thelongitudinalaircraftmodel,(A,B,C,D)hasthetransferfunctionmatrix: G(s) =C(sI−A) −1 B (3.38) sinceinourcasetheD matrixisgivenas: D = [0 0 0 0 0] T (3.39) Arightcoprimefactorizationof G(s) =N(s)M −1 (s) canbeusedtodefine K aw (s) = M(s)−I N(s) (3.40) asthetransferfunctionmatrixfortheanti-windupcompensatorwithsystemmatrices: A aw =A+BF (3.41) B aw =B (3.42) C aw = F C +DF = F C (3.43) D aw = h 0 D T i T = [0 0 0 0 0 0] T (3.44) where F istobechosensuchthat A+BF isaHurwitzmatrix. 43 Asuitablechoicefor F is F =LQ −1 (3.45) with L and Q definedtobethematricessatisfyingthefollowingsetofLMIs: X BU−L T 0 QC T L T U T B T −L −2U I 0 U 0 I −μI 0 −I CQ T 0 0 −W −1 p 0 L U T −I 0 −W −1 r < 0 (3.46) Q> 0 (3.47) U > 0 (3.48) μ> 0 (3.49) where X =QA T +AQ+L T B T +BL, Q∈R 5×5 , U ∈R, μ∈R,and L∈R 1×5 . Note that W p and W r are positive definite weighting matrices by definition, and their inverses are guaranteed to exist. As discussed in [Turner et al., 2004] if the open loop plant is asymptotically stable, a full order anti-windup compensator always exists, and the above set of LMIs are always feasible guaranteeing well-posedness and quadratic stabilityoftheclosed-loopsystem. In [Turner et al., 2004] it is also discussed that the goal should be to optimize per- formance and robustness together, although there is often a trade-off. In our design W r ∈ R is a scalar for reasonable robustness, and W p ∈ R 5×5 is chosen to obtain fast rise time and to avoid oscillations as much as possible. In [Turner et al., 2004] it is furtherexplainedasanotheradvantagethatsimultaneouslyoptimizingperformanceand 44 robustnessusing W p and W r tendstopreventfastpolesappearinginthecompensator dynamics, whereas fast poles would otherwise require a very high sampling frequency forimplementation. The L 2 performance optimizing solution for decision variables Q, U , μ and L isfoundbyemployingMATLABLMIControlToolboxfortherelevantcalculations. Theanti-windupisthenaugmentedintothesystemas: ˙ x aw =A aw x aw +B aw (u act −sat(u act )) (3.50) y aw =C aw x aw +D aw (u act −sat(u act )) where y aw = h y aw 1 y T aw 2 i T , y aw 1 ∈R, y aw 2 ∈R 5 . Note that u act is the unsaturated control,and sat(u act ) istheactualeffectiveinputtothelongitudinalaircraftdynamics whichiscalledthesaturatedcontrolinputsinceitincludesthesaturationnonlinearity. If we neglect the actuator dynamics, 20/(s + 20) in the control design step, the anti-windup modification term y aw 1 changes the control input, u to the actuator by −y aw 1 anddefines u act as: u act =u−y aw 1 (3.51) whereas y aw 2 modifiestheLQcontrollerin(3.37): ˙ x cm =A c x cm +B c (x r −x s )−B c y aw 2 (3.52) u =C c x cm +D c (x r −x s )−D c y aw 2 with A c ,B c ,C c and D c asdefinedbefore. If the initial anti-windup compensator states are set to be zero, as noted in [Kapoor et al., 1997], the extra dynamics do not affect the original controller until sat(u)6=u. 45 Figure 3.1: Adaptive LQ Control Structure with Anti-windup Augmentation (Section 3.2) Similarly, we keep the original controller initially unmodified. If the initial estimate of the system and the relevant initial LQ design do not result in a feasible set of LMIs for the construction of the initial anti-windup compensator, one might also initialize the anti-windup compensator parameters at zero, or the initial estimate of the model might be modified in a way that would result in a feasible set of LMIs for the initial anti-windupcompensatorconstruction. The final form of our system is summarized in Figure 3.1 where the adaptive LQ control structure and the extra feedback compensation at the control implementation stage, i.e., theanti-windupaugmentationareshown. Oneoftheadaptivelawinputs, u isshiftedto sat(u act ) tocompletetheadaptivecontroldesignforthesystemwithinput constraints. Therequiredstepistomodifytheregressorvectorin(3.21)as: φ = 1 s+λ h x s1 x s2 x s3 x s4 x s5 sat(u act ) i T (3.53) Among different anti-windup design techniques, we choose to follow the method proposed in [Turner et al., 2004] not only because the method prescribes a numerical 46 procedurewhichmakesonlinesynthesispossiblewhenafixedorderlinearanti-windup compensator design is aimed, but stability and a level of performance guarantees are also motivating as we combine the adaptive law with the anti-windup scheme based on theCertaintyEquivalencePrinciple. Fulfilling our expectations we do not observe any significant steady-state perfor- mance loss in our simulations with the adaptive anti-windup compensator in the pres- enceofactuatorsaturationnonlinearitiesaspresentedinthenextsection. 3.3 Simulations The efficacy of the Decide-To-Thermal algorithm is first demonstrated in comparison to the pure glide flight where a known glider model with no actuator saturation limits is considered. The adaptive law is then included in the simulation, and the optimal flight is aimed via thermal soaring with an unknown model and no actuator limits. The third simulation allows us to observe the tremendous performance degradation if the original adaptive LQ control design without the anti-windup augmentation is applied to the aircraft subject to input constraints. In the simulation that follows, our adaptive versionoftheLMIbasedanti-windupcompensatordesignisincludedtogetherwiththe adaptiveLQcontroller. Asanalternativebutnotquiteasatisfactorymethod,theresultof includingaproperfiltertoavoidsaturationisalsoillustrated. Finally,anextraadvantage oftheadaptivecontroldesign,actuatorfailureaccommodation,isexemplified. 3.3.1 VerificationoftheDecide-To-ThermalAlgorithm First the aircraft model is assumed to be known, and a fixed LQ controller with distur- bancerejectionisappliedwithnoadaptation. Sincethesaturationlimitsoftheactuator 47 are not included in these simulations, the results demonstrate solely the significance of thermalsoaringinflight-timeoptimizationincomparisontothepureglideflight. Ascenarioischosenwheretherearetwothermalsofequalstrength, W=10ft/sec located at x=2,000ftand x=3,500ft.Thedestinationisat x=5,000ft,andthe maximum V x allowedisassumedtobe V max =210ft/sec. There is also an allowed altitude loss of 400ft defined, and thus the flight at the given V max is not safe for pure glide. To keep the comparison as truthful as possible, by trial and error it is found that the glider can fly as fast as 92ft/sec to lose exactly 400ft andcovertherequireddistanceof 5,000ft,reachingthedestinationin 55sec bypureglide. Ontheotherhand,ifthegliderutilizesthetwothermalsonitspathusingtheDecide- To-Thermal Algorithm, it can reach the destination at t=47sec. Moreover, the total altitudelossismeasuredtobeonly 336ft,whichisstillfarawayfromthegivensafety limitof 400ft. The horizontal distance covered and the altitude change incremental from the trim position are compared for thermal utilization and pure glide motion in Figure 3.2 and Figure3.3respectively. Note that the glider in pure glide flight would need early landing or crash landing right at the destination point if for example the safety limit for altitude represents the presenceofhillsormountainsontheflightpath,whereasthermalsoaringpromisessafe flight. It should also be pointed out that for thermals rising faster, the performance improvementwouldbeevenlargerthanpresentedhere. Oursimulationresultsshowthatwhenadditionalrestrictionssuchasmaximumalti- tude loss limit are defined for safety, thermal soaring might become superior to straight flight not in reduced altitude loss only, but by all means considered here. The glider finishes at a higher altitude which implies an improved flight range, and it takes even 48 lesstimetocovertherequireddistance. Althoughatafirstglanceutilizationofthermals mightappearasawasteoftime,ourresultsconfirmthatinthelongrunitsavestimeand also helps the UAV to arrive at destination points which might not be possible to reach ifthermalsoaringisnotinvolved. 3.3.2 Adaptive LQ Controller for the Unknown Aircraft Model withNoInputConstraints Inordertodemonstratetheperformanceofouroptimalsoaringalgorithmandtheadap- tive LQ controller for the unknown aircraft model, actuator limits are excluded here as intheprevioussubsection. The glider UAV model is initially unknown, but model parameters are estimated throughtheadaptivelawwiththesamplingperiodchosenas T=0.01sec.Ontheother hand, even when the nominal model is known, adaptation is essential since there are unmodeleddynamicsandchangesinstructuresothatparametervariationsareinevitable duringflight. The simulations are based on a simple environmental setup where there is a single thermalofstrength W=5ft/sec thatislocatedat x=3,000ft. The horizontal velocity of the aircraft, the optimal soaring velocity and the given thermal strength function are illustrated in Figure 3.4, whereas Figure 3.5 displays the appliedcontrolinput. Here and for the rest of the simulations in this chapter, only optimal soaring is dis- cussed with no more comparison to the pure glide flight, assuming time, altitude and range advantages thermal soaring provides in the long run are clear so far through the performancecomparisonintheprevioussubsection. 49 3.3.3 Effects of Saturation Limits on the Optimal Autonomous SoaringFlightProblem The simulations in the first two subsections verify that the Decide-To-Thermal Algo- rithmiscapableofachievingtheoptimalperformance,andthattheoptimalautonomous soaring flight problem can be solved in the presence of parametric uncertainties. How- ever, this is only part of the story since the actuator limits do not allow the full use of thecontrolinputatalltimes. Using the previous simulation scenario for the aircraft model subject to saturation nonlinearities demonstrates that it is quite possible for the system to go unstable if the actuatorshavelimitations. TheresultsareshowninFigure3.6andFigure3.7. The Decide-To-Thermal Algorithm can still be applied to generate the optimum flight trajectory, but the adaptive LQ control is short of manipulating the horizontal velocity of the aircraft. Rather than soaring flight path optimization, stabilization of dynamicsbecomesaseriousissuetobeconsidered. 3.3.4 Optimal Soaring Flight Performance with Adaptive Anti- windupCompensatorDesign If our adaptive version of the LMI based anti-windup compensator design is included in the simulation, the system is shown to be stabilized. The properly chosen W p and W r matrices help avoiding oscillations resulting from the anti-windup compensation. The performance of the aircraft and the control input are presented in Figure 3.8 and Figure 3.9 respectively, whereas the modifications on the control input introduced by the anti-windup compensator can be observed from Figure 3.10 and Figure 3.11. Note thatthe 3 rd andthe 4 th elementsofthevector y aw 2 ∈R 5 arescaledbyafactorof 0.4 forconvenienceinthelatterfigure. 50 Consequently,theintegrationoftheanti-windupcompensatorresultsingainingback theoptimalperformancethatwouldotherwisebelostduetotheactuatorsaturation. The behavior of the closed-loop system closely matches the system response in the absence ofactuatorlimitswhichisthepointattheheartoftheoptimalsoaringproblem. 3.3.5 Flight Performance Using an Alternative Method: Including aFiltertoAvoidSaturation Ifanalternativetotheanti-windupcompensatordesignisdesired,thesimplefirstorder prefilter used in the previous simulations can be shifted to a more sophisticated version such that the input saturation is not likely to be triggered in the first place but avoided in an attempt to stabilize the system behavior. Note that such a trade-off to stay with a simplercontroldesignwithoutanti-windupcompensationgivesextremelyconservative resultsingeneral,orintheworstcaseitmaynotwork. Aprefilterdesignedforaspecific case still does not guarantee that the control input at the actuator would always stay in between the saturation limits. Due to different external conditions or changing model parameters, it is possible that a filter designed in advance does not succeed in actuator saturationavoidance. For filter design purposes we employ a second order normalized prototype Butterworth analog low-pass filter with the unity cutoff frequency transformed to 2π0.006rad/sec.Theresultingfilteris 0.00142/(s 2 +0.05331s+0.00142). Inthissimulationweassumethattherearetwothermalsofidenticalstrength, W= 5ft/sec detected at the locations x = 8,000ft and x = 12,000ft. Our choice of the second order filter is not conservative at all for this specific scenario since it allows the control input to the actuator to get quite close to the upper limit, but it does not hit saturation. It is, however, clear that in applications where the goal of the glider UAV is to optimally soar both strong and weak thermals with a broad set of reference values 51 for its horizontal velocity, a more conservative prefilter would be required resulting in furtherperformancedegradation. Distorting the reference signal using the Butterworth prefilter achieves to avoid sat- urationbutalsopresentspoorperformanceasillustratedinFigure3.12andFigure3.13. The prefilter avoids saturation, the glider climbs at the detected thermal locations, and the altitude lost in the inter-thermal glide is recovered. As the aircraft’s response is too slow in tracking the horizontal velocityreference, however, the glider cannotfully ben- efit from thermal soaring in the long run, and thus the significance of the optimization schemeformulatedtominimizethetotalflighttimeislost. Iftheresultingperformance would not present the advantage of utilizing thermals, thermal soaring inevitably loses itsadvantageoverpureglideflight. Theproposedprefiltersolutionisthusnotanaccept- ablealternative. Theseresultsgiveinsightintothesignificanceofanti-windupcompen- sationforoptimalsoaringofaircraftsubjecttoactuatorsaturationnonlinearities. 3.3.6 ActuatorFailureAccommodation Usingouradaptivecontroldesigntheperformancedegradationduetoaclassofactuator failures can be shown to be mitigated to an almost negligible level even in the event of saturationofactuators.[Kahvecietal.,2008] In addition to the simulation scenario where a thermal of strength W =5ft/sec is detectedat x=3,000ft,wealsoconsideracontrolsurfacelossof 20%fort≥15sec. TheresultsarepresentedinFigure3.14andFigure3.15. The simulation results show that our adaptive controller can also tolerate situations withcontrolsurfacefailures. Acontrolsurfacelosscanbemodeledasadeviationinthe B matrix,anditisdemonstratedthatthecorrespondingparameterchangesarehandled properly in real-time by the adaptive control design. Note that a reasonable command- tracking performance is provided, even without the use of smart actuators which can 52 detect and identify the failure size, or a separate post-failure estimator design, i.e., a faultdetectionandidentificationalgorithmisnotincluded. As the system is observed to adapt quickly to the faulty dynamics, multiple failures intimemightpossiblybeaccommodatedtosomeextent. 3.4 Conclusion There are typically unavoidable uncertainties in aircraft modeling where further devia- tions are inevitably introduced as the flight conditions change. In order to achieve the specified control objective under such circumstances we design a robust adaptive flight controller for the UAV model which represents the linearized longitudinal dynamics of theUAVusedforsoaringexperimentsbyNASADrydenFlightResearchCenter. We consider a complex aircraft model including physical limitations on the actu- ator dynamics where conventional tracking controllers are no more sufficient alone to promise stable and optimal performance. To address actuator saturation we include an LMI based anti-windup augmentation. The adaptive anti-windup compensator design combined with the adaptive LQ controller allows optimal autonomous soaring perfor- manceinthepresenceofactuatorlimitsandunknownaircraftdynamics. Our simulations for various scenarios verify the performance of the decision algo- rithm. Itisobservedthattheproposedcontrolschemedoesnotonlymanagetostabilize thesystem,butitisalsowellcapableofaccommodatingaclassofproblemscommonly facedduringflight,namelycontrolsurfacefailure. 53 0 10 20 30 40 50 60 0 1000 2000 3000 4000 5000 6000 7000 8000 time (sec) distance, x (ft) THERMAL SOARING THERMAL SOARING GLIDING GLIDING GLIDING Pure Gliding Gliding and Thermal Soaring Figure 3.2: Distance covered via thermal soaring and pure glide for the known model caseintheabsenceofsaturationlimits(Section3.3.1) 0 10 20 30 40 50 60 −600 −500 −400 −300 −200 −100 0 100 time (sec) altitude change, δh (ft) THERMAL SOARING THERMAL SOARING GLIDING GLIDING GLIDING Pure Gliding Gliding and Thermal Soaring Figure 3.3: Altitude change via thermal soaring and pure glide for the known model caseintheabsenceofsaturationlimits(Section3.3.1) 54 0 10 20 30 40 50 60 −20 0 20 40 60 80 100 120 140 160 180 200 time (sec) velocity (ft/sec) Optimal Velocity, V opt Horizontal Velocity, V x Thermal Strength, W Figure 3.4: Optimal velocity defining the reference, horizontal velocity of the aircraft, andthermalstrengthasfunctionsoftime(Section3.3.2) 0 10 20 30 40 50 60 −1000 −500 0 500 1000 time (sec) elevator deflection (deg) Figure3.5: Controlinputdirectlyappliedtotheactuator(Section3.3.2) 55 0 10 20 30 40 50 60 −600 −400 −200 0 200 400 600 800 time (sec) velocity (ft/sec) Optimal Velocity, V opt Horizontal Velocity, V x Thermal Strength, W Figure 3.6: Instability of the system response for the saturated actuator case (Section 3.3.3) 0 10 20 30 40 50 60 −1 0 1 2 3 4 5 6 x 10 5 time (sec) elevator deflection (deg) Figure 3.7: Control input sent to the actuator for the unknown model where saturation limitsaretobeapplied(Section3.3.3) 56 0 10 20 30 40 50 60 −20 0 20 40 60 80 100 120 140 160 180 200 time (sec) velocity (ft/sec) Optimal Velocity, V opt Horizontal Velocity, V x Thermal Strength, W Figure 3.8: Horizontal velocity, V x of the unknown model using the adaptive anti- windupaugmentationinthepresenceofsaturation(Section3.3.4) 0 10 20 30 40 50 60 −1000 −500 0 500 1000 1500 time (sec) elevator deflection (deg) Control input modified with anti−windup Effective control input including saturation nonlinearity Figure 3.9: Control input modified for the unknown model to recover optimal soaring performance(Section3.3.4) 57 0 10 20 30 40 50 60 −800 −600 −400 −200 0 200 400 600 800 time (sec) Anti−windup modification: y aw1 Figure3.10: Theanti-windupcompensatormodification: y aw 1 (Section3.3.4) 0 10 20 30 40 50 60 −1000 −800 −600 −400 −200 0 200 400 600 800 time (sec) Anti−windup modification: y aw2 elements y aw2 (1) y aw2 (2) y aw2 (3), scaled y aw2 (4), scaled y aw2 (5) Figure3.11: Theanti-windupcompensatormodification: y aw 2 (Section3.3.4) 58 0 100 200 300 400 500 600 700 −20 0 20 40 60 80 100 120 140 160 180 time (sec) velocity (ft/sec) Optimal Velocity, V opt Filtered V opt Horizontal Velocity, V x Thermal Strength, W Figure3.12: Apropersecondorderprefilterchosentoavoidsaturationbysmoothingthe transitionsinthereferencesignalforthecontroloftheunknownmodel(Section3.3.5) 0 100 200 300 400 500 600 700 −20 −10 0 10 20 30 40 50 time (sec) elevator deflection (deg) Figure 3.13: Control input using a proper prefilter so as to avoid the saturation limits (Section3.3.5) 59 0 10 20 30 40 50 60 −20 0 20 40 60 80 100 120 140 160 180 time (sec) velocity (ft/sec) Optimal Velocity, V opt Horizontal Velocity, V x Thermal Strength, W Figure3.14: Controlefficiencylossby20%att=15secwithnosignificantperformance degradationobserved(Section3.3.6) 0 10 20 30 40 50 60 −1000 −500 0 500 1000 1500 time (sec) elevator deflection (deg) Control input modified with anti−windup Effective control input including saturation nonlinearity Figure3.15: Fortheunknownmodelinthepresenceofsaturationandcontrolefficiency loss,controlinputwithanti-windupmodification(Section3.3.6) 60 4. ASTOCHASTICAPPROACHTO OPTIMALSOARINGPROBLEMAND ROBUSTADAPTIVELQGCONTROL 4.1 Introduction Static soaring applications are modeled in general based on the assumption that flight fromonethermallocationtothenextisincalmair. Theinter-thermalglide,however,is more realistically characterized by random gusts acting on the aircraft dynamics. Gust isbydefinitionastrong,abruptrushofwind,anditusuallyreferstoasinglepulse. Fora continuousprofiletheguststructureisgenerallyspokenofasturbulence.[Hoblit,1988] The atmosphere is in a continuous state of motion, and the wind gusts created by the movement of atmospheric air masses can degrade the performance and flying qual- ities of an airplane, whereas wind shears created by thunderstorm systems have been identified as the major contributor to several airline crashes. [Nelson, 1998] From the aircraft design standpoint, atmospheric turbulence can contribute to aircraft structural fatigue and even structural damage or failure if the turbulence spectrum has larger scales.[Fuller,1995]Airplaneshavinglowwingloadingsareparticularlymorerespon- sive to the influence of vertical gusts. In order to reduce the negative effects of these atmosphericphenomenaontheflightperformance,propermodelingisrequired. 61 A sharp-edged gust encountered by an airplane can be modeled by a step function profile. The atmospheric turbulence experienced by an airplane is modeled in [Houbolt et al., 1964] as a combination of discrete patches of gusts, which is replaced later on in a limiting-process sense by a model which has a continuously variable distribution of Root Mean Square (RMS) gust velocity. Alternatively, idealization of a typical gust profile as a stationary Gaussian random process is in wide use. [Hoblit, 1988] The two reasons discussed are the facts that it is vastly more realistic than the simple discrete- gust idealizations and that the statistical characteristics of the airplane response can be determined directly from the statistical description of the gust velocity profile. The equations of motion in a nonuniform atmosphere are summarized in [Nelson, 1998] wherepureverticalorplungingmotionisdiscussedforsinusoidalorsharp-edgedgusts, and the velocity variations in a turbulent flow are described. To study the influence of atmospheric disturbances on the dynamics of the aircraft, we modify the aircraft model accordingly. A zero-mean Gaussian white noise is used to represent the effect of the verticalgustsonthelongitudinalaircraftdynamics. Aircraft sensors are utilized for required measurements during flight. In-flight cali- bration is not applicable at all times, and the inherent inaccuracies in these sensors are thus to be taken into account in the control design phase. Sensors might suffer accu- racy losses and very slow response times. If the sensor data is collected at periodic rates, this also adds to the uncertainty. The accelerations and the angular rates can be measured and used to obtain the aircraft position, velocity and altitude. These mea- surements are likely to be corrupted by noise, and some degree of robustness against the sensor noise is desired. If the data from redundant or complementary sensors are combined,multi-sensordatafusionpossiblyprovidesmorepreciseinformation. Never- theless, implementing such advanced measurement devices might be both complicated 62 and expensive. In our work maximizing the use of the available information is encour- aged,butnoredundantorcomplementarysensorsareincluded. The recent improvements in remote detection methods for thermals are out of the scope of our work. For our design purposes it suffices to make a quantification of the performancedeteriorationduetotheerroneousthermalstrengthdata. Inthischapterwe donotposethesaturationlimitsattheactuatorsinoursimulations,andthuswediscuss thecontroldesignwithnoanti-windupaugmentation.[Kahvecietal.,2007e] Possible performance losses due to incorrect thermal data and the noise effect of verticalgustsactingontheaircraftdynamicsarediscussedinSection4.2andSection4.3 respectively. Inordertoreducetheperformancedeteriorationduetosensorinaccuracies indicated in Section 4.4, the LQR design is shifted to an LQG controller in Section 4.5 wherethestatisticalpropertiesofthepracticalgustdisturbancesandsensorinaccuracies can be exploited to obtain sharper performance. In Section 4.6 comparative simulation resultsarepresented. ConclusionsappearinSection4.7. 4.2 PerformanceLossDuetoIncorrectThermalData Until rather large errors are involved in the thermal strength data, the optimal reference velocity can be shown to remain almost unaffected. However, the corresponding errors might severely affect the overall performance of the UAV if the optimum trajectory generation algorithm is conservative or aggressive in general, and it might become one of the critical issues since the optimum overall soaring performance is desired. On the other hand, the thermal data should also be used to define the maximum allowed horizontal velocity that meets the safety requirements. This implies the fact that the effectofthethermaldataistwofold. 63 Adiscussiononhowtoachieveimprovementsintheavailablethermaldataresources is omitted here. Instead, we quantify the effect of the data error on the optimal velocity calculation as it might indeed be helpful to estimate the resulting performance loss or therisktakenwiththespeed-to-flydecision. If the thermal strength, W is incorrectlyestimated as W +ΔW e , where ΔW e is theadditiveerror,theoptimalhorizontalvelocitywouldbechosenas: V x | W+ΔWe = r c−W−ΔW e a (4.1) which allows us to calculate the effect of any additive error in the thermal strength data onthecorrespondingflighttime. Thetotalflighttimeisoriginallydefinedas: t total =t glide +t thermal = Δx V x + h loss (V x ,V z ) W (4.2) h loss (V x ,V z )=−Δx× V z V x (4.3) V z ∼ =aV 2 x +bV x +c (4.4) where Δx istheinter-thermalglidedistance. Theerrorinthethermalstrengthdatawouldthenaffect(4.2,4.3,4.4): t totale =t glide e +t thermale = Δx V x | W+ΔWe + h loss (V x | W+ΔWe ,V z | W+ΔWe ) W (4.5) h loss (V x | W+ΔWe ,V z | W+ΔWe )=−Δx× V z | W+ΔWe V x | W+ΔWe (4.6) V z | W+ΔWe ∼ =aV x | 2 W+ΔWe +bV x | W+ΔWe +c (4.7) 64 The relative error in the total flight time, which can be confirmed by simulations, is thencalculatedas: η = t total −t totale t total (4.8) The incorrect data such as smaller or larger strength estimation for the thermals on the flight path can be shown to extend the overall flight time needed to reach the destination. Ifthedistancetothenextthermallocationisestimatedwitherror,thisdoes notdirectlyaffecttheoptimumflightvelocity,butthesafetyrequirementsmightneedto be modified. An estimate of the distance to the next thermal location which is smaller than its real value might increase the allowed maximum velocity and result in an early landingoranaircraftcrash. Ontheotherhand,ifanestimatedvaluewhichislargerthan therealdistanceistakenforgranted,thetrajectorygenerationalgorithmmightconclude that the speed-to-fly is not allowed due to safety considerations and define a reference slowerthantherealoptimumflightvelocityresultinginperformanceloss. Asimulationscenarioischosenwherethreeglidersstarttogether,butwithdifferent strengthdataforthesamethermallocatedatadistanceofΔx=3,000ft.Thefirstglider hasthecorrectthermalstrengthdata, W=10ft/sec.Thesecondgliderunderestimates theclimbrateofthethermalas W=5ft/sec,whereasthethirdglideroverestimatesit tobe W =15ft/sec andfliesaggressively. Afterclimbingatthethermallocationand recoveringtheircorrespondingaltitudelosses,theglidersallflyatthesamespeed-to-fly forthenextthermalaheadwhichhasastrengthof W=5ft/sec. ThesimulationresultsinFigure4.1andFigure4.2revealthatthefirstgliderdisplays thebestperformanceasexpected. Itreachesthethermallocationbylosingmorealtitude than the second glider, but it recovers that altitude loss and still finishes earlier. The second glider flies overcautiously and deviates less from the straight flight path. The third glider arrives at the thermal location earlier than the other two as it flies faster. 65 Meanwhile,itlosesmorealtitudeandthusspendsmoretimeforclimbingatthethermal location in comparison with the other two. The shortest total flight time is observed for thefirstgliderwhichfliesoptimallybasedonthecorrectthermaldata. Using the polar curve approximation the relative error, η defined by (4.8) is calcu- lated to be 0.02 for W =5ft/sec, and 0.01 for W =15ft/sec. The simulation results demonstrate higher relative errors of 0.03 and 0.02 for W = 5ft/sec and W=15ft/sec respectively. The discrepancy is partly due to the quadratic polar curve approximation used to represent the aircraft dynamics in calculations and partly due to the transients observed insimulationswhiletrackingthedesiredtrajectory. Therelativeerrorprovidesvaluable informationtoevaluatetheriskbeingtakenintermsofreliabilityassessment. Moreover, the value of the relative error proposed by (4.8) can be shown to match the simulation resultsmorecloselyiflongerglidedistancesareconsideredwherethetransientsbecome negligibleincomparisontothetotalflighttime. 66 0 5 10 15 20 25 30 35 40 45 50 0 1000 2000 3000 4000 5000 6000 7000 time (sec) distance covered, δx (ft) With Correct Thermal Data With Smaller Thermal Strength Estimation With Larger Thermal Strength Estimation Figure4.1: Theflightdistancecoveredfordifferentthermaldataprovided(Section4.2) 0 5 10 15 20 25 30 35 40 45 50 −600 −500 −400 −300 −200 −100 0 100 time (sec) altitude change, δh (ft) With Correct Thermal Data With Smaller Thermal Strength Estimation With Larger Thermal Strength Estimation Figure4.2: Thecorrespondingaltitudelossfordifferentthermaldataprovided(Section 4.2) 67 4.3 PerformanceLossDuetoVerticalGust DisturbancesinInter-thermalGlide We focus on gust fields that create effective pitching gusts rather than rolling. The gust disturbance,ω g ininter-thermalglidecanthusbeincorporatedintothesystemdynamics in (2.9) through the G matrix similar to the thermal strength function, ω. As the G matrixisindependentofthestatevector, x s andthecontrolinput, u,thismodification defines an additive noise problem. We do not discuss multiplicative disturbances and onlytheadditiveonesasgusteffects. The discussion in [Dorato et al., 1995] proves that for the LQR problems in the presence of additive disturbance signals modeled as zero-mean Gaussian white noise processeswherethesystemstatesareavailableforfeedbackuncorruptedbyanydistur- bance, the validity of certainty equivalence can be established and the optimal control for the stochastic problem can be shown to be identical to the optimal control for the deterministiccasewhere ω g (t)=0 atalltimes. The results in Figure 4.3, Figure 4.4, Figure 4.5 and Figure 4.6 are from simula- tions run in the presence of a gust disturbance, ω g represented as a white noise signal with unity power spectral density. For reference input shaping a simple first order low- pass prefilter is included in our simulations, and the resulting optimal horizontal flight velocityisusedasthereference. AsmoothaircraftresponseisobtainedusingtheLQRdesign,andtheprocessnoise doesnotcauseanyconsiderableeffectonthemotionoftheaircraft. Ifthepowerspectral density of the process noise is chosen to be larger, further deviations from the optimal performancewouldbeinevitable. Still,acceptableperformanceisobservedunlessvery significantgustdisturbancesarepresentasprocessnoise. 68 0 5 10 15 20 25 30 35 40 45 50 0 1000 2000 3000 4000 5000 6000 7000 time (sec) distance covered, δx (ft) For No Noise Case In the Presence of Process Noise Figure4.3: TheflightdistancecoveredusingtheadaptiveLQRinthepresenceofverti- calgust(Section4.3) 0 5 10 15 20 25 30 35 40 45 50 −500 −400 −300 −200 −100 0 100 time (sec) altitude change, δh (ft) For No Noise Case In the Presence of Process Noise Figure 4.4: The corresponding altitude loss using the adaptive LQR in the presence of verticalgustdisturbances(Section4.3) 69 0 5 10 15 20 25 30 35 40 45 50 0 50 100 150 200 250 time (sec) horizontal velocity, V x (ft/sec) For No Noise Case In the Presence of Process Noise Figure 4.5: The flight velocity using the adaptive LQR in the presence of vertical gust disturbances(Section4.3) 0 5 10 15 20 25 30 35 40 45 50 −800 −600 −400 −200 0 200 400 600 800 time (sec) control input, u For No Noise Case In the Presence of Process Noise Figure 4.6: The control input using the adaptive LQR in the presence of vertical gust disturbances(Section4.3) 70 4.4 PerformanceLossDuetoSensorInaccuracies If perfect measurements are available where the output views all the states with perfect sensors, the LQR design inherently provides the optimal control in the case of a Gaus- sianwhitenoisedisturbance. Underthemorepracticalassumptionofimperfectsensory information,however,theoptimalcontroladdressedbytheLQRisnomoresatisfactory. When the states and thus the output of the system are corrupted by disturbance sig- nals, those noise-corrupted measurements might severely affect the performance of the aircraftthroughthestatefeedbackandtheadaptivelaw. An additive zero-mean white Gaussian sensor noise with unity power spectral den- sity is applied separately at each state sensor for simulation purposes. The correspond- ingresultsusingtheLQRdesignarepresentedinFigure4.7,Figure4.8,Figure4.9and Figure 4.10 in comparison to the aircraft response where there is no gust disturbance and perfect measurements are available. As can be observed from these results, the oscillationsinthecontrolsignalandthealtitudechangeareunacceptable. 71 0 5 10 15 20 25 30 35 40 45 50 0 1000 2000 3000 4000 5000 6000 7000 time (sec) distance covered, δx (ft) No noise case Process Noise and Measurement Noise Figure4.7: TheflightdistancecoveredusingtheadaptiveLQRinthepresenceofverti- calgustdisturbancesandmeasurementnoise(Section4.4) 0 5 10 15 20 25 30 35 40 45 50 −600 −500 −400 −300 −200 −100 0 100 time (sec) altitude change, δh (ft) No noise case Process Noise and Measurement Noise Figure 4.8: The corresponding altitude loss using the adaptive LQR in the presence of verticalgustdisturbancesandmeasurementnoise(Section4.4) 72 0 5 10 15 20 25 30 35 40 45 50 −50 0 50 100 150 200 250 time (sec) horizontal velocity, V x (ft/sec) No noise case Process Noise and Measurement Noise Figure 4.9: The flight velocity using the adaptive LQR in the presence of vertical gust disturbancesandmeasurementnoise(Section4.4) 0 5 10 15 20 25 30 35 40 45 50 −1000 −800 −600 −400 −200 0 200 400 600 800 1000 time (sec) control input, u No noise case Process Noise and Measurement Noise Figure 4.10: The control input using the adaptive LQR in the presence of vertical gust disturbancesandmeasurementnoise(Section4.4) 73 4.5 PerformanceLossRecovery: AdaptiveLQGDesign Thestatesareavailableformeasurement,andatfirstsightthereisnoneedtodesignan observer to estimate the full state vector, x s (t) in order to make use of the advantages state-feedbackdesignhastooffer. Whenthesensornoiseactstoimpairthesemeasure- ments, however, the states are driven by noise and are thus also random processes. The designofanestimatorisrequiredforsuchstochasticsystems. We use the LQG problem setup based on the stochastic-deterministic dualism and the Separation Principle to design the feedback controller and the observer separately. NotethatiftheestimatesprovidedbyaKalman-Bucyfilterareusedintheadaptivelaw, theupdatemechanismwouldalsobecomemoreaccurate. If the measurement noise, γ is posed at the sensors modifying the output relation, theaircraftresponsecanbediscussedusingthestochasticdynamicalsystem: ˙ x s =Ax s +Bu+ξ (4.9) ξ =Gω g y =Cx s +γ where the gust disturbance, ω g ∈R and the measurement noise, γ∈R 5 are assumed tobezero-meanGaussianwhitenoisewithconstantpowerspectraldensitymatrices. Assumingthat ξ and γ areuncorrelated,wedefine E{ω g (t)ω g (τ)} = Ω g δ(t−τ) (4.10) E{γ(t)γ T (τ)} = Γδ(t−τ) (4.11) E{ξ(t)γ T (τ)} = 0 (4.12) E{γ(t)ξ T (τ)} = 0 (4.13) 74 where Ω g ∈R and Γ∈R 5×5 are spectral density matrices, E{.} is the expectation operator,and δ(t−τ) isadeltafunction. Sincethestochasticprocess,ω g (t) isincorporatedintotheaircraftdynamicsthrough the G matrix,fortheresultingstochasticprocessnoise, ξ weobtain E{ξ(t)ξ T (τ)} =E{Gω g (t)ω g (τ)G T } (4.14) =GG T E{ω g (t)ω g (τ)} =GG T Ω g δ(t−τ) where Ξ =GG T Ω g ∈R 5×5 isthespectraldensitymatrixof ξ. Thelevelsofprocessandmeasurementnoisemightbeenvironmentaldependent. If theprocessandmeasurementnoisecovariancematricesarenotavailable,theycouldbe estimated,andseveralsuchmethodsarediscussedin[Mehra,1970,GutmanandVelger, 1988,Yaz et al., 1999,Hide et al., 2004]. A fuzzy adaptive Kalman filter with weaker reliance on the prior statistical information is presented in [Xiong et al., 2005]. On the otherhand,asstatedin[StevensandLewis,1992],itisnotunreasonabletoassumethat the measurement noise has zero mean since there should be no bias on the measuring instruments. Wecancalculatetheoptimalfiltergainmatrix, K f throughtherelation K f =SC T Γ −1 (4.15) where S satisfiestheFilterAlgebraicRiccatiEquation(FARE): AS +SA T +Ξ−SC T Γ −1 CS = 0 (4.16) TheKalman-Bucyfiltercanthenbeimplementedintheformofastateobserveras: 75 ˙ ˆ x s =Aˆ x s +Bu+K f (y−Cˆ x s ) (4.17) where ˆ x s ∈R 5 isthevectorofstateestimates. Notethatalthoughtheconstantnegative throttlesettingtermisignoredinthefilterdesignphase,itisincludedinthesimulations andasymptoticallyrejectedbytheintegratorinthecontroller. The LQR design is based on the assumption that the system is stabilizable and detectable. To extend the LQR solution to the LQG controller, we further assume that (A T ,C T ) is stabilizable and (A T ,E) is detectable where Ξ=E T E and Γ>0. Theseassumptionsguaranteetheexistenceofthesolution S>0 thatresultsinastable (A−K f C)matrixasdiscussedin[Doratoetal.,1995]. If the gust noise, ω g is not white, the noise-shaping filter proposed in [Lewis and Syrmos,1995]canbeemployedas: ˙ x ω =A ω x ω +B ω n (4.18) ω g =C ω x ω +D ω n where n(t) is a white noise input, and ω g (t) is the filter output. The filter matrices, A ω ,B ω ,C ω ,D ω canbedeterminedbasedonfactoringthespectraldensityofthenoise, ω g (t) asdiscussedin[Lewis,1986]. The aircraft model and the noise-shaping filter can then be combined to obtain the augmenteddynamics ˙ x ωaug =A ωaug x ωaug +B ωaug u+G aug n (4.19) y =C ωaug x ωaug +γ 76 Figure4.11: AdaptiveLQGControlStructure(Section4.5) where x ωaug = x s x ω , A ωaug = A GC ω 0 A ω , B ωaug = B 0 , G aug = CD ω B ω , C ωaug = h C 0 i A nonwhite γ(t) signal can be addressed by a similar method. By solving the FARE for the augmented system, thesame filterdesign methodcan beapplied withthe requiredmodifications. TheadaptiveLQGcontrolstructureisdemonstratedinFigure4.11. Theparameters of the Kalman-Bucy filter along with the parameters of the LQ controller are updated usingthemodelparameterestimatesgeneratedbythesamerobustadaptivelawandthe correspondingsolutionsoftheFAREandtheARE. 77 4.6 AdaptiveLQRandAdaptiveLQGPerformance Comparison TheperformanceoftheadaptiveLQRandthatoftheadaptiveLQGcontrollerarecom- pared via simulations in the presence of vertical gust disturbances and sensor measure- ment noise where Ξ = GG T and identity Γ are the corresponding power spectral densities. Figure 4.12, Figure 4.13, Figure 4.14 and Figure 4.15 display the responses ofthetwosystems. IftheadaptiveKalman-Bucyfilterisincludedinthedynamics,thereisashortdelay observed when the UAV leaves the thermal location to glide towards the next thermal, which is definitely an acceptable compromise. The adaptive LQG controller results in reliable aircraft performance as the oscillations are efficiently suppressed. The system responsebecomesratherpredictableinthepresenceofrandomnoise. 78 0 5 10 15 20 25 30 35 40 45 50 0 1000 2000 3000 4000 5000 6000 7000 time (sec) distance covered, δx (ft) LQR LQG Figure4.12: Theflightdistancecomparisoninthepresenceofverticalgustdisturbances andmeasurementnoise(Section4.6) 0 5 10 15 20 25 30 35 40 45 50 −600 −500 −400 −300 −200 −100 0 100 time (sec) altitude change, δh (ft) LQR LQG Figure4.13: Thecorrespondingaltitudelosscomparisoninthepresenceofverticalgust disturbancesandmeasurementnoise(Section4.6) 79 0 5 10 15 20 25 30 35 40 45 50 −50 0 50 100 150 200 250 time (sec) horizontal velocity, V x (ft/sec) LQR LQG Figure4.14: Theflightvelocitycomparisoninthepresenceofverticalgustdisturbances andmeasurementnoise(Section4.6) 0 5 10 15 20 25 30 35 40 45 50 −1000 −800 −600 −400 −200 0 200 400 600 800 1000 time (sec) control input, u LQR LQG Figure 4.15: The control input comparison in the presence of vertical gust disturbances andmeasurementnoise(Section4.6) 80 4.7 Conclusion Inoptimalstaticsoaringapplicationsthecontrolobjectiveistotracktheoptimalveloc- ity along the flight path. The remote thermal data is in general assumed to be accu- rate, atmospheric disturbances acting on the system during the inter-thermal flight are ignored,andperfectsensoryinformationisrequiredforthecontrollerdesignandimple- mentation. A more realistic flight environment can be described, however, only if these assumptionsarerelaxed. Furthermore,aquantificationoftheeffectoferroneousthermal dataontheperformanceofthegliderUAVcanbesupportiveinreliabilityassessment. Inthischapterverticalgustdisturbancesandsensormeasurementnoisearemodeled as stochastic processes. An adaptive LQG design including an optimal state estimator inadditiontoanoptimalstatefeedbackcontrollerisimplemented. Thesensormeasure- ments are reconstructed using the noise rejection capabilities of an adaptive Kalman- Bucy filter. These enhanced signals are incorporated with the feedback control and the adaptive law. The simulation results demonstrate that the achieved flight performance is comparable to the case with no process or measurement noise, and thus the adaptive control scheme is improved by exploiting the practical statistical properties of random noise. 81 5. OPTIMALSTATICSOARINGOFUAVS USINGVEHICLEROUTING WITHTIMEWINDOWS 5.1 Introduction The performance in soaring flight depends on the experience, aggressiveness and the intuitive decisions of the pilot. The level of autonomy to decide the flight path for optimal soaring can be increased if a more systematic approach possibly involving a shortestpathplanningmethodologyisfollowed. Shortest paths are used in geographical information systems, military mission plan- ning, autonomous robot route planning, terrain following ‘map-of-the-earth’ aircraft route planning, and city regional planning systems. [Saab and VanPutte, 1999] Other areas of shortest path applications include shipping/distribution problem, economics, critical path analysis, circuit board layout and Very Large Scale Integration (VLSI) design, remote sensing, and automated traveling advisory systems. If digital represen- tations of thermals above a region of interest are available in the form of a database, a route planning system can also be developed for gliders in order to fully benefit from thesecurrents. Inordertomeetthedistinctivesetofchallengesinroutegenerationforgroundvehi- cles operated in metropolitan areas a VRP formulation is most often used. The VRP 82 calls for the determination of the optimal set of routes to be performed by a fleet of vehicles to serve a given set of customers, and it is one of the most important and stud- iedcombinatorialoptimizationproblems.[TothandVigo,2002] The VRP was first introduced in [Dantzig and Ramser, 1959]. The models and algorithms proposed for the solution of vehicle routing and scheduling problems can be used effectively not only for the solution of problems concerning the delivery and collection of goods but for the solution of different real-world applications arising in transportationsystemsaswell.[TothandVigo,2002] We propose that optimal static soaring can be considered as another application area exploiting methods already developed for collection and delivery of goods. We show that optimal routing for thermal soaring problem can be studied in more detail by employing a VRP formulation. Our objective is to develop a systematic method to assign thermals in an efficient manner via defining the similarities in the nature of the soaring problem to the VRP. Furthermore, if additional complexities such as thermals with limited lifespan, stochastic thermal locations, or several gliders soaring simulta- neously are to be considered, optimal static soaring problem can still be conveniently stated as one of several extensions of the VRP, and this interpretation allows the use of correspondingmethodsdevelopedfortheVRPsolution. We consider in particular the limited lifespan for thermals which is in analogy with the time window constraints imposed at the customer locations for the VRPTW. [Kahveci et al., 2007c] We present the thermals above a terrain using a graph theoretic approachanddefineanoptimizationproblemtochoosetheflightrouteandthevelocity. Based on the time cost function accordingly determined we develop an exact algorithm forthermalassignmentandobtaintheoptimaltime-feasibleflightpaththatsatisfiesthe timeconstraints. 83 The details of the restatement of the static soaring problem as a VRP are given in Section 5.2. Section 5.3 constructs the static soaring problem on a 2D map. Section 5.4 presents a mathematical model of the problem as a VRPTW. In Section 5.5 the flight arc cost function and the optimal flight velocity for each inter-thermal glide are evaluated. OurexactalgorithmisdescribedindetailinSection5.6. Section5.7presents thesimulationresultsforaspecificproblemscenario. Othersoaringflightissueswhich can be interpreted as various VRP extensions are briefly described in Section 5.8. Our conclusionsformSection5.9. 5.2 AnAnalogy: StaticSoaringProblemasaVRP TheUAVstartsatadesignatedpointtoreachthedestinationinminimumpossibletime. For soaring purposes utilizing thermals only on the straight flight path might be inade- quateasitmightnotbeeitherrealisticornecessarilyoptimal. Itismostoftenpreferable tobenefitfromageographicallydispersedsetofthermalsbydeviatingfromtheoriginal straightflightpath. To keep the problem tractable, we pose an additional requirement at the thermal locations. OncethegliderUAVstartssoaringatacertainthermallocation,neitherdoes the glider quit until the altitude loss due to preceding inter-thermal glide is recovered, nor does it stay inside and climb the thermal any further after gaining back the lost altitude. If we investigate the relation of the optimal static soaring problem with the generic VRP,ananalogycanbedrawn. ThefeaturesdefinethegliderUAVasthesinglevehicle in service, and the thermal locations as customers, terminals or intermodal facilities. One can view the climb time at the thermal locations as the service time for loading or unloading, the set of allowed inter-thermal flight arcs as the road network, and the 84 optimal soaring path as the optimal vehicle route. For a fleet of trucks delivering cus- tomer products to outlets, reduced service delays can be interpreted as better service, whereas for a glider UAV in soaring flight reaching a destination in minimum possible time can be interpreted as the best performance. Consequently, the total time spent for gliding and climbing is chosen as the corresponding travel cost. The fact that thermals mightbeavailableforsoaringonlyforalimitedperiodoftimecanbeinterpretedasthe time window requirement of a customer during which the customer can be served. The crucial maximum working time limit for drivers can also be directly interpreted for the UAVsincetheendofdaylightmeansalsotheendofthermalsoaring. One can formulate a flight in the out-and-return category where the UAV starts at a particularpointonthemapandendsatthesamepointafterithasvisitedapriorideclared turn points in an order decided by the pilot during flight. This would correspond to finding a tour connecting some or all of the customer locations and returning back to the depot in the VRP statement while meeting a budget constraint. In this chapter for simplicitywedonotrequiretheglidertoreturnbacktowhereitstarts,butonlytoarrive ataspecificdestinationpoint. 5.3 StaticSoaringona2DMap The flight-time optimization for a glider UAV in an autonomous static soaring applica- tion can be conveniently stated as a graph problem on a two-dimensional area map of thermalsdetectedoveraspecificarea. Thismaphastwodegreesoffreedom, x and y, lookingdownperpendicularonthemaneuveringareaoftheUAV.Eachthermallocation represents a node which is connected to several adjacent ones under certain distance criteria. TheunderlyingthermalnetworkcanthenberepresentedwithagraphG(N,A) 85 where N ={0,1,...,n} is the set of nodes including the start node and the thermal locations,and A isthesetofedgeswhicharethestraightflightarcsbetweenthenodes. We define two metrics on each arc (i,j) of the graph: A spatial metric, l ij and a temporal metric, c ij . The spatial metric is defined as the Euclidean norm of the flight distanceinbetweeneachpairofthermallocations i and j : l ij =l ji = q (x i −x j ) 2 +(y i −y j ) 2 , {i,j}∈N , (i,j)∈A (5.1) whereas the arc cost function, c ij is defined as the corresponding time that would be spentincasetheUAVincludesthearc(i,j)∈A initsflightroute. We employ the spatial metric to decide which thermal locations would be directly connected on the map. A maximum value, l max is defined, and for any arc (i,j)∈A where l ij >l max , the corresponding temporal metric is assigned to infinity, c ij =∞. Consequently, l max posesarestrictiononthedirectflightrangeas c ij =∞ impliesthat flightfromnode i tonode j isnotallowedunlessotherthermallocationsarevisitedin between. Hereweassumethatthedatabaseofarcssatisfyingthespecifiedrangerestrictionis readily available, or it can be formed by checking each node pair on the graph. If the thermalsaretoolargeinnumber,analgorithmcanbedevelopedtosimplifythedistance check procedure in between the nodes. One might refer to grid method summarized in[Sedgewick,1988]asasimplebuteffectivetechniqueformaintainingproximityrela- tionshipsamongpointsintheplanetoconstructanartificialgridthatdividestheareato besearchedintosmallsquaresandtokeepshortlistsofpointsfallingintoeachsquare. As we model a nonexisting or infeasible flight arc (i,j) ∈ A using an infinitely largecost c ij ,wecanassumethatthethermalnetworkisfullyconnected. Afinitecost c ij of the arc (i,j)∈A is defined to include the optimal pure glide time along the arc 86 (i,j)andthecorrespondingthermalsoaringtimeatnode j.Thestrengthofthethermal atnode j andtheflightdistancecorrespondingtothearc(i,j)address c ij ,whichisa measureoftimeandisthuspositivebydefinition. The thermal map is basically a directed graph (digraph) since c ij for (i,j) ∈ A and c ji for (j,i)∈A are not necessarily equal which in turn implies that the soaring problemisasymmetric. Alsonotethat c ij =l ij = 0 istobedefinedfor i =j. Forthespatialmetriconarcs,atriangleinequalityfollows: l(i,j)<l(i,k)+l(k,j), ∀ i, j, k ∈N (5.2) Suchinequalityisnotrequiredtoholdforthecostfunctionwhichentailsthat c ij ≤c ik +c kj , ∀ i, j, k ∈N (5.3) does not have to be satisfied. Otherwise the shortest path between the nodes i and j wouldalwaysbethesinglearc(i,j)∈A,andtheproblemwouldbecomenonexistent. Thenodesalsohavetimewindowsassociatedwiththemwhichrepresentthethermal lifespanperiodsposedforthermallocations. Atimewindowisdefinedasapair[a i ,b i ] where a i istheearliesttimethethermalatlocation i isavailableforsoaring,and b i is thetimethethermalvanishes. Accordinglythewidthofatimewindow, b i −a i denotes thelifespanofthatthermal. IntheVRPapplicationsitissuggestedtoallowvehiclestowaitatanodeifthevehi- cle arrives earlier than the customer’s availability for service. [Malandraki and Daskin, 1992] In such case with a hard time window constraint if the vehicle is too late rather than being early, a penalty for lateness is incurred in the objective function. However, we do not allow the glider wait at a thermal location as it would meanwhile keep los- ing altitude. Although it is not possible to climb a thermal if the aircraft arrives after 87 the thermal vanishes, we define no penalty if the aircraft does not arrive as soon as the thermalisavailable. The sequence of thermal locations on the graph through which the UAV navigates fromstarttodestinationiscalledarouteonthegraphwhereeachsuccessivepaironthe route defines an arc in A. The flight cost for any route in G is the sum of the costs of thearcsdefiningtheroute. Foraroutetobefeasible,timewindowsforthecorrespondingsetofthermalsshould berespected,andthetotalcostaccumulatedalongthatrouteshouldbefinite. 5.4 Mathematical Formulation of Optimal Static SoaringasaVRPTW Once the problem is represented by a 2D graph as motivated in the previous section, a basic model for the VRPTW can be modified to the optimal static soaring problem describedasfollows: min X i∈N X j∈N c ij x ij (5.4) subjectto X i∈M−{n} x ij = 1 ∀j∈M−{0}, (5.5) X j∈M−{0} x ij = 1 ∀i∈M−{n}, (5.6) 88 X j∈N x 0j = 1, (5.7) X i∈N x in = 1, (5.8) ω i +t thermal si +t glide ij =ω j ∀ i∈M−{0,n}, j∈M−{0,Δ + (0)}, s = Δ − (i), j = Δ + (i), (5.9) a j x ij ≤ω j ≤b j x ij ∀ j∈M−{0}, j = Δ + (i), (5.10) a k ≤ min i∈M−{0} b i , k = Δ + (0), (5.11) b n ≥ max i∈M−{0} a i , (5.12) x ij = 1, ∀ i, j∈M, j = Δ + (i) 0, else (5.13) Intheabovemathematicalmodel[a i ,b i ]isthetimewindowforthethermallocation i. The values of a i and b i represent the earliest and the latest times that a thermal at location i isavailableforsoaring. x ij istheflowvariablefortheflightarc(i,j)∈A whichtakesvalue 1 ifthearc(i,j)∈A belongstotheoptimalrouteand 0 otherwise. Δ − (i) is the node that precedes and Δ + (i) is the node that follows node i along the optimalroute. t thermal ij isthetimerequiredforsoaringatnode j inordertogainback thealtitudelostduringtheglidealongthearc(i,j)where t glide ij istheglidetimealong thatarc. Thevariable ω i ,i∈M−{0} specifiesthetimeaircraftstartsthermalsoaring atthethermallocation i. M−{0} isthesetofclimbedthermalswhere i = 0∈M denotesthestartnode,and i =n∈M isthedestination. The constraints (5.5) and (5.6) guarantee that a UAV which enters a thermal for soaring also leaves that thermal after its climb unless it is a thermal located at the des- tination. They further imply that no circles are allowed in the optimal flight route and thattherouteisconnected. Theconstraints(5.7)and(5.8)areimpliedbytheconstraints (5.5)and(5.6)andthuscanberemoved. Theconstraint(5.9)enforcessoaringstarttime 89 conditiononconsecutivethermals. Theconstraint(5.10)ensuresthateachtimewindow isrespectedalongtheoptimalroute. In(5.11)itisimpliedthatthefirstthermallocation visited should have its time window’s start time earlier than the end time of any of the time windows thatare posedatthethermallocationsvisitedintheoptimalroute. Simi- larly,in(5.12)itisguaranteedthatthetimewindowposedonthethermallocatedatthe destination does not end before the thermals on the optimal route become available for soaring. Thebinaryconstraint,alsocalledtheintegralityconstraint,isgivenin(5.13). NotethatMcanbeanysubsetofN aslongasitincludesthestartandthedestination nodes. The glider UAV, however, has to visit all the nodes in that subset. The task is then to determine the set M such that the UAV visits each node in the set M exactly onceinanordertominimizethetotalflighttimeanddoesnotviolatethetimewindows. As a result, it is required that the vehicle travels between two distinguished nodes, 0 and n,whichisin[GuptaandKrishnamurti,1997]referredtoas S−T Routing. A model based on time-dependent travel speeds is developed in [Hill and Benton, 1992] where the average speeds in certain areas during different time periods are con- sidered. Inthestaticsoaringproblemthatwediscuss,theglideraimstoflyataspecific constant horizontal velocity from one thermal location to the other. The optimal flight velocityforthatglideiscalculatedinthenextsection. 5.5 FlightArcCostandOptimalFlightVelocity The faster the aircraft glides from the thermal location i towards j, the more it sinks duringitsglidealongthearc(i,j)andthelongeritneedstoclimbatthethermallocation j ifitistorecoverthataltitudeloss. Thisimpliesthatthemaximumallowedvelocityis ingeneralnottheoptimalglidespeedtobechosen. 90 The glide time is a function of the glide velocity, V x ij and the distance, l(i,j). The time spent for thermal soaring, on the other hand, depends on the altitude loss, h loss (V x ij ,V z ij ) during the glide and the strength of the thermal ahead, W j . The total costforthearc(i,j)isthendefinedas: c ij =t total ij =t glide ij +t thermal ij = l(i,j) V x ij + h loss (V x ij ,V z ij ) W j (5.14) h loss (V x ij ,V z ij )=l(i,j)× |V z ij | V x ij =l(i,j)× −V z ij V x ij (5.15) V z ij ∼ =aV 2 x ij +bV x ij +c (5.16) In the above equations V x ij is to be chosen for each arc (i,j) such that the corre- sponding c ij isminimized. TheflightvelocityoptimizationproblemintroducedinChapter2canthenbemodi- fiedandsummarizedaccordingly: V x ij = arg min Vx ij ∈V 0 t total ij (5.17) subjectto t glide ij = Z l(i,j) 0 1 V x ij dl (5.18) ω j =ω i +t thermal si +t glide ij (5.19) h(ω j +t thermal ij ) =h(ω i +t thermal si ) (5.20) ⇒ t thermal ij = h(ω i +t thermal si )−h(ω j ) W j (5.21) ˙ x s =Ax s +Bu+Gω +Hδ T (5.22) 91 ω(t) = 0, ω i +t thermal si <t≤ω j W j , ω j <t≤ω j +t thermal ij (5.23) s = Δ − (i), j = Δ + (i) (5.24) V 0 def = {V x ij ∈R | V min ≤V x ij ≤V max } (5.25) As the quadratic approximation to the polar curve is assumed to be a truthful rep- resentation of the performance of the aircraft, the optimal horizontal velocity for the inter-thermalflightcanbewritteninclosedformas: V opt ij = r c−W j a , c−W j a > 0 (5.26) Notethattheoptimalhorizontalvelocitytoapproachathermallocationj fromnode i, V opt ij depends on the thermal strength, W j and the quadratic approximation to the polarcurveoftheglider,butnotontheflightdistance, l(i,j). Solving this particular VRP is hence equivalent to deciphering the shortest time- feasible path from the start point to the destination where V opt ij is the horizontal flight velocity to be tracked along each inter-thermal flight arc (i,j) in the optimal route. In thenextsectionwedescribeanexactmethodologythatprovidesthesequenceofthermal locationstobevisitedinordertoformtheoptimalroute. 5.6 AnAlgorithmicFrameworkBasedonShortestPath Search Three standard algorithms for the shortest path problems can be listed as the Bellman- Ford Algorithm, Dijkstra’s Algorithm, and the Floyd-Warshall Algorithm. The first two find the shortest path between one specific node and the others, whereas the third 92 onefindstheshortestpathsbetweenallthenodesinthenetwork. However,asdiscussed in[BertsekasandGallager,1987],itturnsoutthatifonesimplywantstofindtheshortest path from one node i to another node j, there is no known way of doing this without essentially finding the shortest paths from all nodes to j or finding the shortest paths fromnode i toallothernodes. WeemploytheFloyd-WarshallAlgorithmwhichcarries out n steps of iterations where n is the number of nodes and provides a predecessor matrix, Π showingthenodestobevisited. The Floyd-Warshall Algorithm takes single arc costs as starting estimates. It then calculates shortest paths under the constraint that only node 1 can be used as an inter- mediatenode,andthenwiththeconstraintthatonlynodes 1 and 2 canbeused,andso forth.[BertsekasandGallager,1987] Thealgorithmcanbesummarizedasfollows: c (k) ij = c ij , k = 0 min(c k−1 ij ,c k−1 ik +c k−1 kj ), k≥ 1 (5.27) π (0) ij = nil, i =jorc ij =∞ i, i6=jandc ij <∞ (5.28) π (k) ij = π (k−1) ij , c (k−1) ij ≤c (k−1) ik +c (k−1) kj π (k−1) kj , c (k−1) ij >c (k−1) ik +c (k−1) kj ThelaststepofiterationsdefinesthecorrespondingentriesofmatricesC f =C (n) ∈ R nxn and Π f = Π (n) ∈R nxn whichdeterminetheshortestpathsforall i, j∈N . We use the above iterative implementation of the traditional Floyd-Warshall Algo- rithm. Yet, in [Park et al., 2004] cache-friendly optimizations to the Floyd-Warshall Algorithm are developed where a recursive implementation is shown to reduce the processor-memory traffic, and a tiled implementation is proven to be asymptotically 93 optimal among all implementations of the Floyd-Warshall Algorithm. We do not dis- cuss the details of these improvements. However, the reader should be aware of the recentadvancesinrealexecutiontimeforthegraphalgorithms. Aswehavetimewindowconstraintsamethodologyprovidingatime-feasiblesolu- tion is needed. For the problem of ranking shortest paths, one might first aim to rank k number of paths in order to check them for time-feasibility later on. Algorithms for the k shortestpathsproblemarediscussedin[Eppstein,1998]wherealistofprevious work on the subject is also given and the previously known bounds are shown to be improved by an order of magnitude. By checking the k shortest paths list for feasi- bility, however, one might not be able to reach the optimal solution as the value of k cannot be directly specified in advance. This method thus causes a serious drawback sinceitactuallyrequiresrankingallpossiblepathsinourcasewherethecorresponding computationaleffortrisesrapidlywithincreasing k. Among the shortest path planning techniques are the two extremes: To precompute andstoretheentirefinitesetofallpossiblequerieswithaverylargestorageoverhead,or toeliminatetheprecomputationstep. Theoptimalsolutionisreportedtobesomewhere inbetweenthetwoextremeswhichisusingsomeinitialprecomputationefforttoobtain along-termgaininsearchtime.[PreparataandShamos,1985] To solve constrained shortest path problems, state of the art solvers compute lower and upper bounds on the problem and then close the duality gap. [Sellmann, 2003] Amidthesetighteningstrategies,cost-basedpreprocessingtechniquesaredescribedand a repeated problem reduction procedure is proposed for hard constrained problems in [DumitrescuandBoland,2003]. Building on these ideas we start with a preprocessing filter step and keep reducing thesizeofthenetworkbyeliminatingflightarcswhichareshowntobeinfeasible. Our 94 methodologycanbesummarizedinthreesteps: Preprocessing,RouteOptimization,and RouteValidation. Preprocessingaimstoimprovethesearchtimebyidentifyingasetofinfeasibleflight arcs which are shown to cause inevitable violations of the time constraints if they are includedinaflightroutenomatteratwhatstageoftheroutetheyareemployed. This action space reduction might discover much about the structure of the graph and reduce the computational time for the subsequent steps. We base this filtering step onthreesetsofinitialconditionsrequiringthetimecost, c ij tobesettoinfinity: a) The thermal at location j vanishes early, and the aircraft that glides along arc (i,j)cannotcompleteitsrequiredclimbatnode j beforethethermal’sendtime: b j <a i +t glide ij +t thermal ij =a i +c ij (5.29) b)Thethermalatlocation j isnotyetavailableforsoaringwhentheaircraftarrives: a j >b i +t glide ij (5.30) c) The lifespan of the thermal at location j is not long enough, and the aircraft cannotfullybenefitfromsoaringatthatlocation: b j −a j <t thermal ij (5.31) Thearcswithfinitecostshavethecostfunctiongivenbytheanalyticalexpressionin (5.14)wheretheoptimalvalueof V x ij isdefinedin(5.26). Thesetofarcswhichmight be involved in the feasible solution are thus defined. This Preprocessing step can be demonstratedtobesurprisinglyeffectiveinreducingtheproblemsizefordensethermal mapswithoutdegradingthesolutionquality. 95 In Route Optimization step we relax all temporal constraints. The arc costs from the Preprocessing step can be used to find the shortest path by employing the Floyd- WarshallAlgorithm. Anoptimalrouteisthusproposed. To validate the proposed optimal route, in Route Validation step we check the fea- sibility of each arc (i,j) in that route with respect to the associated time constraints. If theassignedsoaringtimeatthethermallocation j doesnotfitthetimeintervalduring which that thermal is available for soaring, the arc cost, c ij is set to infinity, and thus the arc (i,j) ∈ A is eliminated from further consideration. If there are two or more violations determined, we only reconstruct the cost of the arc corresponding to the first violationintime,andtheothercostsinthatroutearekeptattheirfinitevalues. Although a time feasibility check is provided in Preprocessing step, Route Opti- mization might still give an infeasible output. This is the reason why a validation step is still required. The Preprocessing step is not route-specific in the sense that the time to start and the time to end soaring at a thermal location actually depend on the path taken before arriving at that node. The conditions given in (5.29,5.30,5.31) thus might beinsufficientfordetectingseveraldifferentinfeasibilities. Theconditionin(5.29)does notillustratetheinfeasibilitywhentheaircraftcompletesitssoaringatnode i attime: t =ω i +t thermal si >a i , s = Δ − (i) (5.32) andthusfinishesitsclimbatthesuccessivenode j attime: t =ω i +t thermal si +t glide ij +t thermal ij >a i +t glide ij +t thermal ij , s = Δ − (i) (5.33) The condition in (5.30) does not provide any guarantees that the thermal located at node j would be available for soaring otherwise. It is highly probable that the aircraft leavesthethermallocation i atsometime: 96 t =ω i +t thermal si <b i , s = Δ − (i) (5.34) andthusstartssoaringatthethermallocation j atarelativelyearliertime: t =ω i +t thermal si +t glide ij <b i +t glide ij , s = Δ − (i) (5.35) On the other hand, the condition in (5.31) cannot detect the mismatch between the thermal lifespan and the time interval required for thermal soaring if the glider starts soaringatnode j withthesame t thermal ij attime: t =ω j >a j (5.36) If the time feasibility is shown in Route Validation step, the optimal solution is defined as that particular shortest path with no further calculations. If there are any temporal constraint violations detected, one returns back to Route Optimization step after the cost reassignment. The search is likely to intensify if the temporal constraints appear to be more restrictive. Additionally, at each recursive step the shortest path calculateddefinesalowerboundonthecostoftheoptimalsolution. In Preprocessing and Route Validation steps, we eliminate the arc (i,j) if the time window at node j is not respected. This is essentially based on the assumption that theoptimalpathshouldnotincludeanycyclesasotherwisethecycleeliminationwould result in a more optimal solution. If cycles are to be allowed, the VRP formulation shouldberestatedtolettheflowvariable, x ij tobegreaterthan1. For the static soaring problem we consider scenarios with relatively small number of thermals where an exact method is capable of finding the optimal path. However, as the number of thermals under consideration gets larger, the computational time would increaseasexpected. Itisstatedin[TothandVigo,2002]thatthelargestVRPinstances that can be consistently solved by the most effective exact algorithms proposed so far 97 containabout50customers,whereaslargerinstancesmaybesolvedtooptimalityonlyin particular cases. Many of theproficientVRP solutiontechniquesthat efficientlyhandle arbitrary numbers of vehicles and waypoints are heuristic which can typically provide highqualitysolutionswithmoderatecomputationalcosts.[Stacketal.,2004] Heuristic algorithms provide simplicity, potentially at the cost of accuracy or preci- sion where in general a feasible solution is built first and then improved performing a sequenceofedgeorvertexexchanges. In[LinandKernighan,1973]onebasicapproach toheuristicsforcombinatorialoptimizationproblemsissaidtobetheiterativeimprove- ment of a set of randomly selected feasible solutions. Metaheuristic algorithms such as SimulatedAnnealing,AntColonyOptimization,TabuSearchandGAskeeptheempha- sis on performing a deep exploration of the most promising regions of solution space. In [Desrosiers et al., 1995] a review of the heuristic methods for the VRPTW is given, and an extensive computational study of their performance is conducted in [Solomon, 1987]whereitisreportedthatseveralheuristicsperformwellindifferentproblemenvi- ronments with different percentage of customers having time windows and different tightnessorpositioningofthesewindows. Likewise, one might return to approximate methods such as heuristics to solve the staticsoaringproblemifanexactmethodisnotrequired. Inthatcase,onceinfeasibility ofarouteunderconsiderationisshown,thenextroutetobeconsideredmightbeallowed to differ from the previous one by the inclusion of a random arc on the graph, and the subsequentfeasibilitycheckmightbeappliedsimilarly. 5.7 ASimulationScenario The glider UAV starts at the location given by the coordinates (x,y) = (0,0) which correspondstonode0onthegraph. Atotalofseventhermalsaredetectedtobeavailable 98 for soaring where the strength and the lifespan of each thermal are known. The goal of the aircraft is to reach the destination (x,y)=(16,000,−7,000) represented by node 7 inminimumtime. ThesimulationscenarioissummarizedinFigure5.1. The maximum flight arc distance, l max is given as 10,000ft which limits the number of finite-cost arcs in A to 36. The maximum allowed horizontal velocity is defined as V max =250ft/sec. If the optimal horizontal flight velocity is evaluated to begreaterthan V max ,thegliderUAVisassignedtoflyat V x ij =V max . Each flight arc is associated with a cost using the analytical expression given in (5.14). In Table 5.1, Table 5.2 and Table 5.3 numerical values of the time costs eval- uated for all arcs on the graph are listed. These time costs, t calc are compared with the measurements from simulations, t meas . The two are verified to match within an acceptable error range. The error is vastly due to the polar curve approximation used in (5.14) to evaluate t calc and also due to the transients of the tracking response in the simulationswhichinturnaffectthecorresponding t meas values. Figure 5.1: The UAV assigned to autonomously navigate from start to destination and thesetofthermalsdetectedtobeavailableforsoaringintheflightregion(Section5.7) 99 Table5.1: FlightTimeTable(I):SimulationScenarioinSection5.7 W=10ft/sec and V opt =204ft/sec Distance t glide t thermal t calc t meas (ft) (sec) (sec) (sec) (sec) 3,400 16.7 28.8 45.5 42.7 4,000 19.6 33.9 53.5 49.9 4,400 21.6 37.2 58.8 54.6 5,000 24.5 42.3 66.8 61.7 5,200 25.5 44.0 69.5 64.1 6,800 33.4 57.6 91.0 83.1 7,100 34.8 60.1 94.9 86.7 7,400 36.3 62.6 98.9 90.3 Table5.2: FlightTimeTable(II):SimulationScenarioinSection5.7 W=15ft/sec and V opt =233ft/sec Distance t glide t thermal t calc t meas (ft) (sec) (sec) (sec) (sec) 2,500 10.7 15.4 26.1 26.7 4,000 17.2 24.6 41.8 40.7 5,000 21.5 30.7 52.2 50.0 5,200 22.3 31.9 54.2 51.9 5,300 22.7 32.6 55.3 52.8 5,400 23.2 33.2 56.4 53.8 6,200 26.6 38.1 64.7 61.2 6,800 29.2 41.8 71.0 66.9 7,400 31.8 45.5 77.3 72.5 7,700 33.0 47.3 80.3 75.3 8,000 34.3 49.1 83.4 78.1 8,500 34.0 41.2 75.2 71.4 8,800 37.8 54.1 91.9 85.6 9,200 39.5 56.5 96.0 89.3 9,500 40.8 58.4 99.2 92.1 9,700 41.6 59.6 101.2 94.0 100 Table5.3: FlightTimeTable(III):SimulationScenarioinSection5.7 W=20ft/sec and V opt =259ft/sec ⇒ 250ft/sec Distance t glide t thermal t calc t meas (ft) (sec) (sec) (sec) (sec) 2,500 10.0 12.1 22.1 23.7 4,400 17.6 21.3 38.9 38.8 7,100 28.4 34.4 62.8 60.3 7,700 30.8 37.3 68.1 65.1 8,500 41.7 71.9 113.6 103.4 Figure5.2illustratesthe2Dgraphconstructedforthisscenarioincludingthecorre- spondingtimewindowforeachthermal. Theflightarcsfromthestartnodeandtheonesintothedestinationnodeareforward arcs inthe sense thatthey arealwaysawayfromthestart nodeand towardsthedestina- tion. Thetimecostsforarcs(1,0), (2,0), (7,3), (7,4),(7,5)and(7,6)arethusallset tobeinfiniteaswellasthearcsomittedfromdisplayonthegraph. There are a total of eight nodes on the graph including the seven thermal locations andthestartnode 0 forwhich[a 0 ,b 0 ]=[0,0]sec canbedefined. IftheFloyd-Warshall Algorithm is iterated n=8 times, final asymmetric cost matrix, C f =C (8) and the predecessor matrix, Π f =Π (8) defining the nodes on the shortest paths are obtained. Accordingly,theoptimalroutewouldbe 0−2−4−5−7 withatotalcostof 215.3 sec if the temporal constraints are ignored. Preprocessing step, on the other hand, assigns infinitely large costs to arcs (4,5), (5,1) and (5,4), and the route 0−2−4−5−7 is notfeasible. RouteOptimizationstepproposesthepath 0−1−5−7 withatotalcostof 216.0 sec as the optimal route. However, Route Validation step points at the temporal constraint violation at node 5 which is due to the aircraft’s arrival after the thermal vanishes. Notethattheaircraftisassignedtosoarthethermalatnode 5 from 94.8 sec to 151.3 sec, whereas the thermal lifespan is defined by the associated time window, 101 Figure5.2: StaticSoaringProblemona2Dgraph(Section5.7) 102 Figure 5.3: The thermals assigned by the solution method in order to achieve optimal soaringflightperformance(Section5.7) [a 5 , b 5 ]=[160, 400]sec. The algorithm thus sets c 15 to infinity and goes back to Route Optimization step. The Floyd-Warshall Algorithm is reapplied to obtain a new optimal cost of 216.4 sec which corresponds to the route 0−2−4−7. This route is shownto befeasible withrespectto thetimeconstraintsthroughRouteValidationstep, andFigure5.3isthe3Drepresentationoftheoptimalsolution. 5.8 OtherVRPInterpretationsforStaticSoaring The basic VRP formulation is an integral part of all its extensions. Similar to the VRPTW other variants such as the Capacitated VRP (CVRP), the Multi-Vehicle Rout- ing Problem, the VRP with Pickup and Delivery, the DVRP and the Stochastic Vehicle Routing Problem (SVRP) have their specific constraints which should be added to the problemformulation. 103 Thermals do not have demands relevant to cargo as trucks do. However, if there are anylimitsonthenumberofthermalsthatcanbesoared,theymightaswellbeinterpreted as defining the capacity of the vehicle. Solutions for the CVRP can hence be modified tosolvethestaticsoaringproblem. Ifavehiclecooperationisavailable,Multi-VehicleRoutingProblemisdiscussed. If there is a fleet of gliders in flight which are not allowed to soar inside the same thermal simultaneously due to safety restrictions, air traffic issues should also be considered. We develop an environment that corresponds to a single VRP, but it can potentially be generalizedtothemulti-vehiclecase.[KahveciandIoannou,2008] The VRP with Pickup andDeliveryposesintrinsicprecedenceconstraints. Itisdis- cussedin[SavelsberghandSol,1995]thatmanypracticalpickupanddeliverysituations are demand responsive, i.e., new transportation requests become available in real-time and are immediately eligible for consideration; as a consequence, the set of routes has tobereoptimizedatsomepointtoincludethenewtransportationrequests. Thoseissues allhavetheirequivalenceinthestaticsoaringproblemformulation. In real-world applications, traffic conditions change dynamically, possibly causing theinitialbestassignmenttobeinvalid. Frameworksthatarecapableofreactiveroutings arediscussedinliteratureinthecontextofaDVRP.Ifthestrengthofathermalistime- varying, the optimal route might need to be adjusted accordingly. The problem can then be discussed as a DVRP under varying levels of dynamism using relevant solution methods. We suppose that the travel cost for a glider’s flight along a specific arc is deter- ministic which is a basic assumption for most VRPs where network parameters are known a priori. If remote thermal detection data is not reliable, and the exact thermal strength is known only upon arrival at the thermal location, the static soaring problem could as well be formulated as an SVRP. The closed-form expressions and algorithms 104 are found to compute the expected length of an a priori sequence under general proba- bilistic assumptions in [Bertsimas, 1992], and they might also be useful for the thermal soaringproblem. 5.9 Conclusion In today’s soaring flight, intuitive decisions of pilots play the major role in flight path determination,whereasamoresystematicwayisdesiredtogeneratetheoptimaltrajec- tory. WepointatthefactthattheoptimalstaticsoaringproblemfitstheVRPstatement. Due to other practical considerations of soaring applications, flight path decision is in generalmorecomplicatedthansolvingadeterministicVRPinitsmostbasicform. The optimal static soaring problem with extra constraints, however, can still benefit from availablesolutionmethodologiesifextensionsoftheVRPareconsidered. Starting from a designated point we define the goal of the glider UAV as reaching some specific destination in minimum time via thermal soaring. Since the limited ther- mal lifespan is posed as an extra restriction at each thermal location, the aircraft also aims to climb the thermals while they are available for static soaring flight. The formu- lationstructureishencechosenasaVRPTW,anextensionoftheVRPwheretheservice mustbeprovidedwithinanassociatedtimewindow. The problem is constructed on a 2D graph, and a solution methodology based on shortest path algorithms is developed. If a route does not respect at least one of the time constraints, not all arcs in that route can be used in the optimal solution, and the problemspacecanbereducedbyproperarcelimination. Followingtheseideastheopti- mal time-feasible path is found by a recursive exact algorithm, the three steps of which arePreprocessing,RouteOptimizationandRouteValidation. Simulationresultsdemon- stratethedeterminedoptimalroutetobeconsistentwiththethermallifespanconstraints. 105 Besides, if very large scenarios are involved and the method becomes computationally cumbersome,ourtechniquecanbemodifiedusingheuristicalgorithms. Asaresult,byexploitingtheanalogybetweenthestaticsoaringapplicationsandthe VRP,thegliderUAVevaluatesitsoptimalflightpathautonomously. Theapproachalso appearstobevaluabletoassistpilotsinmannedsoaringflightswheretheycanmakeuse ofanonboardcomputationalsystemforrouteselectioninsteadofpurehumandecision. 106 6. AHEURISTICSEARCHALGORITHM FORMANEUVERINGOFUAVS ACROSSDENSETHERMALAREAS 6.1 Introduction InChapter5weformulatetheoptimumsoaringperformancerequirementsofaUAVas aVRP,arelevantpopularproblemintransportationresearch. Inthischapterwefurther proposeasystemizedoptimalflightpathdecisionprocedureforagliderUAVtocovera densethermalregion.[Kahvecietal.,2007b] The VRP requires the efficient use of a fleet of vehicles that make a number of stops to pick up and/or deliver passengers or products. Thermal soaring problem in analogy requires the glider to visit a number of thermal locations in the flight region. The similarities can be listed further if extended versions of the VRP are considered. The reader might refer to Chapter 5 for a variety of theVRP interpretations for thermal soaring problem and a discussion on how to exploit the corresponding VRP solution methodologiesforthesoaringpathdecision. In this chapter we consider a glider flying over a region dense with thermals and applyaheuristicsearchmethodinordertosolvethethermalsoaringprobleminvolving some extra flight constraints. Thermal soaring is defined as a constrained optimization problem, and the mathematical formulation is given in the form of a CVRP in Section 107 6.2. In Section 6.3 a flight path decision algorithm is developed. Simulation results for aproblemscenarioarepresentedinSection6.4inordertodemonstratetheefficiencyof theproposedmethodology. Finally,Section6.5drawsconclusions. 6.2 Soaring Assigned Set of Thermals and MathematicalFormulationasaCVRP The thermal soaring problem is adapted to the VRPTW in Chapter 5 where an analogy between the two problems is drawn. Extra constraints called time windows that corre- spond to thermal lifespan periods in soaring applications are considered, and an exact algorithmisdevelopedfortheUAVtoarriveatthedestinationnodeinminimumtime. Inthischapterweremovethethermallifespanconstraintsassumingthatthermalsare available for soaring flight at all times once they are detected. We additionally require that the UAV returns back to where it starts after climbing at all the assigned thermal locations. ThiscanbeinterpretedasservingallcustomersintheVRPformulation. The autonomous thermal soaring is considered as a graph problem as motivated in Chapter5wherethetwometrics,namelythespatialmetric, l ij andthetemporalmetric, c ij aredefinedonthethermalmap. Using the same mathematical relations, optimal arc costs and corresponding hori- zontalflightvelocitiescanbelistedas: c ij = t glide ij +t thermal ij (6.1) = l(i,j) V opt ij + h loss (V opt ij ,V z ij ) W j 108 h loss (V opt ij ,V z ij )=l(i,j)× |V z ij | V opt ij =l(i,j)× −V z ij V opt ij (6.2) V z ij ∼ =aV 2 opt ij +bV opt ij +c (6.3) V opt ij = r c−W j a , c−W j a > 0 (6.4) where l(i,j) istheflightdistancebetweenthermallocations i and j, W j isthethermal strengthatthearrivalnode j,and a, b, c arethecoefficientsofthepolarcurve. Thethermalsoaringproblemona2DgraphG(N,A)canthenbedefinedinitsmore generalformasaCVRP: min X i∈N X j∈N c ij x ij (6.5) subjectto X i∈N x ij = X k∈N x jk ∀ j∈N (6.6) X j∈N x 0j =K (6.7) x ij ∈{0,1} ∀ i, j∈N (6.8) In the above formulation x ij is the flow variable, the value of which is 1 if arc (i,j)∈A is included in the optimal route and 0 otherwise. N is the set of thermals assignedtobesoaredwhere i=0∈N denotesboththestartandtheendlocation. Constraint(6.6)guaranteesthatiftheUAVentersathermalforsoaring,italsoleaves that thermal after its climb. Constraint (6.7) defines the number of tours to be included 109 andisacapacityconstraint. Iftheassignmentiscompletedinasingletourasweassume here, the parameter K is set to 1. A binary constraint is given in (6.8) which is also calledtheintegralityconstraint. InthethermalsoaringproblemscenarioofChapter5becauseoftherelativelysmall number of available thermals, an exact method can be effectively employed to find the optimalpathreachingthedestinationfromthestartnode. In this chapter we consider a specific set of assigned thermals in the flight region. Formulating a VRP as a set-covering problem is suggested in [Balinski and Quandt, 1964], and some exact algorithms are given in [Christofides et al., 1981,Laporte et al., 1986,Agarwaletal.,1989,Hadjiconstantinouetal.,1995]. Ontheotherhand,theVRP in most instances can be solved approximately using near-optimal algorithms. [Lin and Kernighan, 1973,Cullen et al., 1981,Toth and Vigo, 2002] This chapter focuses on dense thermal areas where an approximate solution is acceptable. The details of our routegenerationmethodologyaredescribedinthenextsection. 6.3 AParallelSavingsBasedHeuristicasaSolution Methodology The glider UAV is assigned to visit all thermal locations in the area of interest which mightbeviewedasaset-coveringproblem. Aninstance(N,F)ofaset-coveringproblemconsistsofafiniteset N andafamily F ofsubsetsof N ,suchthateveryelementof N belongstoatleastonesubsetin F asstatedin[Cormenetal.,2001]whereonecanwrite N = [ S∈F S (6.9) 110 A minimum-size set cover can be considered a good match for the thermal soaring problemwhereN isthesetofassignedthermalsonthemapG(N,A).However,weaim toconstructaflightroutewhereeachthermallocationotherthanthestartnodeisvisited only once. We start by dividing the thermal map into nonoverlapping and equal-sized partitionsunlessweknowinadvancethataparticularkindofpartitioningleadstobetter solutions. In order to cover all nodes in these sectors, thermal locations are assigned to initial routes which are defined by applying a minimal spanning tree algorithm. The collection of the routes in all sectors forms a feasible solution for the thermal soaring problem. The problem requirements are then relaxed, and the routes are connected to one another usinga heuristic approachtoselectanear-optimalroute. Thedetailsofthe procedurearedescribedinthefollowingfivesubsections. 6.3.1 Sectorization The flight region is first partitioned into several sections as part of an initialization pro- cedure. In our applications we divide the thermal map into four areas which are East, West,NorthandSouthregions. Itshouldbenotedthatthenumberofinitialroutesinasectordependsonhowdense thethermalsareandhowtheyarelocatedrelativetoeachotherinthatsector. 6.3.2 Kruskal’sAlgorithm Once theflight region isdividedintosectors, aminimalspanningtreeisfoundforeach sector. Aminimalspanningtreeofaweightedgraphisacollectionofedgesconnecting all the vertices on the graph such that the sum of the weights of the edges is at least as small as the sum of the weights of any other collection of edges connecting all the vertices.[Sedgewick,1988] 111 It is explained in [Sedgewick, 1988] that one of the approaches to find the minimal spanning tree is simply to add edges one at a time, at each step using the shortest edge thatdoesnotformacycle. Thealgorithmstartswithan N-treeforest,andfor N steps it combines two trees (using the shortest edge possible) until there is just one tree left. ThisalgorithmisgenerallyattributedtoKruskal. Kruskal’s algorithm computes the minimal spanning tree of a graph in O(ElogE) steps. For a discussion on the running time comparisons of the three minimal spanning treealgorithms,namelythePriority-FirstSearch,Kruskal’sAlgorithmandPrim’sAlgo- rithm, the interested reader might refer to [Sedgewick, 1988]. We implement Kruskal’s Algorithm in MATLAB and employ it for the required minimal spanning tree calcula- tionsinoursimulations. For simplicity we assume that the given map consists of a geographically dispersed set of thermals of the same strength. This assumption results in a symmetric cost def- inition for the edges connecting the nodes on the thermal map which then represents a symmetric weighted graph. Thermals with different climb rates would in general result in flight arcs with asymmetric costs. Such an asymmetric weighted graph is likely to appearinapplications,andtheflightarccostscanstillbedefinedusingtheformulation given in (6.1). However, the solution algorithm we develop in this chapter would then requiresomefurthermodifications. 6.3.3 RouteConstruction Thethermallocationsinthespanningtreeofeachsectorareexamined. Thenodesother thantheonesconnectedtoatleasttwoothernodesinthesamesectorareusedtodefine new edges connecting them to node 0. This step constructs a set of loops on the map whichcollectivelyvisitsallnodes. 112 The Route Construction step in our algorithm presents similarities with the petal heuristics [Renaud et al., 1996] in the sense that there are possibly more than one route in any sector. However, we do not generate a fixed number of petals for each sector, andthenumberofroutesratherdependsonthestructureofthespanningtreedefinedby Kruskal’sAlgorithmappliedinthepreviousstep. 6.3.4 RelaxationandParallelSavingsBasedHeuristics We apply a Relaxation step such that the aircraft does not necessarily visit the thermal locationat node 0,whereas i ∈N−{0}representsthenodesrequiredtobeincluded in the route of the UAV. We involve a parallel savings criterion in order to modify and mergetheconstructedroutesintoasingleflightpathwhichinturndefinesanear-optimal solutionforthisparticularrelaxedproblem. A survey of the classical and modern heuristics for the VRP is given in [Laporte et al., 2000]. The classical heuristics for the CVRP are classified in [Toth and Vigo, 2002] as Constructive Methods, Two-Phase Methods, and Improvement Heuristics. Constructive Methods are further classifed as Clarke and Wright Savings Algorithm [Clarke and Wright, 1964], Matching-Based Savings Algorithms [Desrochers and Ver- hoog, 1989,Altinkemer and Gavish, 1991,Wark and Holt, 1994], and Sequential Inser- tionHeuristics[MoleandJameson,1976,Christofidesetal.,1979]. We inspire from the Clarke and Wright Savings Algorithm and apply a heuristic method similar to its parallel version. Furthermore, we incorporate some additional but simpleimprovementstoenhancethequalityoftheresultingsolution. A parallel and a sequential version of the Clarke and Wright Savings Algorithm are available. The parallel version of the algorithm, as discussed in [Toth and Vigo, 2002], proposesthebestfeasiblemergesuchthatthesesavingsarecheckedinanonincreasing fashion, and given a saving s ij , it can be determined whether there exist two routes, 113 onecontaininganedge(0,j)andtheothercontaininganedge(i,0)thatcanfeasiblybe merged. Thesetworoutesarethencombinedbydeletingtheedges(0,j)and(i,0),and addingtheedge(i,j). Ifthetwonodes i and j areconnected,thesaving s ij canbecomputedas: s ij =c i0 +c 0j −c ij , i6=j, i, j∈S (6.10) whereS inourcaseisthesetofnodesdirectlyconnectedtonode 0 uponcompletionof theRouteConstructionstep. Thenodes i and j in(6.10)arealsoimpliedtobenodes indifferentinitialrouteswhicharenotconnectedotherthanatnode 0. In order to merge the initial routes using such a savings criterion, the edges (i,0) and (0,j) are interchanged with (i,j), and the costs c i0 and c 0j are saved no matter whichsetofnewedges(i,j)theParallelSavingsBasedHeuristicsstepintroduces. The performanceofthesavingsalgorithmisthusdeterminedsolelybythecosts c ij ,where i6=j, i, j∈S, i and j arenodesindifferentinitialroutes. Our savings strategy reduces to pairing these nodes i and j such that the sum of the costs c ij is minimum. It should be noted that the heuristic step in our algorithm is not necessarily repeated for all nodes on the map, i.e., we do not calculate the savings forallnodepairs i and j where i, j ∈N , i6=j. Afurtherimprovementinthealgorithmcanbestraightforwardlyachievedincasethe solution for the relaxed problem is self-intersecting. The intersecting edges determine thenodeswhichmightpossiblybereconnectedinanotherorderresultinginaflightcost reduction. Wealsocheckthenodesdefiningtheseself-intersectionsalongwiththeothernodes visitedbetweenthemandevaluateallpossibleconnectionsinordertoconcurrentlyreop- timizetheroutesproducedintheRouteConstructionphase. 114 Althoughthisimprovementprocessisnotasophisticatedneighborhoodsearchrule, it can be implemented with a low calculation cost. Since it is basically a deeper explo- ration of the most promising regions of the solution space, the solution of the problem islikelytobeimproved. 6.3.5 Near-OptimalRouteSelection Node 0 is not included in the route determined by the Parallel Savings Based Heuris- tics due to the Relaxation phase involved prior to applying the heuristic step. Feasible solutionstotheoriginalproblemcanbegeneratedsimplybyinsertingnode 0 between anytwonodesintheroutedeterminedusingtheParallelSavingsBasedHeuristics. Let j = Δ + (i) be the node that follows node i along the route proposed for the relaxed problem. The nodes i and j can be disconnected by deleting the edge (i,j). If node 0 is then inserted between these two nodes, two new edges(i,0) and(0,j) are addedtotheroute. While guaranteeing the transformation of the solution for the relaxed problem to a feasibleroutefortheoriginalproblem,thisinsertionproducesanextracostof e ij which isgivenby e ij =c 0i +c j0 −c ij , j = Δ + (i), i, j ∈N , i6=j, i6= 0, j6= 0 (6.11) The two nodes i and j are thus chosen to minimize the value of e ij defined in (6.11). The resulting route is feasible and is an approximate solution to the optimal soaringproblem. Sincethethermalmapunderconsiderationisnotadirectedgraph,visitingthesame nodes in the reverse order results in the same overall cost, and is not considered as a differentsolutionforthescenariosunderdiscussion. 115 The thermal map is initially divided into sectors with no a priori information on the solutionoftheproblem. Thealgorithmthereforemightnotbeabletoreachtheoptimal solutionandisingeneralanapproximatemethodonly. However,thepathsinthesectors allocating thermal locations are eventually modified and merged so that the sectors are notisolated. Anear-optimalsolutionisreachedwhichrepresentsanupperboundonthe costoftheoptimalroute. Yet,theproposedmethodologyiscomposedofasetofsimple steps,andissignificantlyefficientintermsofcomputationsinvolved. Our solution algorithm might reach a set of solutions in the end if the overall costs determinedfortheseroutesareequallylow. 6.4 Simulations A specific problem scenario is chosen to illustrate the performance of our approximate algorithm. A total of twenty thermals are detected to be available for soaring within the flight range including the one located at node 0 where the glider UAV starts. The aircraftisassignedtovisitallthermallocationsandreturnbacktonode 0 inminimum possibletime. The maximum flight arc distance is defined as l max =20,000ft, and the thermals locatedfurtherawayfromeachotherarenotconnectedbyadirectarconthemap. Each thermalisassumedtohaveastrengthof W=10ft/sec definingtheoptimalhorizontal flight velocity of the glider as V opt =204ft/sec. The problem scenario is summarized inFigure6.1. Each flight arc on the map is associated with a cost using the analytical expression givenin(6.1). ThenumericalvaluesofthetimecostsareevaluatedforallarcsinFigure 6.1 and listed in Table 6.1 where the calculated values, t calc are compared with the measurements, t meas fromsimulations. Theerrorinbetweenthetwosetsisvastlydue 116 Figure 6.1: The thermal map with the direct flight distances between the thermal loca- tionsdepictedsothatthetimecostsofthecorrespondingflightarcscanbelisted(Section 6.4) to the polar curve approximation used to evaluate t calc and the transients observed in thesimulationswhichinturnaffectthevalueof t meas asdiscussedinChapter5. Theflightregionisdividedintosmallerpartitionstoinitializeoursolutionalgorithm, andtheresultingsectorsareshowninFigure6.2. Theminimalspanningtreesarefound forallfoursectorsseparatelyandaredemonstratedbysinglesolidlinesonthemap. The doublesolidlinesinthesamefigurerepresenttheedgesaddedtoconnecttheendpoints ofthesespanningtreestonode 0.Thesolidsingleanddoublelinestogetherdisplaythe initialflightroutes. 117 Table6.1: FlightTimeTable: SimulationScenarioinSection6.4 W=10ft/sec and V opt =204ft/sec Distance t glide t thermal t calc t meas (ft) (sec) (sec) (sec) (sec) 4,000 19.6 33.9 53.5 49.9 5,000 24.5 42.3 66.8 61.7 6,000 29.4 50.8 80.2 71.7 7,000 34.3 59.3 93.6 85.6 8,000 39.2 67.7 106.9 98.5 9,000 44.1 76.2 120.3 109.6 10,000 49.0 84.7 133.7 125.8 11,000 53.9 93.1 147.0 138.2 12,000 58.8 101.6 160.4 148.2 13,000 63.7 110.1 173.8 161.2 14,000 68.6 118.5 187.1 178.4 15,000 73.5 127.0 200.5 189.4 16,000 78.4 135.4 213.9 201.3 17,000 83.3 143.9 227.2 209.7 18,000 88.2 152.4 240.6 224.0 19,000 93.1 160.8 253.9 241.4 20,000 98.0 169.3 267.3 254.2 A total of seven routes are determined where the North routes are 0−2−1−0 and 0−2−3−4−0,theSouthroutesare 0−10−11−9−0 and 0−12−16−15−0,theEastroute is 0−7−6−5−8−0,andtheWestroutesare 0−17−19−18−0 and 0−17−13−14−0. The East route 0−7−6−5−8−0 has two intersecting edges which are (7,6) and (8,0), whereas the same situation holds for the edges (17,19) and (18,0) of the West route 0−17−19−18−0. All possible combinations for reconnecting the set of nodes {0,5,6,7,8} fortheEastrouteand {0,17,18,19} fortheWestrouteareexamined. As a result, the East route is improved into 0−6−5−8−7−0 which provides the mostsavings,whereastheWestrouteisshowntobeanoptimalconnectionofthenodes involvedandhencekeptunchanged. TheRouteConstructionstepiscompleted. 118 Figure6.2: Thesolutionalgorithmiteratedonthethermalmap(Section6.4) In the Relaxation step the edges connecting the nodes { 1,2,4,6,7,9,10,12,14, 15,17,18} tonode 0 aredeleted. Parallelsavingsareconsidered,andtheedges(4,6), (7,10), (9,15), (12,14) and (1,18) are added. Note that the two routes in the North remainconnectedatnode 2 althoughtheedge(0,2)isdeleted. Thesameholdsforthe tworoutesintheWestatnode 17.Thetwoedges(10,11)and(9,15)canbeshowntobe intersecting. Afterrelevantcalculationsthesetofnodes{9,10,11,15} arereconnected as 10−9−11−15. 119 Figure6.3: Anear-optimalflightpathreachedusingtheproposedsolutionmethodology (Section6.4) In the Near-Optimal Route Selection step two possible solutions are determined to be S1 and S2 where S1 : 0−8−5−6−4−3−2−1−18−19−17−13−14−12−16−15−11−9−10−7−0 S2 : 0−7−8−5−6−4−3−2−1−18−19−17−13−14−12−16−15−11−9−10−0 The solution S1 is demonstrated in Figure 6.3. For this specific route the overall flight time is calculated to be t calc =2,218.9sec, whereas the total cost is measured throughsimulationsas t meas =2,044.1sec. 120 6.5 Conclusion We propose an approximate algorithm based on the VRP interpretation of the thermal soaringproblemforanautonomousUAV.Weconsiderproblemscenarioswhereanexact solution algorithm appears to be demanding in computational time. Assuming that a near-optimal solution is acceptable, a methodology involving Sectorization, Kruskal’s AlgorithmandParallelSavingsBasedHeuristicsisdeveloped. The simulation results confirm the fact that the physical requirements and the flight constraintsofsoaringUAVscanbeproperlytranslatedintomathematicallanguage,and efficientalgorithmscanbedevelopedtodefinetheflightpath. Similar heuristic methods can also be useful for powered UAVs in search flight. Our methodology is promising for surveillance purposes where energy requirements inevitablydeterminetheflightperformanceoftheaircraft. 121 7. GENETICALGORITHMSFOR SHORTESTPATHROUTING OFAUTONOMOUSGLIDERS 7.1 Introduction Several criteria are established in [Garey and Johnson, 1979] in order to measure the extent to which a heuristic’s solution is near to or far from the optimal one, and as discussedin[Falkenauer,1998]thesecriteriaareusuallycastasperformanceguarantees inthattheymeasurethequalityofthesolutionfortheinstanceswhicharethehardestto solvebytheheuristicmethod. GAsaredefinedin[Goldberg,2002]assearchproceduresbasedonthemechanicsof natural selection and genetic. They are heuristic search techniques which demonstrate highperformancequalitywhenproperlyemployed. The research studies on GAs were initiated in 1975 through [Holland, 1975] where Holland presents the theoretical foundations and explores applications of the field. As discussed in [Man et al., 1997] another quantum leap in GA research and development was witnessed in 1990s by academics and engineers worldwide, the reason of which is explained through the advances in solid-state microelectronics fabrication, and the correspondingcostreductionandimprovementsincomputertechnology. 122 A GA is a randomized global search technique used to solve problems by imitat- ing processes observed during natural evolution where a population of chromosomes is evolved by mimicking natural phenomena. [Toth and Vigo, 2002] The mechanism of generationofnewcreatureshighlyskillfullyexploitstheknowledgeaccumulatedinthe current pool of living organisms in order to produce new ones as capable of surviving as their ancestors. [Falkenauer, 1998] A set of possible feasible solutions for a problem canbeinterpretedasapopulationofindividualsrepresentedbychromosomes. For an optimization problem the objective function takes its optimal value for the feasible solution which is called an optimal solution. This function can be defined in termsofthefitnessvaluewhichreflectsthequalityofthechromosomes. Theideaofthe survival of the fittest from theevolutionarytheorythusleadstotheoptimalsolutionfor theproblemunderconsideration. TodayGAshaveawiderangeofapplicationsinverydiversefields. Inanotherclas- sic work of Holland’s [Holland, 1992] a mathematical model for complex systems is presented and applied to economics, physiological psychology, game theory, and arti- ficial intelligence. The GA applications in the areas of control and signal processing receive specific attention in [Man et al., 1997] where some advanced GA applications suchasGAsintimedelayestimation,GAsinactivenoisecontrol,andGAsinautomatic speechrecognitionarestudied. TheGAsareappliedtorobottrajectorygenerationprob- lemsin[Davidor,1991]. MorerecentworkontheuseofGAsinvehicleautomationcan be found in [Hernandez et al., 1995,Marin et al., 1999,Nikolos et al., 2003,Alvarez etal.,2007]. In this chapter we propose that GAs can also be incorporated as optimization tools forautonomoussoaringflightapplications. Thetraditionalterminologyforconventional GAs is given in Section 7.2. A specific thermal soaring problem with no temporal constraintsisinvestigatedinSection7.3wheretheformulationasaVRPisrevisited. A 123 GA-based SP algorithm is applied to choose the soaring flight path, and the details of thealgorithmaredescribedinSection7.4. AsimulationscenarioispresentedinSection 7.5 where the result of the employed GA is compared to the optimal solution reached by an exact algorithm, demonstrating the efficiency of the formulation and the solution strategy. OurconclusionsaresummarizedinSection7.6. 7.2 OptimizationinTermsofBiologicalEvolution InGAapplicationsthesearchspaceofpossiblesolutionsfortheproblemismappedonto a set of finite strings called chromosomes. [Falkenauer, 1998] An initial population of thechromosomesrepresentingtheindividualsisselected,andafitnessvalueisevaluated for each individual in the population. Depending on their fitness values, the individuals are selected to form the parent set. The selection process is noisy in the sense that the selected individuals are not necessarily the fittest ones. Still, the probability of an individualtobeselectedishigherforabetterfitnessvalue. Thechromosomesintheparentsetarecrossedovertoproducenewindividuals. The less fit individuals in the population are replaced with these new chromosomes where theselectionoflessfitindividualsisalsoanoisyprocess. In order to increase the variety in the population, the resulting set of chromosomes allows mutations which are small random changes on a randomly selected set. In some GAapplications,asdiscussedin[Falkenauer,1998],asmallrandomlychosenportionof thepopulationisalsosubjecttoanothergeneticoperatorcalledinversionwhichchanges thepositionsofthegenesonthechromosomes. Onceanewpopulationisgenerated,theterminationconditionoftheGAischecked. IftheGAisnotterminated,itstartsworkinginparallelonthenewsetofchromosomes again for the next iteration of the algorithm. The termination condition can be defined 124 asalimitonthenumberofallowedgenerationsoracriterionforthefitnesslevelofthe currentpopulation. The design of GAs, as discussed in [Goldberg, 2002], requires both clarity of aims andtheapplicationofeffectivedesignmethods,whereastheGAswhichcansolvehard problemsquickly,accuratelyandreliablyarecalledcompetent. WefocusontheautomationofagliderUAVinthermalsoaringapplicationsandfor- mulate the problem as a VRP involving several flight constraints which can be reduced to an SP routing problem. A genetic algorithmic approach suitable for solving SP rout- ingproblemsisthenfollowedasareliabletrajectorygenerationmethodology. 7.3 Mathematical Formulation of Thermal Soaring Problem: TheVRPRevisited ThetwodimensionalgraphG(N,A)isusedtodescribethethermalmapasinChapter5 andChapter6where N isthesetofthermallocationsand A isthesetofstraightflight arcs between these thermals. We use the same metrics and pose the same restrictions on the direct flight range in between thermal locations i and j where each arc cost c ij on the graph represents the particular time required to include that arc in the flight route, involving both the glide time along the arc (i,j) and the climb time at the arrival thermallocation, j. Usingthesamemathematicalrelations,theoptimalarccostsc ij andthecorrespond- ing horizontal flight velocities V opt ij can be listed. In order to formulate the optimal static soaring problem where the UAV starts at a node and is guided towards another, thebasicmodelfortheVRPcanthusbemodifiedasfollows: 125 min X i∈N X j∈N c ij x ij (7.1) subjectto X i∈M−{n} x ij = 1 ∀j∈M−{0} (7.2) X j∈M−{0} x ij = 1 ∀i∈M−{n} (7.3) X j∈N x 1j = 1 (7.4) X i∈N x in = 1 (7.5) x ij = 1, ∀ i, j∈M, j = Δ + (i) 0, else (7.6) where x ij istheflowvariableofvalue 1 ifthearc(i,j)∈A isintheoptimalrouteand 0 otherwise. Along the optimal route the node Δ + (i) represents the node that follows node i.Thesetofthermalssoaredaredenotedby M−{0} where i=0∈M isthe startnode,and i=n∈M isthedestination. ThisproblemformulationisdifferentfromtheoneinChapter6inthesensethatwe nomorediscussasetcoveringproblem,butratherfocusonflightfromthestartpointto the destination across a region which is dense with thermals where there is no specific set of thermal locations to be visited other than the start and the destination nodes. The problem definition in this chapter is thus similar to the one in Chapter 5 except the fact that the thermal lifespan constraints are ignored. This is in accordance with the assumption of Chapter 6 that thermals are available for soaring flight at all times once theyaredetected. 126 In Chapter 5 we employ the Floyd-Warshall Algorithm which can solve a given SP routing problem to optimality. Since we consider much denser thermal maps in this chapter, the solutions are encouraged to be computed using algorithms which possibly resultinreducedcalculationtime. Weemployaroutegenerationmethodologybasedon GAs assuming that an approximate solution is acceptable. The details of the methodol- ogyaredescribedinthenextsection. 7.4 A GA-Based SP Algorithm Adapted to Thermal SoaringProblem Several VRPs are solved using GAs in [Grefenstette et al., 1985,Blanton and Wain- wright,1993,Hwang,2002]. Itisdiscussedin[GillandZomaya,1998]thatGAsarealso highlysuitedtosolvingpathplanningproblemswherereasonablesolutionsarerequired to be obtained within short computational time. Different applications of GAs to the SProutingproblemcanbefoundin[Munetomoetal.,1998,Leungetal.,1998,Inagaki et al., 1999,Ahn and Ramakrishna, 2002]. We follow the algorithm proposed in [Ahn and Ramakrishna, 2002] since it promises high quality solutions and an enhanced rate ofconvergence. Forpathplanningpurposestheroutesfromthestartnodetothedestinationareused to define the chromosomes of the GA. The genes of the chromosomes are the variables denotingthenodesonthemapwhicharevisitedalongthatroute. Alocusistheposition of a gene on the chromosome. The first gene on a chromosome defining a feasible path is 0 whichisthestartnode. Thesecondgenecanbeanyofthenodesdirectlyconnected to node 0, and so on. The last gene of each chromosome representing a feasible route is the destination node, n, and has the locus m for a chromosome of m genes. The length of a chromosome is limited in our case by the total number of thermal locations 127 detected in the flight region since the loops are necessarily avoided in a shortest path routingscheme. The GA in [Ahn and Ramakrishna, 2002] uses variable length chromosomes and employsacrossoveroperatoratpositionallyindependentcrossingsites. Anoverviewof themainphasesofthisGAisgiveninthefollowingsubsections. 7.4.1 PopulationInitialization Therearedifferentmethodsthatcanbeusedtoinitializethepopulationofchromosomes in GAs, whereas the one proposed in [Ahn and Ramakrishna, 2002] and used in this chapterisrandominitialization. Thepopulationsizeistakentobethenumberofnodes onthethermalmapinourapplication. Thereaderisreferredto[AhnandRamakrishna, 2002]foradiscussiononpopulationsizeanditsinfluenceonqualityofsolution. 7.4.2 Selection(Reproduction) A fitness function that can measure the quality of chromosomes in the population is definedin[AhnandRamakrishna,2002]as: f r = 1 P lr−1 k=1 c gr(k)gr(k+1) (7.7) where f r denotes the fitness value of the r th chromosome which has a length l r . In (7.7) g r (k) is the gene in the k th locus of the r th chromosome, or equivalently the gene representing the k th node in the route number r with l r −2 other thermal locations in between the nodes 0 and n. By definition g r (1)=0 and g r (l r )=n for anyhealthychromosomeinthepopulation. In order to focus on the most promising regions in the solution space two chromo- somes are selected from the population of chromosomes, and a pairwise tournament 128 selection without replacement is applied where the fitter among the two chromosomes are selected from each pair by comparing their fitness values through (7.7). A chromo- someselectedonceisnotaddedbacktothepopulation,andthusitcannotbechosenas aparenttwice. 7.4.3 Crossover(Recombination) ThecrossoveroperatorintheGAcombinestogetherpiecesofinformationcomingfrom different individuals in the population [Falkenauer, 1998]. In the proposed scheme of [Ahn and Ramakrishna, 2002] two chromosomes chosen for crossover should have at least one common gene except for the source and destination nodes, but the gene is not necessarily located at the same locus of the two chromosomes, thus the crossing points of the two chromosomes are not necessarily the same. In the partial route exchange in between these chromosomes a restriction can be included such that one partial route connectstonode 0,whereastheotheronecrossoveredconnectstonode n.Thishasthe interpretationofgeneratingoffspringsfromparents. Notethatfromeachpairofparents, apairofoffspringsarecreated. 7.4.4 Mutation Mutation can be applied to modify the generated offsprings and to maintain genetic diversityofthepopulationwhileavoidinglocaloptima. If a chromosome undergoes mutation, one of its genes is randomly selected as the mutation point. An alternative partial route is then generated by randomly connecting thenodeatthemutationpointtothedestinationwithnonoderepetitionsallowed. ThemutationprobabilityintheproposedGAissettoatypicalvalueof 0.05.Once the population of chromosomes converges, however, mutation is not applied to avoid uninvitedevolution. 129 7.4.5 RepairFunction The crossover might generateloops andthus infeasiblechromosomes whichneed tobe repaired. Asimplerepairfunctionisusedin[AhnandRamakrishna,2002]toeliminate the loops in the generated route which corresponds to eliminating lethal genes from a chromosomeandturningitintoahealthyone. Thelessfitchromosomesdeterminedintheselectionphasearethenreplacedbythe offsprings, and the algorithm is iterated until the population convergence is achieved. For more details on the described GA, the interested reader might refer to [Ahn and Ramakrishna,2002]. 7.5 Simulations Wechooseaspecificproblemscenariowith20thermallocationsdetectedtobeavailable forsoaringwithintheoverallflightrangeincludingthestartnode 0 andthedestination node 19. The aircraft is assigned to define its path from node 0 to node 19 possibly visiting several thermal locations such that the total flight time is minimized, and thus thecross-countryspeedofthevehicleismaximized. The maximum flight arc distance is defined as l max = 2,000ft to determine the thermal locations which are connected by a direct flight arc. If a strength of W = 10ft/sec is detected for each thermal, the optimal horizontal flight velocity of the glider is calculated as V opt =204ft/sec along the inter-thermal flight arcs. Note that theterrainfeaturesdonotnecessarilyallowtheUAVtoflyfromonethermallocationto thenextalthoughitiswithin 2,000ft flightdistance,andthusnotallthenodeswithin thedefinedmaximalrangearedirectlyconnected. In a more realistic problem scenario the UAV might also be required to fly through downdrafts in between detected thermal locations where the air currents in fact cause 130 Figure 7.1: The feasible direct flight distances between the thermal locations shown on the thermal map so that the time costs of the respective flight arcs can be calculated (Section7.5) the vehicle to lose altitude. We consider various downdrafts of strength in a range of −1ft/sec to −15ft/sec, and they are also denoted on the relevant thermal map in Figure 7.1. The circles are the nodes of the graph representing the thermal locations, whereastrianglesareusedtolocatethedowndrafts,thesinkratesofwhicharespecified onthethermalmapin ft/sec. The directly connected thermal locations are linked with arcs satisfying the triangle inequalityforflightdistance,however,thelengthofthearcsisnotnecessarilyscaledin thefigure. 131 In order to find the optimal path connecting node 0 to node 19, each flight arc on themapisassignedtohaveacostwhichisderivedinpreviouschapters. Intherelevant analyticalexpression,however,theeffectsofthedowndraftsareyettobeincluded. The respective numerical values are evaluated for all arcs in Figure 7.1, and are listed in Table7.1andTable7.2. Foreachdirectflightdistancerepresentedbyanarc, t glide istheglidetimerequired to cover the given distance, t thermal is the thermal climb time at the arrival node, and t downdraft is the time cost of the downdrafts located along that flight arc. Note that a downdraft sinking at a rate of 1ft/sec is said to have a strength of −1ft/sec and is assumed to result in an increase of 10sec in the cost of the corresponding flight arc. ThenumericalvalueofCostisthusthetotalcostofthearcincludingtheeffectsof downdraftsaswell. IftheflightdistanceandthedowndraftdatagiveninFigure7.1arereplacedwiththe costs of the arcs, the map is shifted to the one in Figure 7.2. The problem scenario thus reducestothenetworkstructuredepictedin[AhnandRamakrishna,2002]. The size of the population for the GAs is chosen as the number of nodes in the network,whichis n+1=20.TheGAin[AhnandRamakrishna,2002]thenproposes theflightpathcorrespondingtotheroute S 1 : 0−2−7−13−19 (7.8) whichhasatotalcostof 142sec. As discussed in [Ahn and Ramakrishna, 2002] for this particular problem scenario Munetomo’sAlgorithm[Munetomoetal.,1998]wouldsettlefortheroute S 2 : 0−1−2−7−13−19 (7.9) 132 Table7.1: FlightTimeTable(I):SimulationScenarioinSection7.5 W=10ft/sec and V opt =204ft/sec Distance t glide t thermal t downdraft Cost (ft) (sec) (sec) (sec) (sec) 156 0.8 13.2 0 14 168 0.8 14.2 0 15 179 0.9 15.1 0 16 190 0.9 16.1 0 17 223 1.1 18.9 0 20 246 1.2 20.8 0 22 268 1.3 22.7 0 24 279 1.4 23.6 0 25 290 1.4 24.6 0 26 302 1.5 25.5 0 27 324 1.6 27.4 0 29 335 1.6 28.4 0 30 357 1.8 30.2 0 32 357 1.8 30.2 40 72 391 1.9 33.1 0 35 424 2.1 35.9 50 88 447 2.2 37.8 0 40 480 2.4 40.6 0 43 503 2.5 42.5 0 45 559 2.7 47.3 0 50 559 2.7 47.3 200 250 570 2.8 48.2 10 61 581 2.8 49.2 10 62 603 3.0 51.0 0 54 648 3.2 54.8 0 58 670 3.3 56.7 0 60 670 3.3 56.7 160 220 681 3.4 57.6 0 61 726 3.5 61.5 0 65 782 3.8 66.2 20 90 782 3.8 66.2 80 150 793 3.9 67.1 0 71 133 Table7.2: FlightTimeTable(II):SimulationScenarioinSection7.5 W=10ft/sec and V opt =204ft/sec Distance t glide t thermal t downdraft Cost (ft) (sec) (sec) (sec) (sec) 804 3.9 68.1 0 72 849 4.1 71.9 60 136 860 4.2 72.8 0 77 894 4.4 75.6 30 110 894 0.9 15.1 150 230 994 4.9 84.1 0 89 1005 4.9 85.1 40 130 1050 5.1 88.9 50 144 1050 5.1 88.9 100 194 1307 6.4 110.6 30 147 1798 8.8 152.2 0 161 2000 9.8 169.2 0 179 whereasInagaki’sAlgorithm[Inagakietal.,1999]wouldproposetheflightpath S 3 : 0−3−8−13−19 (7.10) whichwouldcost 187sec and 234sec respectivelyinourcasefortherelevantthermal soaringproblem. One might refer to [Ahn and Ramakrishna, 2002] for a comparison of the perfor- mance of the proposed GA to Munetomo’s Algorithm and Inagaki’s Algorithm which settle only for suboptimal paths. These algorithms are also shown to have higher route failureratios. In order to verify the performance of the GAs discussed in this section, we also compare the resulting flight path to the one iterated by an exact algorithm. We employ Floyd-Warshall Algorithm as in Chapter 5 to obtain the optimal solution by carrying out as many steps of iterations as the number of nodes. Iterating the Floyd-Warshall 134 Figure 7.2: The time costs of the feasible flight arcs displayed on the thermal map (Section7.5) Algorithm n + 1=20 times, the final cost matrix, C f =C (20) and the predecessor matrix, Π f =Π (20) definingthevisitednodesareobtained. Theoptimalrouteis S 4 : 1−3−8−14−20 (7.11) whichisexactlythesamerouteas S 1 withatotalcostof 142sec,andisdemonstrated inFigure7.3. TheGA-basedSPalgorithmisthusshowntoconvergetotheoptimalflightpath. Itis claimedin[AhnandRamakrishna,2002]thattheproposedGAalgorithmisinsensitive 135 Figure7.3: TheoptimalsolutionobtainedbyiteratingtheGA(Section7.5) to variations in network topologies in respect of both route optimality and convergence speed. For regions with much denser thermal locations the GAs can then be employed with no extra computational time requirements, whereas the exact algorithms become less and less useful. The employed GA-based SP algorithm is hence a reliable method- ologythatcanbeusedtoguideUAVsacrossgiventhermalregionsingeneral. TheGAemployedforthethermalsoaringprobleminthischaptercanbemodifiedto provide an approximate solution to the VRPTW formulation considering thermals with limited lifespan as previously discussed. The GA can replace the Route Optimization step of the exact algorithm proposed in Chapter 5 and used together with the proposed 136 PreprocessingandRouteValidationstepswherethetemporalconstraintsarealsoincor- porated. 7.6 Conclusion We reduce the VRP interpretation of a thermal soaring problem to the corresponding SP routing formulation where a route planning algorithm can be developed based on the natural evolution paradigm. In simulations the path proposed by the GA-based SP algorithmcanbeshowntoconvergewithinreasonablecomputationaltimetotheoptimal solutioncalculatedusingtheFloyd-WarshallAlgorithm. Inthermalsoaringapplicationsadditionalthermalsmightbedetectedoncethevehi- cle starts navigating across the thermal region, or the already detected thermals might be no longer available for soaring flight due to unexpected changes in their lifespan. In order to achieve optimal flight performance in dynamic environments, online compu- tations based on the most current thermal data are desired although they increase the computational load. However, the computational time of the proposed GA is reported nottoincreaseextensivelywiththenetworksizeunlikecorrespondingexactalgorithms whichmightnotsupportthereal-timerequirementsforlargenetworksduetotheirpro- hibitivecomputationalcosts. TheGA-basedSPalgorithmcanalsobeemployedfordifferentnetworktopologies. It can thus be used to guide UAVs between designated start and destination locations across dense thermal regions with different thermal location distributions. Although the algorithm does not guarantee convergence to the optimal path in general, the route failureratiocanbeshowntoberelativelysmallincomparisontoseveralotherheuristic 137 searchalgorithms. Consequently,wesuggestinthischapterthattheGAunderconsider- ationcanbeimplementedforeffectiveuseinintelligentflightpathplanningapplications soastocontributetotheautomationofUAVs. 138 8. STABILITYANALYSISFOR INDIRECTADAPTIVECONTROLDESIGN WITHANTI-WINDUP COMPENSATION 8.1 Introduction Themainchallengewithalinearcontrolsystemsubjecttoinputsaturationistomaintain stabilityasthesystemgoesintothesaturationphase. Implementingaproperlydesigned anti-windup compensator guarantees a swift recovery of the linear behavior so that the original control objective canbeachievedwithnosignificantperformancedegradation. The anti-windup compensation or conditioning is essentially an a posteriori approach where an additional compensator is designed to cope with control constraints. [Grimm etal.,2001]Suchacompensatordesign,however,requiresingeneralaplantmodelwith aconsiderabledegreeofaccuracy. Several approaches for the control of uncertain systems with bounded inputs are investigatedin[TarbouriechandGarcia,1997]includingdynamicoutputfeedbackcom- pensation and piecewise-linear LQ control design techniques. Robust anti-windup con- trol design methods have also been recently considered by the anti-windup community. Different configurations for anti-windup control of uncertain systems are summarized in [Villota et al., 2006], and uncertainty is discussed in LMI based anti-windup syn- thesis for continuous-time plants in [Grimm et al., 2002,Turner et al., 2004,Wu and 139 Jayasuriya,2006]. Anonlinearschedulingtechniquewithhysteresisswitchingamonga family of linear gains for anti-windup construction is proposed in [Zaccarian and Teel, 2004]. LinearParameterVarying(LPV)controltheoryisappliedtodesignanti-windup controllers for parameter varying plants with input constraints in [Wu and Grigoriadis, 1999]. Simulation results in [Lu et al., 2005] confirm that an LPV anti-windup com- pensator yields enhanced reliability and an expanded envelope for aircraft under large maneuver operations. The stability of systems involving LMI based anti-windup aug- mentation is still an open issue for plants with unknown parameters which we aim to addressinthischapter. The adaptive control design with anti-windup compensator is applied to control a glider UAV model with unknown parameters in Chapter 3. In this chapter we provide thestabilityanalysisasdevelopedin[KahveciandIoannou,2007]. We introduce some basic definitions in Section 8.2. A nominal LQ control design with plant-order anti-windup augmentation for the tracking performance of an asymp- toticallystableplantisdescribedindetailforcompletenessinSection8.3. Plantparam- eters are assumed to be unknowninSection8.4wherean indirectadaptivecontrolleris proposed. The stability properties of the design are developed in Section 8.5. Conclu- sionsaregiveninSection8.6. 8.2 Preliminaries The same notation as in [Turner et al., 2004] is used for basic norm definitions. The induced L 2 normorthefinite L 2 gainofanoperator H isdefinedas: kHk i,2 = sup 06=x∈L 2 kHxk 2 kxk 2 (8.1) 140 wherekxk 2 = q R ∞ 0 kxk 2 dt istheL 2 normofthevectorx(t),andkxk isitsEuclidean norm. If kxk 2 isfinite, x∈L 2 . Wedefinetheexponentiallyweighted L 2 normas: kx t k 2δ = Z t 0 e −δ(t−τ) x T (τ)x(τ)dτ 1 2 (8.2) where δ≥ 0 isaconstant. Similarly,if kx t k 2δ isfinite, x∈L 2δ . Thedefinitionof L ∞ normisgivenby kxk ∞ = sup t≥0 |x(t)| (8.3) where x∈L ∞ if kxk ∞ exists. A decentralized saturation function for a vector u ∈ R nu is described in [Grimm etal.,2003]as: sat(u) = [sat 1 (u 1 ) sat 2 (u 2 ) ... sat nu (u nu )] T (8.4) sat i (u i ) = u i max n 1, |u i | M i o (8.5) where M i ∈R, M i > 0 for i = 1,...,n u .Wefocusonascalarinput u∈R anduse sat(u) =sign(u)min(|u|,¯ u), ¯ u∈R, ¯ u> 0 (8.6) for the saturation function where ¯ u and −¯ u are the upper and the lower saturation limitsrespectively. 141 8.3 NominalControlDesign We consider a linear plant with a saturation nonlinearity imposed on the control input u act whichcanbedescribedusing ˙ x =A p x+B p sat(u act ) (8.7) y =C p x If it is assumed that all plant states in (8.7) are available for measurement, we can define C p as the identity matrix. Note that we also assume there is no direct path of nonzero gain between the input and the output which is a valid assumption for most physicalsystemsasdiscussedin[IoannouandSun,1996]. Thecontrolobjectiveistoregulatethetrackingerror e =x−r (8.8) tozerowhere r isaconstantreferenceinput. Wechoosetoapplyatwo-stepdesignpro- cedure where the saturation limits are first ignored and the control for the linear system isdeveloped,andthisinitialdesignisthenaugmentedwithananti-windupcompensator. Our control design represents the transfer function matrix from −e = (r−x) to u, where u is the output signal of the controller. If system matrices are accurately knownandaretimeinvariant,wecanbenefitfromthestraightforwardnessinsystematic implementationofanLQcontroldesign. Itcanbeshownthatiftheplantiscontrollable andobservable,thefollowingARE: A T aug P +PA aug +Q z −PB aug R −1 z B T aug P = 0, (8.9) 142 A aug = A p 0 I 0 , B aug = B p 0 , Q z = 0 0 0 I ×Q, Q> 0, R z > 0 canbesolvedtoobtainthecontrolparameters K = h K 1 K 2 i where K =R −1 z B T aug P (8.10) TheresultisaPIcontrollerintheform: u =−K 1 e−K 2 Z t 0 e(τ)dτ (8.11) andapossiblerealizationofthecontrolleris ˙ x c =A c x c +B c (r−x) (8.12) u =C c x c +D c (r−x) with A c = 0, B c =I, C c =K 2 , D c =K 1 . Thetrackingproblemisthussolvedforalinearplantwithnoinputconstraintswhere it is assumed that the output of the controller, u can be directly applied to the control input of the plant. In applications involving the saturation nonlinearity, however, the correspondinglimitsneedtobeconsidered. Havingthecontrollerforthelinearplantathand,wedefinethe L 2 gain γ fromthe externalreferencesignal r toourperformanceoutput e,whichwouldsatisfy kek 2 ≤γkrk 2 (8.13) 143 The goal of the anti-windup design in [Grimm et al., 2003] is to construct an anti- windup compensator that guarantees a performance level γ as small as possible. For a similar anti-windup synthesis, we follow the method in [Turner et al., 2004] and con- structaplant-orderLMIbasedanti-windupcompensator. If the anti-windup compensator is chosen to be static, the main limitation of the synthesis technique, as discussed in Chapter 3, is the fact that the corresponding LMI constraints are not always feasible, whereas it is proven in [Grimm et al., 2003] that a plant-order linear compensator is always feasible for large enough L 2 gain if the plant is asymptotically stable. Hence, we assume that the plant with unknown parameters is Hurwitz. Giventhesystemmatrices (A p ,B p ,C p ,0) with C p =I,theplanttransferfunction matrixisdefinedas: G(s) = (sI−A p ) −1 B p (8.14) If G(s)=N(s)M −1 (s) is a full-order right coprime factorization of the transfer functionmatrixofthelinearplant,theanti-windupcompensatorcanbedescribedby K aw (s) = M(s)−I N(s) (8.15) Apossiblerealizationoftheanti-windupcompensatoris (A aw ,B aw ,C aw ,D aw ) with A aw =A p +B p F , B aw =B p , C aw = F I , D aw = 0 where F istobechosensuchthat A p +B p F isaHurwitzmatrix. Itisdiscussedin[Turneretal.,2004]thatthefeasibilityofthefollowingsetofLMIs 144 X B p U−L T 0 L T U T B T p −L −2U I U 0 I −μI −I L U T −I −I < 0, (8.16) Q> 0, (8.17) U > 0, (8.18) μ> 0, (8.19) with X definedas: X =QA T p +A p Q+L T B T p +B p L (8.20) and U asadiagonalmatrix,canbeusedtoobtainasuitablechoiceof F givenby F =LQ −1 (8.21) with L and Q definedthroughtheLMIsolution. As A p is Hurwitz and the problem is a convex optimization problem, F is guar- anteed to exist. Stabilizability and detectability are the common assumptions for the design of an LQ controller and the synthesis of the LMI based plant-order anti-windup compensatorweemploy. In order to reflect the relative importance of performance and robustness, the designer might choose positive definite weighting matrices W p and W r as proposed in[Turneretal.,2004],whereasthedetailsofthediscussiononperformanceandrobust- nesstrade-offareomittedhere. 145 Onecanevaluate F bysolvingthecorrespondingsetofLMIslistedas: X B p U−L T 0 Q L T U T B T p −L −2U I 0 U 0 I −μI 0 −I Q T 0 0 −W −1 p 0 L U T −I 0 −W −1 r < 0, (8.22) Q> 0, (8.23) U > 0, (8.24) μ> 0, (8.25) with X asdefinedin(8.20). Theanti-windupcompensatorcanbeaugmentedusing ˙ x aw =A aw x aw +B aw (u act −sat(u act )) (8.26) y aw =C aw x aw +D aw (u act −sat(u act )) where y aw = h y aw 1 y T aw 2 i T , y aw 1 ∈R, y aw 2 ∈R n ,and n istheorderoftheplant. Given u act as the unsaturated control input, sat(u act ) represents the actual effective inputtotheplantincludingthesaturationnonlinearity. Theterm y aw 2 modifiesthecontroldesignin(8.12)as: ˙ x cm =A c x cm +B c (r−x)−B c y aw 2 (8.27) u =C c x cm +D c (r−x)−D c y aw 2 with A c , B c , C c and D c asdefinedbefore. Theterm y aw 1 furthermodifies u by −y aw 1 anddefines u act : 146 Figure8.1: Thenominalcontrolstructurewithanti-windupcompensator(Section8.3) u act =u−y aw 1 (8.28) where u is the modified control input in (8.27). The nominal control structure with anti-windupcompensatorissummarizedinFigure8.1. It is discussed in [Weston and Postlethwaite, 2000] that the augmented anti-windup compensator has no effect when the actuators are operating linearly, but it is used to modifythesystem’sbehaviorduringandfollowingasaturationeventtoensurestability and eventual escape from saturation, restoring the desired linear behavior reasonably quickly. It is also shown that the closed-loop system can be represented by a combina- tionoftheintendedlinearsystem,anonlinearloopincludingthesaturationfunctionand adisturbancefilter. It is assumed in [Weston and Postlethwaite, 2000] that the system behaves linearly mostofthetime,andthattherecoverybacktotheoriginallyintendedlinearbehavioris requiredifthecontrolleroutputcausessaturation. Thisunderlyingassumptionformost anti-windupsynthesistechniquesisalsosharedinourwork. Analternativemethodtoconstructtheplant-orderanti-windupcompensatorforsta- ble plants is proposed in [Sofrony et al., 2007] where the LMIs are eliminated from 147 the design process by developing a synthesis approach which relies on the solution to a singleARE. The adaptive law and the corresponding stability analysis we establish can be extendedandshowntobevalidforthisalternativemethodaswell. 8.4 An Adaptive LQ Control Design with Anti-windup Augmentation As the control design and the anti-windup augmentation are based on the explicit plant model,theplantparameterscanbeestimatedonlineandusedtoevaluatetheparameters oftheadaptiveLQcontrollerandtheanti-windupcompensator. All the terms in (8.7) are first filtered to avoid high frequency sensor noise amplifi- cation by the derivative term. We introduce a tuning parameter λ > 0 to define a first orderfilter 1/(s+λ). Foranyfixedsetofplantparametersweobtain s s+λ x | {z } def = z =A p 1 s+λ x+B p 1 s+λ sat(u act ) (8.29) ⇒z = h A p B p i | {z } def = θ ∗T 1 s+λ x 1 s+λ sat(u act ) | {z } def = φ (8.30) Using the n th order plant state vector x= h x 1 x 2 ... x n i T and the saturated 148 controlinput sat(u act ),theregressorvectorisdefinedas: φ = 1 s+λ h x 1 x 2 ... x n sat(u act ) i T (8.31) Theestimationmodelisgivenby ˆ z i =θ T i φ, i = 1,2,...,n (8.32) wheretheestimateof θ ∗ i attime t is θ i (t) = h ˆ a i1 ˆ a i2 ... ˆ a in ˆ b i i T (8.33) with ˆ a ij astheestimateof a ij and ˆ b i astheestimateof b i .Notethat a ij isthe ij th entryof A p ∈R n×n ,and b i , x i arethe i th entriesof B p ∈R n , x∈R n respectively. Weconstructthenormalizedestimationerror, i : i = z i −θ T i φ m 2 , i = 1,2,...,n (8.34) where m 2 isdefinedas: m 2 = 1+φ T φ (8.35) sothat φ m ∈L ∞ canbeguaranteed. In Chapter 3 for the adaptive law the discrete version of the LSA given in [Ioannou and Sun, 1996] is employed, whereas we can as well use the Gradient Algorithm for generating θ i (t) where ˙ θ i = Γ i i φ, Γ i = Γ T i > 0 (8.36) 149 Wecanthenguaranteeforanyconstant θ ∗ i that k ˙ θ i k =k ˙ ˜ θ i k≤kΓ i k i,2 k i mk kφk m (8.37) where ˜ θ i =θ i −θ ∗ i (8.38) is the vector of parameter errors for i = 1,2,...,n. Correspondingly, we define the matricesofparametererrorsas: ˜ A p = ˆ A p −A p (8.39) ˜ B p = ˆ B p −B p (8.40) Together with kφk m ∈L ∞ , and i m∈L 2 T L ∞ , it is implied by (8.37) that ˙ θ i ∈ L 2 T L ∞ .Asaresultwehave θ i ∈L ∞ and i , i m, ˙ θ i ∈L 2 T L ∞ for i = 1,2,...,n. To eliminate possible parameter drift due to modeling errors, estimates of param- eters are forced to lie inside compact sets by using parameter projection so that their boundedness is guaranteed. It is assumed that the lower and upper limits of each entry θ ∗ ij of θ ∗ i , j = 1,2,...,n+1, i = 1,2,...,n are known so that L ij ≤θ ∗ ij ≤U ij for j = 1,2,...,n+1, i = 1,...,n for some L ij , U ij ∈R. The projection sets are thus cross products of closed sets on the real line. The orthogonal term-by-term projection oneachsetisgivenasinChapter3by ¯ θ ij (t) =Pr{θ ij (t)} = θ ij (t), L ij ≤θ ij ≤U ij L ij , θ ij (k)<L ij U ij , θ ij (k)>U ij (8.41) where θ ij (t) is calculated online using the adaptive law given in (8.36) at time t for i = 1,2,...,n. 150 Figure8.2: Theadaptivecontrolstructurewithanti-windupcompensator(Section8.4) Wecandescribetheerrordynamicsforanyconstant r as: ˙ e =A p e+A p r+B p sat(u act ) (8.42) where the tracking error can be regulated if the LQ control structure with anti-windup modificationisimplemented. SinceA p andB p areunknown,K 1 andK 2 in(8.11)cannotbeevaluatedbysolving the ARE in (8.9), and the LQ controller in (8.12) cannot be constructed. Likewise, the set of LMIs cannot be solved to obtain F using (8.21), thus the anti-windup controller in(8.26,8.27)cannotbeimplemented. Instead,weproposeanadaptivecontroldesignmethodologybasedontheCertainty EquivalencePrincipleandimplement u act =− ˆ K 2 1 s (e+y aw 2 )− ˆ K 1 (e+y aw 2 )−y aw 1 (8.43) 151 wheretheanti-windupmodificationterms y aw 1 and y aw 2 aregeneratedas: y aw 1 = ˆ F(sI− ˆ A p − ˆ B p ˆ F) −1 ˆ B p [u act −sat(u act )] (8.44) y aw 2 = (sI− ˆ A p − ˆ B p ˆ F) −1 ˆ B p [u act −sat(u act )] (8.45) TheadaptiveLQcontrolgaincanbedenotedby ˆ K = h ˆ K 1 ˆ K 2 i (8.46) andevaluatedthroughthesolutionoftheAREin(8.9)usingtheparameterestimates ˆ A p and ˆ B p fromtheadaptivelaw. ˆ F canbefoundthroughtheLMIsolutionagainusingthe parameterestimates. Figure8.2demonstratestheblockdiagramfortheadaptivecontrol structureincludingtheanti-windupscheme. In order to reduce the computational complexity, replacing the continuous-time adaptive law with a hybrid one is proposed in [Ioannou and Sun, 1996] which we do notdiscusshereanyfurther. Theapplicationofthecontroldesignonanaircraftmodelwithinputconstraintsand unknown parameters is described in Chapter 3 where simulation results for different scenariosaredemonstrated. 8.5 StabilityAnalysis Thepurposeofanadaptivetrackingcontrolschemeingeneralistodesignthecontroller and parameter adjustment mechanism so that all signals in the closed-loop plant are bounded and the plant output tracks the reference signal as closely as possible. Any additional nonlinearities such as actuator saturation constraints introduce further chal- lengesinadaptivesystemsandneedtobeaddressed. 152 Foradaptivesystemssubjecttoinputsaturation,itisdiscussedin[Polycarpouetal., 2003] that one possible approach is to completely stop adaptation during saturation of thecontrolinputwhereitisalsomentionedthatthisad-hocmethodpreventsthetracking errorinducedbyactuatorconstraintsfromcorruptingparameterestimations,butthatthe stabilitypropertiesoftheclosed-loopsystemcannotbeestablished. Usingourindirectadaptivecontrolstructurewithanti-windupaugmentation,where the LQ control design assumptions coincide with the LMI based anti-windup feasibil- ity, we can establish the stability properties of the closed-loop system without stopping adaptation. Ifweaddandsubtracttheterms ˆ A p e and ˆ B p sat(u act ) tothetrackingerrordynamics in(8.42),weobtaintheequation ˙ e = ˆ A p e− ˆ B p sat(u act )+A p r− ˜ A p e− ˜ B p sat(u act ) (8.47) where − ˜ A p e and − ˜ B p sat(u act ) are the terms due to parameter errors. In order to express the saturated plant input sat(u act ), the state vector x, and thus the tracking error e intermsofthenormalizedestimationerror = h 1 2 ... n i T (8.48) we first write m 2 using the estimation model in (8.32) and the normalized estimation errorin(8.34,8.48)as: m 2 = − ˜ A p 1 s+λ x− ˜ B p 1 s+λ sat(u act ) (8.49) = − ˜ A p 1 s+λ (e+r)− ˜ B p 1 s+λ sat(u act ) Operatingoneachsideof(8.49)with (s+λ),weobtain 153 (s+λ)m 2 = − ˜ A p e− ˜ B p sat(u act )− ˙ ˜ A p 1 s+λ e− ˙ ˜ B p 1 s+λ sat(u act ) (8.50) − ˜ A p r− ˙ ˜ A p 1 s+λ r whichcanbecombinedwith(8.47)toderive ˙ e = ˆ A p e+(s+λ)m 2 + ˙ ˜ A p 1 s+λ e+ ˙ ˜ B p 1 s+λ sat(u act )+ ˜ A p r (8.51) + ˙ ˜ A p 1 s+λ r− ˆ B p sat(u act )+A p r Anewvectorvariableisnowdefinedas: ¯ e =e−m 2 (8.52) andissubstitutedinto(8.51)toobtain ˙ ¯ e = ˆ A p ¯ e+(λI + ˆ A p ) | {z } ∈L∞ m |{z} ∈L 2 m+ ˙ ˜ A p |{z} ∈L 2 1 s+λ e + ˙ ˜ B p |{z} ∈L 2 1 s+λ sat(u act ) | {z } ∈L∞ (8.53) + ˆ A p r |{z} ∈L∞ + ˙ ˜ A p |{z} ∈L 2 1 s+λ r | {z } ∈L∞ − ˆ B p sat(u act ) | {z } ∈L∞ To establish exponential stability of the homogeneous part of (8.53), all the roots of det(sI− ˆ A p (t)) should be in the negative half-plane. For first order plants where A p ∈R,thiscanbesatisfiedbysimplysettingtheupperprojectionlimitfortheestimate ˆ A p tobenegativeandinitializingtheparameterestimationproperly. Whileguaranteeing this property simply by projection is not trivial for higher order systems, the projection operatorcanbemodifiedtoallow ˆ A p tohaveeigenvalueswithnegativerealpartsonly aslongastheassumptioniskeptthattheplantisstable. The adaptive control input u act implemented in (8.43) does not affect the homoge- neouspartof(8.53)whenthesystemisinthesaturationphasesincetheactualeffective 154 control input to the plant is sat(u act ) which can be shown to be a constant, either positiveornegative. The parameter estimates, ˆ A p and ˆ B p are bounded by projection. The saturated control input, sat(u act ) and the reference signal, r are bounded by definition. The adaptivelawguaranteesthat m, ˙ ˜ A p , ˙ ˜ B p ∈ L ∞ . Nowwedefineafictitiousnormalizingsignal m f suchthat m f 2 =c+ke t k 2 2δ (8.54) where c> 0 issomeconstant. Using ¯ e definitionin(8.52)wecanwrite ke t k 2δ ≤ k¯ e t k 2δ +k(m 2 ) t k 2δ (8.55) andobtainaboundfor m 2 f usingthe L 2δ normproperties: m 2 f ≤ c+ck¯ e t k 2 2δ +ck(m 2 ) t k 2 2δ (8.56) Ontheotherhand,(8.53)impliesthat k¯ e t k 2 2δ ≤ ck(m 2 ) t k 2 2δ +ck( ˙ ˜ A p e) t k 2 2δ + ck( ˙ ˜ B p sat(u act )) t k 2 2δ (8.57) + ck( ˆ A p r) t k 2 2δ + ck( ˙ ˜ A p r) t k 2 2δ + ck( ˆ B p ) t k 2 2δ wherek( ˙ ˜ B p sat(u act )) t k 2 2δ andk( ˙ ˜ A p r) t k 2 2δ areknowntobeboundedbythepropertiesof theadaptivelaw,andtheboundednessof k( ˆ A p r) t k 2 2δ and k( ˆ B p ) t k 2 2δ canbeestablished byprojection. Using the definition of the fictitious normalizing signal m f , k( ˙ ˜ A p e) t k 2 2δ can be showntobeboundedby k( ˙ ˜ A p m f ) t k 2 2δ . 155 Wepartition φ in(8.31)as φ = h φ T 1 φ 2 i T where φ 1 = 1 s+λ h x 1 x 2 ... x n i T (8.58) φ 2 = 1 s+λ sat(u act ) (8.59) The state vector x and the error e bound φ 1 , and φ 2 is bounded by definition. Thus, m = p 1+φ T φ canbeshowntobeboundedby m f introducedin(8.54). If(8.57)issubstitutedin(8.56),theinequality m 2 f ≤ c+ck(mm f ) t k 2 2δ +ck( ˙ ˜ A p m f ) t k 2 2δ (8.60) ≤ c+c Z t 0 e −δ(t−τ) g 2 (τ)m 2 f (τ)dτ isestablishedwhere g 2 = 2 m 2 +k ˙ ˜ A p k 2 i,2 (8.61) Since m, ˙ ˜ A p ∈ L 2 , it follows that g ∈ L 2 . Bellman-Gronwall Lemma (see AppendixB),whichallowsustoobtainanexplicitboundfromanimplicitonegivenby anintegralinequality,canthenbeappliedtoshowthat m f ∈L ∞ .Boundednessof m f implies that m, ¯ e, x, φ∈L ∞ . Thus, we establish the boundedness of all signals in theclosed-loopevenwhenthesystemgoesintothesaturationphase. If we could show that ˆ A p and ˆ B p are in L 2 , as r and sat(u act ) are bounded by definition, this would imply that ¯ e ∈ L 2 . Together with m 2 ∈ L 2 , ¯ e ∈ L 2 would guarantee that e∈L 2 . Barbalat’s Lemma (see Appendix C) could then be used to establish that the tracking error e converges to zero as we also know that ˙ e∈L ∞ . However, ˆ A p and ˆ B p cannotbeshowntobein L 2 ,andalthoughsignalboundedness 156 canbeproved,trackingerror e cannotbeguaranteedtoconvergetozerowheneverthe adaptivecontrolleractsontheplantduringthesaturationphase. Eveninthecasewherethesystemparametersareknown,thenominalcontrollerwith anti-windup augmentation does not promise that the tracking error converges to zero in general while the system is in saturation. The anti-windup augmentation is employed to return the control input back to its linear region where the tracking response can be observed. After the linear performance recovery, the adaptive controller, too, guaran- tees that the tracking error converges to zero asymptotically, since in that case we have sat(u act )=u act =u, and the adaptive control scheme with anti-windup compensator reducestothestandardadaptiveLQcontrolstructurewherethesaturationlimitscanbe neglecteduntiltheactuatorshittheselimitsagain. Fortheexistenceofthesolution P(t) =P T (t)> 0 totheARE,thepair ( ˆ A, ˆ B) has tobestabilizableforall t≥ 0 asdiscussedbefore. Therelevantstabilizabilityproblem in indirect Adaptive Pole Placement (APP) control is addressed in [Ioannou and Sun, 1996]. Itisdiscussedthatifoneisdealingwithtime-varyingestimates,uniformitywith respect to time is guaranteed by requiring the level of stabilizability to be greater than some constant ∗ > 0. It is proposed that the level of stabilizability can be defined as theabsolutevalueofthedeterminantoftheSylvestermatrixoftheestimatedplantpoly- nomials. The same modification can also be added to our adaptive LQ control design withanti-windupaugmentation. Stabilizabilitycanbeassumedinthehigherordercase, where loss of stabilizability occurs at certain isolated manifolds in the parameter space when visited by the estimated plant parameters; therefore, one can easily argue that the lossofstabilizabilityisnotafrequentphenomenonthatoccursintheimplementationof APPcontrolschemes.[IoannouandSun,1996] It is stated in [Turner et al., 2004] that an anti-windup scheme is not expected to yield greater robustness margins than the linear system upon which it is constructed. 157 Similarly,iftheplantparametersareunknown,wecouldnotexpectabetterperformance fromtheadaptivecontroldesignwithanti-windupaugmentationincomparisonwiththe casewithknownplantparameters. Ontheotherhand,ourdesignproposesasystematic way to deal with plants with unknown parameters and inputs prone to saturation while guaranteeingclosed-loopstability. 8.6 Conclusion We consider the tracking performance of asymptotically stable systems with unknown plant parameters and input saturation constraints. Our analysis on the adaptive control design with anti-windup augmentation demonstrates that the perturbation terms due to parameter errors in the adaptive scheme do not destroy the stability properties of the system. Weestablishboundednessofallsignalsintheclosed-loopandmaintainstability undersaturation. Oncethecontrolinputreturnstoitslinearregion,theadaptivecontrol structure withthe anti-windupcompensatorreducestoanadaptivecontrollerwherethe inputconstraintscansimplybeignored,andthusthesystemresponsetrackstheexternal referencesignalasymptotically. 158 9. FURTHERRESULTS: INDIRECTADAPTIVECONTROLFOR SYSTEMSWITHINPUTRATESATURATION 9.1 Introduction Althoughrelativelylesscommonthanlimitationsoncontrolinputmagnitude,inputrate saturationisalsoobservedinmanypracticalcontrolsystems. Theeffectsofthemagnitudeandratelimitsoftheactuatorsontheperformanceofa combustion chamber pressurecontrol systemaresimulated in[Krsticet al., 1999]. The dominating effect of rate saturation in the control of jet engine compressors is recog- nizedin[Hardtetal.,2000]. Themagnitudeandratelimitsinpracticalimplementation of bleed valves as actuators for active stall control in aircraft engines are discussed in[Wangetal.,2002]. Thetimedelaysinanaircraftcontrolsystemcansharplydegradethehandlingquali- ties.[Hess,1984]Thedestabilizingeffectofthetimedelayassociatedwithratelimiting isconsideredin[Bergetal.,1996]. The actuator rate saturation in aircraft control is known to lead Pilot Induced Oscil- lations(PIOs)whichintherecentpasthavebeenrecognizedasthecauseoftheYF-22A 159 prototype advanced tactical fighter crash in April, 1992 [Dornheim, 1992] and subse- quent JAS 39 Gripen fighter crash at an air show in August, 1992 [Shifrin, 1993], dis- playing the significance of addressing actuator rate saturation in order to avoid relevant unstableoscillationsandtheriskofcrashlanding. In the quest for high performance, or in the event of a control surface failure, it is reasonable to expect additional complexities if the rate saturation is dominant. [Pachter et al., 1995,Pachter and Miller, 1998] Several control methodologies to handle input magnitudeandrateconstraintscanbefoundin[Lauvdaletal.,1997,Farrelletal.,2003, ChengandWang,2003,Tarbouriechetal.,2006]. There has been considerable attention focused on handling magnitude saturation through anti-windup compensation. The rate saturation nonlinearity, which remains inadequatelyaddressedincontroldesign,hasalsobeenconsideredasanactiveresearch topic by the anti-windup community in the past decade. [Teel and Buffington, 1997, Wu and Grigoriadis, 1999,Barbu et al., 2000,Wu and Soto, 2003,Barbu et al., 2005, Galeani et al., 2006,Brieger et al., 2007] The robustness characteristics of these design methodologiesareyettobeinvestigated. Thecontrolofinputmagnitudeconstrainedsystemswithunknownplantparameters is addressed using an LMI based anti-windup augmentation in Chapter 3, and the sta- bility analysis of the adaptive control design is provided in Chapter 8. In this chapter weproposeanadaptive controlstructurewithanti-windupcompensationtoremedythe rate saturation nonlinearities imposed on the control input, in particular if plants with unknownparametersareinvolved. Therate-limiterisintroducedasaninputnonlinearityinSection9.2wheretheactu- ator dynamics are absorbed partly into the plant and partly into the unconstrained con- troller chosen as an LQ design for state tracking. In the same section a Riccati based anti-windup compensator is constructed based on the augmented nominal system. In 160 Figure9.1: Arate-limiterwithfirst-orderlag(Section9.2) Figure9.2: Aninfinitelyfastactuatormodel(Section9.2) Section 9.3 the plant parameters are assumed to be unknown, and an indirect adaptive control design technique is developed using the estimates of the original linear plant parameters and the Certainty Equivalence Principle. Our conclusions are presented in Section9.4. 9.2 InputRateSaturation Arate-limitedactuatorcanbemodeledasafirst-orderlagandasymmetricrate-limiting nonlinearityasshownintheblockdiagraminFigure9.1. The time constant of the lag, τ determines the cut-off frequency of the equivalent linearactuator.[Briegeretal.,2007] Iftheactuatordynamicsareveryfast,τ isverysmall,andthesystemcanbemodeled as an ideal relay instead of a saturation element. A pure rate limiter is shown in Figure 9.2.[Bergetal.,1996] 161 Figure9.3: Therate-limitedsystem(Section9.2) The system matrices, (A p ,B p ,C p ,0) with C p = I form a state space realization withthestatevector x∈R n foralinearplantoforder n andtheinput u act ∈R.Note also that A p ∈R n×n , B p ∈R n , C p ∈R n×n , and the plant transfer function matrix is G(s)=(sI−A p ) −1 B p . Basedonthelinearplantmodel,theLQcontrollerwithparameters K = h K 1 K 2 i , K 1 ∈R 1×n , K 2 ∈R 1×n (9.1) canbeconstructedtomanipulatethestatesoftheplant. Usingthedefinitionofthesaturationnonlinearity, sat(u r ) canbewrittenas: sat(u r ) =sign(u r )min(|u r |,¯ u r ) (9.2) where ¯ u r is the rate saturation level posed at the actuator such that ¯ u r > 0, ¯ u r ∈R. In the presence of rate limits on the control actuators the system can be summarized as depictedinFigure9.3. We follow the absorption method in [Brieger et al., 2007] and subsume the time constantintothecontrollerandtheintegratorblockintothelinearplant. 162 Theaugmentedplant,whichistheoriginalplantmodifiedbyabsorbingtheintegra- torintheactuatordynamics,hasthetransferfunctionmatrix: ˜ G(s) = ˜ C p (sI− ˜ A p ) −1 ˜ B p (9.3) where ˜ A p = A p B p 0 0 , ˜ B p = 0 I , ˜ C p = I 0 0 I (9.4) Theaugmentedlinearplant’sstatevector, x a ∈R n+1 isdefinedas x a =[x T u act ] T , anditsinputis sat(u r )∈R.Thesystemmatricesare ˜ A p ∈R (n+1)×(n+1) , ˜ B p ∈R n+1 and ˜ C p ∈R (n+1)×(n+1) whichisanidentitymatrix. Theaugmentedcontrollawiswrittenas: u r =−(1/τ)K 1 e−(1/τ)K 2 Z t 0 e(τ)dτ−(1/τ)u act (9.5) where u act = 1 s sat(u r ).Apossiblerealizationoftheaugmentedcontrolis ˙ x c = ˜ A c x c + ˜ B c −e −u act (9.6) u r = ˜ C c x c + ˜ D c −e −u act with ˜ A c = 0, ˜ B c = h I 0 i , ˜ C c = (1/τ)K 2 , ˜ D c = h (1/τ)K 1 (1/τ) i (9.7) where ˜ A c ∈R n×n , ˜ B c ∈R n×(n+1) , ˜ C c ∈R 1×n and ˜ D c ∈R 1×(n+1) . 163 Once the system subject to input rate saturation is interpreted as an augmented sys- temwithinputmagnitudesaturationconstraintsasabove,therelevantanti-windupcom- pensatorcanbeconstructed. If ˜ G(s)= ˜ N(s) ˜ M −1 (s) is a full-order right coprime factorization of the transfer function matrix of the augmented plant, the anti-windup compensator can be described bythetransferfunctionmatrix: ˜ K aw (s) = ˜ M(s)−I ˜ N(s) (9.8) A possible realization of the anti-windup compensator is then given through the systemmatrices,( ˜ A aw , ˜ B aw , ˜ C aw , ˜ D aw ) : ˜ A aw = ˜ A p + ˜ B p F , ˜ B aw = ˜ B p , ˜ C aw = F I , ˜ D aw =0 where ˜ A aw ∈R (n+1)×(n+1) , ˜ B aw ∈R n+1 , ˜ C aw ∈R (n+2)×(n+1) , ˜ D aw ∈R n+2 , and F istobechosensuchthatthat ˜ A p + ˜ B p F isaHurwitzmatrix. Two equivalent methods to construct the anti-windup compensator guaranteeing similar L 2 performanceareLMIandRiccatibasedsynthesistechniques. Eitheroneof them can be chosen to define F , and the implementation of the compensator into the systemfollowssimilarlyafterwardsregardlessofthechoicemade. Instead of the LMI based anti-windup compensator design methodology followed in previous chapters, here we use the alternative method to construct the plant-order anti-windupcompensatorforstableplants. Accordingly, F canbechosenasdiscussedin[Briegeretal.,2007]using: F =− R+ W −1 ˜ B p T P (9.9) 164 Figure9.4: Thecontrolstructurewithanti-windupcompensatorforrate-limitedsystems (Section9.2) where P > 0 satisfiestheARE: ˜ A T p P +P ˜ A p −P ˜ B p R ˜ B T p P + ˜ C T p ˜ C p = 0 (9.10) and R, W , are design parameters that can be properly chosen as in [Brieger et al., 2007] for a similar L 2 performance as promised by the respective LMI based design methodology. Once F isdefinedusing(9.9),theanti-windupcompensatormatricescanbecalcu- latedandtheanti-windupcompensatorcanbeaugmentedas: ˙ x aw = ˜ A aw x aw + ˜ B aw (u r −sat(u r )) (9.11) y aw = ˜ C aw x aw + ˜ D aw (u r −sat(u r )) with y aw = h y aw 1 y T aw 2 i T , y aw 1 ∈R and y aw 2 ∈R n+1 . Theanti-windupmodificationterm, y aw 2 isincorporatedintotheaugmentedcontrol structurein(9.6): 165 ˙ x cm = ˜ A c x cm + ˜ B c −e −u act − ˜ B c y aw 2 (9.12) u m = ˜ C c x cm + ˜ D c −e −u act − ˜ D c y aw 2 with ˜ A c , ˜ B c , ˜ C c and ˜ D c asdefinedbefore. Theterm y aw 1 furthermodifiesthecontrolinput u m in(9.12)as: u r =u m −y aw 1 (9.13) The nominal control structure with anti-windup compensation and the augmented linearplantsubjecttoinputmagnitudesaturationaredisplayedinFigure9.4. 9.3 AnIndirectAdaptiveControlDesign As the control structure is based on the explicit plant model, online estimations can be used to design an adaptive control with anti-windup augmentation for the case with unknownplantparameters. To avoid high frequency sensor noise amplification by the derivative term, we use a firstorderfilter 1/(s+λ) where λ> 0. Foranyfixedsetofplantparameters,weobtain s s+λ x | {z } def = z =A p 1 s+λ x+B p 1 s+λ 1 s sat(u r ) (9.14) ⇒z = h A p B p i | {z } def = θ ∗T 1 s+λ x 1 s+λ 1 s sat(u r ) | {z } def = φ (9.15) 166 where sat(u r ) is the time derivative of u act , the control input which includes the rate saturationnonlinearityandisdirectlyappliedtotheoriginalplant. As long as the control input does not hit the rate saturation limits ¯ u r or −¯ u r , and remains in the region where sat(u r )=u r , the relation u act = 1 s sat(u r ) reduces to u act = 1 s u r . Since the anti-windup compensator remains inactive, the output of the augmentedcontroller, u m isdirectlyappliedtotheaugmentedplantas u r . Theregressorvectorisdefinedusingtheoriginalplantstatevector, x andthesatu- ratedcontrolinput, sat(u r ) : φ = 1 s+λ h x T 1 s sat(u r ) i T = 1 s+λ h x T u act i T = 1 s+λ x a (9.16) whichisthestatevectoroftheaugmentedlinearplantinitsfilteredform. Theestimationmodelisgivenby ˆ z i =θ T i φ, i = 1,2,...,n (9.17) Theestimateof θ ∗ i attime t isdefinedas: θ i (t) = h ˆ a i1 ˆ a i2 ... ˆ a in ˆ b i i T (9.18) with ˆ a ij as the estimate of a ij and ˆ b i as the estimate of b i , whereas a ij is the ij th entryof A p ∈R n×n ,and b i isthe i th entryof B p ∈R n . We can then construct the normalized estimation error, i and use the Gradient Algorithmfortheadaptivelawtogenerate θ i (t). Inordertoguaranteetheboundednessofparameterestimatesaprojectionalgorithm canalsobeincluded. 167 Figure 9.5: The adaptive control structure for rate-limited systems with unknown plant parameters(Section9.3) As the plant parameters are unknown, K 1 and K 2 in (9.5) cannot be evaluated by solving the ARE for the corresponding system. The modified controller in (9.12) thus cannotbeconstructed. Similarly, F cannotbeobtainedusing(9.9),andtheanti-windup controllerin(9.11)cannotbeimplemented. Based on the Certainty Equivalence Principle we implement an adaptive control design instead. The adaptive control gain ˆ K is evaluated through the solution of the ARE using the parameter estimates from the adaptive law. ˆ F is also found through the Riccati solution using the parameter estimates. Since we work on an artificial plant augmentation,wedonotneedtoestimateallelementsofthematrices ˜ A p , ˜ B p and ˜ C p butonlytheoriginalplantparameters, A p and B p inordertoconstructtheestimateof ˜ A p .Wedonotestimate ˜ B p and ˜ C p sincetheyareknown. The block diagram of the adaptive control structure including the anti-windup scheme is illustrated in Figure 9.5 where both the augmented control and the anti- windupcompensatorparametersareevaluatedonline. 168 9.4 Conclusion In this chapter we consider systems prone to input rate saturation and implement an extensionoftheindirectadaptivecontrolstructurewithanti-windupcompensation. For systems with unknown plant parameters, the oscillations in the system response caused by the control input rate saturation and the corresponding performance degradation can potentially be suppressed by employing the proposed indirect adaptive control design. Thestabilitypropertiesoftheschemeareunderinvestigation. 169 10. SUMMARYANDCONCLUDINGREMARKS Wediscussthechallengesandadvancesinstaticsoaringflightwhereverticallyrisingair currentsareemployedinordertoimprovetheaircraftperformance. Usingtheformula- tionofstaticsoaringscenariosweproposeoptimumtrajectorygenerationalgorithmsfor UAVs in autonomous soaring applications. For simulation purposes we use an aircraft model which is modified for powerless flight. The flight range, endurance and overall velocityenhancementsthroughtheuseofourautonomouspathplanningmethodologies areillustratedforgliders. Recent soaring studies assume the aircraft dynamics to be known and to remain unaffected under changing flight conditions. The physical limitations such as actuator saturation nonlinearities are ignoredin the control design which isobserved to result in seriousperformancedegradationinrelevantflighttests. Avoidingtheassumptionsthatthesystemistime-invariantandthemodelparameters are known, we use ideas from robust adaptive control and anti-windup design method- ologies to develop an adaptive control scheme based on LQ control with disturbance rejectionandanti-windupaugmentation. Thesaturation-typenonlinearitypresentinthe systemisaddressedbycombininganadaptivelawwiththecontroldesignusingtheCer- taintyEquivalencePrincipleforasymptoticallystableplantswithmagnitudelimitations onthecontrolinput. 170 We are motivated by the recent LMI based anti-windup compensator design strate- gies for linear time-invariant models as the corresponding feasibility, stability and per- formance guarantees appear to be promising for the design of their adaptive counter- parts. Our results encourage the development and implementation of similar adaptive anti-windupschemesingeneral. Simulations for various scenarios verify the performance of the adaptive control design with anti-windup augmentation. The results validate that the flight requirements canbesatisfiedevenunderconditionsofgreatuncertaintyandambiguity. Itisobserved that the proposed adaptive control design also accommodates possible control surface failures to some extent during flight despite the presence of significant saturation limits ontheactuators. Perfectsensoryinformationandaccurateremotethermaldataareassumedingeneral in static soaring applications where the atmospheric disturbances acting on the system duringtheinter-thermalflightarealsoignored. Inanattempttorelaxtheseassumptions anddescribeamorerealisticflightenvironment,wemodeltheverticalgustdisturbances and the sensor measurement noise asstochastic processes. An adaptiveLQG controller isimplementedwiththerelevantmodificationsinthesystem. Theperformanceachieved by exploiting the practical statistical properties of random noise is shown to be compa- rabletothecasewithnoprocessormeasurementnoiseeffects. We suggest alternative routing strategies for UAVs which can benefit from a geo- graphically dispersed set of thermals in the flight region by deviating from the original flight path. The optimal soaring problem is investigated in detail by employing a VRP statement. The features of this formulation define the glider UAV as the single vehicle inservice,andthethermallocationsascustomers,terminalsorintermodalfacilities. We further consider duration constraints such as limited thermal lifespan and formulate the problem in the form of a VRPTW where the aircraft aims to climb the thermals when 171 theyareavailableforclimbing. Wedetermineatimecostfunctionanddeveloparecur- siveexactalgorithminvolvingPreprocessing,RouteOptimizationandRouteValidation stepstodefinetheoptimaltime-feasibleflightpath. If the thermal data is not deterministic, an SVRP formulation can be used, whereas time-varying thermal strength has its interpretation in the form of a DVRP statement. Several practical issues relevant to static soaring applications can thus be conveniently formulated by the VRP extensions, and one can benefit from the available solution methodologiesfortheseformulations. We broaden the application of our results to scenarios with dense thermal regions wherethecomputationaltimewouldseverelyincreaseifanexactalgorithmisemployed todeterminetheoptimalflightroute. Thecontrollerdesignremainsvalidinsuchasce- nario,however,thepracticalityofthetrajectorygenerationalgorithmissuggestedtobe improved. Heuristic methods are used to provide simplicity in terms of computations possiblyatthecostofaccuracyoftheproposedsolution. Asanotherapproximatesolu- tion methodology GA-based SP algorithm is effectively applied for the navigation of UAVsacrossrelativelyhighnumberofthermals. The tools of stability analysis are used to investigate the control design with anti- windup augmentation in the adaptive context. Upon combining the control structure with an adaptive law, we establish the closed-loop stability. Consequently, a systematic waytodealwithplantswithunknownparametersisprovidedusingcontrolinputswhich arepronetosaturation. Aratelimiterisadynamicnonlinearityandanothercommoncontrolinputconstraint deserving of attention in the system performance discussion. We further consider the anti-windupsynthesismethodologiesforrate-saturatedsystems. Asanextensionofthe currentresultsweconstructanindirectadaptivecontroldesignforplantssubjecttoinput ratesaturation. 172 In this research we investigate practical flight applications regarded as sources of competitive advantages for autonomous UAVs. However, the routing approach also appearstobevaluabletoassistpilotsinmannedsoaringflightwheretheycanmakeuse ofanonboardcomputationalsystemforrouteselectioninsteadofpurehumandecision. The scope of this research can be expanded to the soaring performance of UAVs in powered flight. The corresponding trajectory optimization problem would require a quantitativecomparisonofsoaringflighttotherelevantexpenditureofanonboardpower supplyandchangethecharacteroftheoptimizationproblemincreasingthepopularityof therelevantguidancestrategies. Usingproperflighttechniquestheenergyconsumption canbereducedwithsignificantimpactonfueleconomy. Soaringthereforestandsoutas anenvironmentfriendlytechnology. Future work in this area possibly involves the flight performance optimization con- cept for UAV formations flying across thermal regions. Multiple aircraft can more effi- ciently detect thermal locations. However, there are additional cooperative planning requirementstobeidentified. AUAVformationforinstancedemandsaddressingsafety issues such as collision avoidance. Moreover, the dynamics of different gliders in the formation might be represented by polar curves which are not identical, and thus might needtobediscussedseparately. It might also be of interest to researchers to investigate the use of adaptive anti- windup compensator design mechanisms in the presence of output constraints such as sensorsaturationnonlinearities. 173 REFERENCES [Abreuetal.,1992] Abreu,V.J.,Barnes,J.E.,andHays,P.B.(1992). Observationsof windswithanincoherentlidardetector. AppliedOptics,31(22):4509–4514. [Agarwaletal.,1989] Agarwal, Y., Mathur, K., and Salkin, H. M. (1989). A set- partitioning-based exact algorithm for the vehicle routing problem. Networks, 19(7):731–749. [AhnandRamakrishna,2002] Ahn, C. W. and Ramakrishna, R. S. (2002). A genetic algorithm for shortest path routing problem and the sizing of populations. IEEE TransactionsonEvolutionaryComputation,6(6):566–579. [Allen,2005] Allen, M. (2005). Autonomous soaring for improved endurance of a smalluninhabitedairvehicle. In43rdAIAAAerospaceSciencesMeetingandExhibit, Reno,NV. [Allen,2006] Allen,M.(2006). Updraftmodelfordevelopmentofautonomoussoaring uninhabited air vehicles. In 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno,NV. [AltinkemerandGavish,1991] Altinkemer,K.andGavish,B.(1991). Parallelsavings basedheuristicsforthedeliveryproblem. OperationsResearch,39(3):456–469. [Alvarezetal.,2007] Alvarez, A., Garau, B., and Caiti, A. (2007). Combining net- works of drifting profiling floats and gliders for adaptive sampling of the ocean. In 2007 IEEE International Conference on Robotics and Automation, pages 157–162, Roma,Italy. [AndersonandMoore,1990] Anderson, B. D. O. and Moore, J. B. (1990). Optimal Control: LinearQuadraticMethods. PrenticeHall,EnglewoodCliffs,NJ. [BalinskiandQuandt,1964] Balinski, M. L. and Quandt, R. E. (1964). On an integer programforadeliveryproblem. OperationsResearch,12(2):300–304. 174 [Barbuetal.,2005] Barbu, C., Galeani, S., Teel, A. R., and Zaccarian, L. (2005). Nonlinear anti-windup for manual flight control. International Journal of Control, 78(14):1111–1129. [Barbuetal.,2000] Barbu, C., Reginatto, R., Teel, A. R., and Zaccarian, L. (2000). Anti-windup for exponentially unstable linear systems with inputs limited in mag- nitude and rate. In Proceedings of the 2000 American Control Conference, pages 1230–1234,Chicago,IL. [Bergetal.,1996] Berg, J. M., Hammett, K. D., Schwartz, C. A., and Banda, S. S. (1996). Ananalysisofthedestabilizingeffectofdaisychainedrate-limitedactuators. IEEETransactionsonControlSystemsTechnology,4(2):171–176. [BertsekasandGallager,1987] Bertsekas, D. and Gallager, R. (1987). Data Networks. PrenticeHall,EnglewoodCliffs,NJ. [Bertsekas,1999] Bertsekas,D.P.(1999). NonlinearProgramming. AthenaScientific, Belmont,MA. [Bertsimas,1992] Bertsimas, D. J. (1992). A vehicle routing problem with stochastic demand. OperationsResearch,40(3):574–585. [BlantonandWainwright,1993] Blanton,J.L.andWainwright,R.L.(1993). Multiple vehicle routing with time and capacity constraints using genetic algorithms. In Pro- ceedingsofthe5thInternationalConferenceonGeneticAlgorithms,pages452–459, Champaign,IL. [Boslough,2002] Boslough, M. (2002). Autonomous dynamic soaring platform for distributedmobilesensorarrays. Technicalreport,SandiaNationalLaboratories. [Boydetal.,1994] Boyd, S., Ghaoui, L. E., Feron, E., and Balakrishnan, V. (1994). Linear Matrix Inequalities in System and Control Theory. Society for Industrial and AppliedMathematics,Philadelphia,PA. [Briegeretal.,2007] Brieger, O., Kerr, M., Leißling, D., Postlethwaite, I., Sofrony, J., and Turner, M. C. (2007). Anti-windup compensation of rate saturation in an exper- imental aircraft. In Proceedings of the 2007 American Control Conference, pages 924–929,NewYorkCity. [Caoetal.,2002] Cao, Y.-Y., Lin, Z., and Ward, D. G. (2002). H ∞ antiwindup design for linear systems subject to input saturation. AIAA Journal of Guidance, Control, andDynamics,25(3):455–463. [ChadwickandGossard,1983] Chadwick, R. B. and Gossard, E. E. (1983). Radar remote sensing of the clear atmosphere - review and applications. Proceedings of theIEEE,71(6):738–753. 175 [ChenandGuo,1991] Chen, H.-F. and Guo, L. (1991). Identification and Stochastic AdaptiveControl. Birkh¨ auser,Boston,MA. [Chen,1981] Chen, M. K. (1981). From paleoaeronautics to altostratus - a technical history of soaring. In AIAA Aircraft Systems and Technology Conference, Dayton, OH. [ChengandWang,2003] Cheng, J.-W. J. and Wang, Y.-M. (2003). Adaptive one-step- ahead control for non-minimum phase systems with input magnitude and rate con- straints. InternationalJournalofControl,76(15):1508–1515. [Christofidesetal.,1979] Christofides,N.,Mingozzi,A.,andToth,P.(1979). Thevehi- cle routing problem. In Combinatorial Optimization, pages 315–338. John Wiley & Sons,Inc.,Chichester,UK. [Christofidesetal.,1981] Christofides, N., Mingozzi, A., and Toth, P. (1981). Exact algorithms for the vehicle routing problem, based on spanning tree and shortest path relaxations. MathematicalProgramming,20(1):255–282. [ClarkeandWright,1964] Clarke,G.andWright,J.W.(1964). Schedulingofvehicles fromacentraldepottoanumberofdeliverypoints. OperationsResearch,12(4):568– 581. [Cochrane,1999] Cochrane, J. H. (1999). MacCready theory with uncertain lift and limitedaltitude. TechnicalSoaring,23(3):88–96. [Cone,1964] Cone,C.D.(1964).Thedesignofsailplanesforoptimumthermalsoaring performance. Technicalreport,NASALangleyResearchCenter. [Cormenetal.,2001] Cormen, T. H., Leiserson, C. E., Rivest, R. L., and Stein, C. (2001). IntroductiontoAlgorithms. TheMITPress,Cambridge,MA. [Cullenetal.,1981] Cullen, F. H., Jarvis, J. J., and Ratliff, H. D. (1981). Set partition- ingbasedheuristicsforinteractiverouting. Networks,11(2):125–143. [DantzigandRamser,1959] Dantzig, G. B. and Ramser, J. H. (1959). The truck dis- patchingproblem. ManagementScience,6(1):80–91. [Davidor,1991] Davidor, Y. (1991). Genetic Algorithms and Robotics: A Heuristic Strategy for Optimization. World Scientific Series in Robotics and Automated Sys- tems.WorldScientific. [deMatteisandMenichetti,2000] de Matteis, G. and Menichetti, P. (2000). Handling qualitiesofatowedsailplane. In38thAIAAAerospaceSciencesMeetingandExhibit, Reno,NV. 176 [DerrandLittle,1970] Derr, V. E. and Little, C. G. (1970). A comparison of remote sensing of the clear atmosphere by optical, radio, and acoustic radar techniques. AppliedOptics,9(9):1976–1992. [DesrochersandVerhoog,1989] Desrochers, M.andVerhoog, T.W.(1989). Amatch- ingbasedsavingsalgorithmforthevehicleroutingproblem. Technicalreport,Ecole desHautesEtudesCommerciales,Montreal,Canada. [Desrosiersetal.,1995] Desrosiers, J., Dumas, Y., Solomon, M. M., and Soumis, F. (1995). Time constrained routing and scheduling. In Network Routing, volume 8 of Handbooks in Operations Research and Management Science. North Holland Press, Amsterdam,TheNetherlands. [Doratoetal.,1995] Dorato,P.,Abdallah,C.,andCerone,V.(1995). Linear-Quadratic Control: AnIntroduction. PrenticeHall,EnglewoodCliffs,NJ. [Dornheim,1992] Dornheim,M.A.(1992). Reportpinpointsfactorsleadingto YF-22 crash. AviationWeekandSpaceTechnology,137(19):53–54. [DumitrescuandBoland,2003] Dumitrescu, I. and Boland, N. (2003). Improved pre- processing, labeling and scaling algorithms for the weight-constrained shortest path problem. Networks,42(3):135–153. [DuncanandPasik-Duncan,1992] Duncan,T.E.andPasik-Duncan,B.,editors(1992). Stochastic Theory and Adaptive Control. Lecture Notes in Control and Information Sciences.Springer-Verlag,NewYork. [Eppstein,1998] Eppstein, D. (1998). Finding the k shortest paths. SIAM Journal on Computing,28(2):652–673. [Falkenauer,1998] Falkenauer,E.(1998). GeneticAlgorithmsandGroupingProblems. JohnWiley&Sons. [Farrelletal.,2003] Farrell, J., Polycarpou, M., and Sharma, M. (2003). Adaptive backstepping with magnitude, rate, and bandwidth constraints: Aircraft longitude control. InProceedingsofthe2003AmericanControlConference,pages3898–3904, Denver,CO. [FerrariandStengel,2004] Ferrari, S. and Stengel, R. F. (2004). Online adaptive critic flightcontrol. AIAAJournalofGuidance,Control,andDynamics,27(5):777–786. [Fuller,1995] Fuller,J.R.(1995).Evolutionofairplanegustloadsdesignrequirements. AIAAJournalofAircraft,32(2):235–246. 177 [Galeanietal.,2006] Galeani,S.,Onori,S.,Teel,A.R.,andZaccarian,L.(2006). Fur- ther results on static linear anti-windup design for control systems subject to magni- tude and rate saturation. In Proceedings of the 45th IEEE Conference on Decision andControl,pages6373–6378,SanDiego,CA. [GareyandJohnson,1979] Garey, M. R. and Johnson, D. S. (1979). Computers and Intractability: A Guide to the Theory of NP-completeness . Series of Books in the MathematicalSciences.W.H.Freeman,SanFrancisco,CA. [Geetal.,2001] Ge,S.S.,Hang,C.C.,Lee,T.H.,andZhang,T.(2001). StableAdap- tiveNeuralNetworkControl. KluwerAcademic,Boston,MA. [Geetal.,1998] Ge, S. S., Lee, T. H., and Harris, C. J. (1998). Adaptive Neural Net- workControlofRobotManipulators. WorldScientific,RiverEdge,NJ. [GillandZomaya,1998] Gill,M.A.andZomaya,A.Y.(1998). ObstacleAvoidancein Multi-Robot Systems: Experiments in Parallel Genetic Algorithms . World Scientific SeriesinRoboticsandIntelligentSystems.WorldScientific,RiverEdge,NJ. [Goldberg,2002] Goldberg,D.E.(2002). TheDesignofInnovation: Lessonsfromand forCompetentGeneticAlgorithms. KluwerAcademicPublishers,Norwell,MA. [Grefenstetteetal.,1985] Grefenstette, J. J., Gopal, R., Rosmaita, B. J., and Gucht, D.V.(1985). Geneticalgorithmsforthetravelingsalesmanproblem. InProceedings of the 1st International Conference on Genetic Algorithms, pages 160–168, Pitts- burgh,PA. [Grimmetal.,2003] Grimm, G., Hatfield, J., Postlethwaite, I., Teel, A. R., Turner, M. C., and Zaccarian, L. (2003). Antiwindup for stable linear systems with input saturation: An LMI-based synthesis. IEEE Transactions on Automatic Control, 48(9):1509–1525. [Grimmetal.,2001] Grimm,G.,Postlethwaite,I.,Teel,A.R.,Turner,M.C.,andZac- carian, L. (2001). Linear matrix inequalities for full and reduced order anti-windup synthesis.InProceedingsofthe2001AmericanControlConference,volume5,pages 4134–4139,Arlington,VA. [Grimmetal.,2002] Grimm, G., Teel, A. R., and Zaccarian, L. (2002). Robust LMI- basedlinearanti-windupdesign: Optimizingtheunconstrainedresponserecovery. In Proceedingsofthe41stIEEEConferenceonDecisionandControl,volume1,pages 293–298,LasVegas,NV. [GuptaandKrishnamurti,1997] Gupta, A. and Krishnamurti, R. (1997). Parallel algo- rithmsforvehicleroutingproblems. InProceedingsoftheFourthInternationalCon- ferenceonHighPerformanceComputing,pages144–151,Bangalore,India. 178 [GutmanandVelger,1988] Gutman,P.O.andVelger,M.(1988). Trackingtargetswith unknown process noise variance using adaptive kalman filtering. In Proceedings of the27thConferenceonDecisionandControl,volume1,pages869–874,Austin,TX. [Hadjiconstantinouetal.,1995] Hadjiconstantinou,E.,Christofides,N.,andMingozzi, A. (1995). A new exact algorithm for the vehicle routing problem based on q-paths andk-shortestpathsrelaxations. AnnalsofOperationsResearch,61(1):21–43. [Hardtetal.,2000] Hardt, M., Helton, J. W., and Kreutz-Delgado, K. (2000). Numeri- cal solution of nonlinear H 2 and H ∞ control problems with application to jet engine compressors. IEEETransactionsonControlSystemsTechnology,8(1):98–111. [Hendriks,1974] Hendriks, F. (1974). Dynamic soaring in shear flow. In 2nd Interna- tionalSymposiumontheTechnologyandScienceofLowSpeedandMotorlessFlight, Cambridge,MA. [Hernandezetal.,1995] Hernandez, F., LeTraon, P.-Y., and Barth, N. H. (1995). Opti- mizing a drifter cast strategy with a genetic algorithm. Journal of Atmospheric and OceanicTechnology,12(2):330–345. [Hess,1984] Hess, R. A. (1984). Analysis of aircraft attitude control systems prone to pilot-induced oscillations. AIAA Journal of Guidance, Control, and Dynamics, 7(1):106–112. [Hideetal.,2004] Hide, C., Moore, T., and Smith, M. (2004). Adaptive kalman fil- tering algorithms for integrating GPS and low cost INS. In Position Location and NavigationSymposium,pages227–233. [HillandBenton,1992] Hill,A.V.andBenton,W.C.(1992).Modelingintra-citytime- dependenttravelspeedsforvehicleschedulingproblems. JournaloftheOperational ResearchSociety,43(4):343–351. [Hoblit,1988] Hoblit, F. M. (1988). Gust Loads on Aircraft: Concepts and Applica- tions. AIAAEducationSeries. [Holland,1975] Holland, J. H. (1975). Adaptation in Natural and Artificial Systems. TheUniversityofMichiganPress,AnnArbor,MI. [Holland,1992] Holland, J. H. (1992). Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intel- ligence. TheMITPress,Cambridge,MA. [Houboltetal.,1964] Houbolt, J. C., Steiner, R., and Pratt, K. G. (1964). Dynamic response of airplanes to atmospheric turbulence including flight data on input and response. Technicalreport,NASALangleyResearchCenter. 179 [Hwang,2002] Hwang, H.-S. (2002). An improved model for vehicle routing problem withtimeconstraintbasedongeneticalgorithm. ComputersandIndustrialEngineer- ing,42(2-4):361–369. [Inagakietal.,1999] Inagaki, J., Haseyama, M., and Kitajima, H. (1999). A genetic algorithm for determining multiple routes and its applications. In Proceedings of the 1999 IEEE International Symposium on Circuits and Systems, volume 6, pages 137–140,Orlando,FL. [IoannouandKokotovic,1983] Ioannou, P. A. and Kokotovic, P. V. (1983). Adaptive Systems with Reduced Models. Lecture Notes in Control and Information Sciences. Springer-Verlag,NewYork. [IoannouandSun,1996] Ioannou, P. A. and Sun, J. (1996). Robust Adaptive Control. PrenticeHall,UpperSaddleRiver,NJ. [Itoetal.,1989] Ito,Y.,Kobori,Y.,Takehisa,M.,andMitsuta,Y.(1989). Development ofwindprofilingsodar. JournalofAtmosphericandOceanicTechnology,6(5):779– 784. [KahveciandIoannou,2007] Kahveci, N. E. and Ioannou, P. A. (2007). An indirect adaptive control design with anti-windup compensation: Stability analysis. In 46th IEEEConferenceonDecisionandControl,pages1294–1299,NewOrleans,LA. [KahveciandIoannou,2008] Kahveci, N. E. and Ioannou, P. A. (2008). Flight trajec- tory optimization for multiple autonomous glider unmanned aerial vehicles. In 46th AIAAAerospaceSciencesMeetingandExhibit,Reno,NV. [Kahvecietal.,2007a] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2007a). Adaptive control of unmanned aerial vehicles in atmospheric flight with reduced models. In2007SAEAeroTechCongressandExhibition,LosAngeles,CA. [Kahvecietal.,2007b] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2007b). A heuristic search algorithm for maneuvering of UAVs across dense thermal areas. In AIAA Guidance, Navigation and Control Conference and Exhibit, Hilton Head, SC. [Kahvecietal.,2007c] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2007c). Optimal static soaring of UAVs using vehicle routing with time windows. In 45th AIAAAerospaceSciencesMeetingandExhibit,Reno,NV. [Kahvecietal.,2007d] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2007d). Arobustadaptivecontroldesignforgliderssubjecttoactuatorsaturationnonlineari- ties. InProceedingsofthe2007AmericanControlConference,pages492–497,New YorkCity. 180 [Kahvecietal.,2007e] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2007e). Astochasticapproachtooptimalsoaringproblemandrobustadaptive LQG control. InInfotech@Aerospace2007ConferenceandExhibit,RohnertPark,CA. [Kahvecietal.,2008] Kahveci, N. E., Ioannou, P. A., and Mirmirani, M. D. (2008). Adaptive LQ control with anti-windup augmentation to optimize UAV performance inautonomoussoaringapplications. IEEETransactionsonControlSystemsTechnol- ogy(toappear). [Kapooretal.,1997] Kapoor, N., Teel, A. R., and Daoutidis, P. (1997). An invariant subspace technique for anti-windupsynthesis. In Proceedingsofthe1997American ControlConference,volume5,pages3083–3087,Albuquerque,NM. [Khalil,1996] Khalil,H.K.(1996). Adaptiveoutputfeedbackcontrolofnonlinearsys- tems represented by input-output models. IEEE Transactions on Automatic Control, 41(2):177–188. [Khalil,2002] Khalil, H. K. (2002). Nonlinear Systems. Prentice Hall, Upper Saddle River,NJ. [Krsticetal.,1999] Krstic, M., Krupadanam, A., and Jacobson, C. (1999). Self-tuning controlofanonlinearmodelofcombustioninstabilities. IEEETransactionsonCon- trolSystemsTechnology,7(4):424–436. [Kyleetal.,2005] Kyle, J., Evans, K., and Costello, M. (2005). Atmospheric wind energy extraction by a small autonomous glider. In AIAA Atmospheric Flight MechanicsConferenceandExhibit,SanFrancisco,CA. [Laporteetal.,2000] Laporte, G., Gendreau, M., Potvin, J.-Y., and Semet, F. (2000). Classicalandmodernheuristicsforthevehicleroutingproblem. InternationalTrans- actionsinOperationalResearch,7(4-5):285–300. [Laporteetal.,1986] Laporte,G.,Nobert,Y.,andArpin,D.(1986). Anexactalgorithm for solving a capacitated location-routing problem. Annals of Operations Research, 6(9):293–310. [Lauvdaletal.,1997] Lauvdal, T., Murray, R. M., and Fossen, T. I. (1997). Stabiliza- tion of integrator chains in the presence of magnitude and rate saturations; a gain scheduling approach. In Proceedings of the 36th IEEE Conference on Decision and Control,volume4,pages4004–4005,SanDiego,CA. [Lazos,2005] Lazos,B.S.(2005). Biologicallyinspiredfixed-wingconfigurationstud- ies. AIAAJournalofAircraft,42(5):1089–1098. 181 [Leungetal.,1998] Leung, Y., Li, G., and Xu, Z.-B. (1998). A genetic algorithm for themultipledestinationroutingproblems. IEEETransactionsonEvolutionaryCom- putation,2(4):150–161. [Lewis,1986] Lewis,F.L.(1986). OptimalEstimation. JohnWiley&Sons,Inc.,New York. [LewisandSyrmos,1995] Lewis, F. L. and Syrmos, V. L. (1995). Optimal Control. JohnWiley&Sons,Inc. [LinandKernighan,1973] Lin,S.andKernighan,B.W.(1973). Aneffectiveheuristic algorithmforthetraveling-salesmanproblem. OperationsResearch,21(2):498–516. [Lissaman,2005] Lissaman,P.(2005). Windenergyextractionbybirdsandflightvehi- cles. In43rdAIAAAerospaceSciencesMeetingandExhibit,Reno,NV. [Luetal.,2005] Lu, B., Wu, F., and Kim, S. (2005). Linear parameter-varying anti- windup compensation for enhanced flight control performance. AIAA Journal of Guidance,Control,andDynamics,28(3):494–505. [MacCready,1998] MacCready,P.B.(1998). Regenerativebattery-augmentedsoaring. InSelf-LaunchingSailplaneSymposium ,Elmira,NY. [MalandrakiandDaskin,1992] Malandraki,C.andDaskin,M.S.(1992). Timedepen- dent vehicle routing problems: Formulations, properties and heuristic algorithms. TransportationScience,26(3):185–200. [Manetal.,1997] Man, K. F., Tang, K. S., Kwong, S., and Halang, W. A. (1997). Genetic Algorithms for Control and Signal Processing. Advances in Industrial Con- trol.Springer-Verlag,London. [Marinetal.,1999] Marin, J. A., Radtke, R., Innis, D., Barr, D. R., and Schultz, A. C. (1999). Usingageneticalgorithmtodeveloprulestoguideunmannedaerialvehicles. In1999IEEEInternationalConferenceonSystems,Man,andCybernetics,volume1, pages1055–1060,Tokyo,Japan. [Maughmer,2003] Maughmer,M.D.(2003). Designofwingletsforhigh-performance sailplanes. AIAAJournalofAircraft,40(6):1099–1106. [Mehra,1970] Mehra, R. K. (1970). On the identification of variances and adaptive kalmanfiltering. IEEETransactionsonAutomaticControl,15(2):175–184. [MetzgerandHedrick,1974] Metzger, D. E. and Hedrick, J. K. (1974). Optimal flight paths for soaring flight. In 2nd International Symposium on the Technology and ScienceofLowSpeedandMotorlessFlight,Cambridge,MA. 182 [Milne,2001] Milne, G. (2001). Glider aerotow training - feedback theory interpre- tation and PC-based simulation. In AIAA Modeling and Simulation Technologies ConferenceandExhibit,Montreal,Canada. [MoleandJameson,1976] Mole,R.H.andJameson,S.R.(1976). Asequentialroute- buildingalgorithmemployingageneralisedsavingscriterion. OperationalResearch Quarterly,27(2):503–511. [Mooij,1998] Mooij, E. (1998). Linear Quadratic Regulator Design for an Unpow- ered,WingedRe-entryVehicle . DelftUniversityPress. [Mulderetal.,2001] Mulder, E. F., Kothare, M. V., and Morari, M. (2001). Multivari- able anti-windup controller synthesis using linear matrix inequalities. Automatica, 37(9):1407–1416. [Munetomoetal.,1998] Munetomo, M., Takai, Y., and Sato, Y. (1998). A migration schemeforthegeneticadaptiveroutingalgorithm. In1998IEEEInternationalCon- ference on Systems, Man, and Cybernetics, volume 3, pages 2774–2779, San Diego, CA. [Nelson,1998] Nelson,R.C.(1998). FlightStabilityandAutomaticControl. McGraw- Hill,NY. [Nikolosetal.,2003] Nikolos, I. K., Valavanis, K. P., Tsourveloudis, N. C., and Kostaras, A. N. (2003). Evolutionary algorithm based offline/online path planner for UAV navigation. IEEETransactionsonSystems,Man,andCybernetics-PartB: Cybernetics,33(6):898–912. [Noonkester,1976] Noonkester,V.R.(1976). Theevolutionoftheclearairconvective layer revealed by surface-based remote sensors. Journal of Applied Meteorology, 15(6):594–606. [Pachteretal.,1995] Pachter, M., Chandler, P. R., and Mears, M. (1995). Control reconfiguration with actuator rate saturation. In Proceedings of the 1995 American ControlConference,pages3495–3499,Seattle,WA. [PachterandMiller,1998] Pachter, M. and Miller, R. B. (1998). Manual flight control withsaturatingactuators. IEEEControlSystemsMagazine,18(1):10–20. [PaitandKassab,1997] Pait,F.M.andKassab,F.(1997). Parallelalgorithmsforadap- tivecontrol: Robuststability.InControlUsingLogic-BasedSwitching ,LectureNotes inControlandInformationSciences.Springer. [PanandBasar,1998] Pan, Z. and Basar, T. (1998). Adaptive controller design for trackinganddisturbanceattenuationinparametricstrict-feedbacknonlinearsystems. IEEETransactionsonAutomaticControl,43(8):1066–1083. 183 [Parketal.,2004] Park, J.-S., Penner, M., and Prasanna, V. K. (2004). Optimizing graph algorithms for improved cache performance. IEEE Transactions on Parallel andDistributedSystems,15(9):769–782. [Parle,2004] Parle,J.(2004). Preliminarydynamicsoaringresearchusingaradiocon- trolglider. In42ndAIAAAerospaceSciencesMeetingandExhibit,Reno,NV. [Parrott,1970] Parrott,G.C.(1970). Aerodynamicsofglidingflightofablackvulture coragypsatratus. JournalofExperimentalBiology,53:363–374. [PiersonandChen,1979] Pierson, B. L. and Chen, I. (1979). Minimum landing- approachdistanceforasailplane. AIAAJournalofAircraft,16(4):287–288. [Piggott,2002] Piggott, D. (2002). Gliding: A Handbook On Soaring Flight. A&C Black,London. [Polycarpouetal.,2003] Polycarpou, M., Farrell, J., and Sharma, M. (2003). On-line approximation control of uncertain nonlinear systems: Issues with control input sat- uration. In Proceedings of the 2003 American Control Conference, volume 1, pages 543–548,Denver,CO. [Popov,1973] Popov, V. M. (1973). Hyperstability of Control Systems. Springer- Verlag,NewYork. [PreparataandShamos,1985] Preparata, F. P. and Shamos, M. I. (1985). Computa- tional Geometry: An Introduction. Texts and Monographs in Computer Science. Springer-Verlag. [QiandZhao,2005] Qi, Y. C. and Zhao, Y. J. (2005). Energy-efficient trajectories of unmannedaerialvehiclesflyingthroughthermals.JournalofAerospaceEngineering, 18(2):84–92. [Reichmann,1993] Reichmann,H.(1993). Cross-CountrySoaring . SoaringSocietyof America. [RenandKumar,1994] Ren, W. and Kumar, P. R. (1994). Stochastic adaptive pre- diction and model reference control. IEEE Transactions on Automatic Control, 39(10):2047–2060. [Renaudetal.,1996] Renaud, J., Boctor, F. F., and Laporte, G. (1996). An improved petal heuristic for the vehicle routing problem. The Journal of the Operational ResearchSociety,47(2):329–336. [RosenandHedenstrom,2001] Rosen, M. and Hedenstrom, A. (2001). Gliding flight in a jackdaw: A wind tunnel study. Journal of Experimental Biology, 204(6):1153– 1166. 184 [SaabandVanPutte,1999] Saab, Y. and VanPutte, M. (1999). Shortest path planning on topographical maps. IEEE Transactions on Systems, Man, and Cybernetics, Part A:SystemsandHumans,29(1):139–150. [Saberietal.,1999] Saberi, A., Stoorvogel, A. A., and Sannuti, P. (1999). Control of LinearSystemswithRegulationandInputConstraints. Springer,NewYork. [SachsanddaCosta,2003] Sachs, G. and da Costa, O. (2003). Optimization of dynamic soaring at ridges. In AIAA Atmospheric Flight Mechanics Conference and Exhibit,Austin,TX. [SafonovandTsao,1995] Safonov,M.G.andTsao,T.-C.(1995). Theunfalsifiedcon- trolconcept: Adirectpathfromexperimenttocontroller. InFeedbackControl,Non- linear Systems and Complexity, Lecture Notes in Control and Information Sciences. Springer-Verlag. [SavelsberghandSol,1995] Savelsbergh, M. W. P. and Sol, M. (1995). The general pickupanddeliveryproblem. TransportationScience,29(1):17–29. [Sedgewick,1988] Sedgewick,R.(1988). Algorithms. Addison-Wesley. [Sellmann,2003] Sellmann, M. (2003). Cost-based filtering for shorter path con- straints. In 9th International Conference on Principles and Practice of Constraint Programming,pages694–708,Kinsale,Ireland. [Shannonetal.,2002] Shannon, H. D., Young, G. S., Yates, M. A., Fuller, M. R., and Seegar, W. S. (2002). Measurements of thermal updraft intensity over complex ter- rain using American white pelicans and a simple boundary-layer forecast model. Boundary-LayerMeteorology ,104(2):167–199. [Shifrin,1993] Shifrin, C. A. (1993). Sweden seeks cause of Gripen crash. Aviation WeekandSpaceTechnology,139(7):78–79. [Short,2005] Short, S. (2005). Birth of American soaring flight: A new technology. AIAAJournal,43(1):17–28. [Sima,1996] Sima,V.(1996). AlgorithmsForLinear-QuadraticOptimization . Marcel Dekker,NY. [Sofronyetal.,2007] Sofrony, J., Turner, M. C., and Postlethwaite, I. (2007). Anti- windup synthesis using Riccati equations. International Journal of Control, 80(1):112–128. [Solomon,1987] Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research, 35(2):254–265. 185 [Stacketal.,2004] Stack,J.R.,Smith,C.M.,andHyland,J.C.(2004). Efficientreac- quisition path planning for multiple autonomous underwater vehicles. In Techno- Ocean2004,pages1564–1569,Kobe,Japan. [StevensandLewis,1992] Stevens,B.L.andLewis,F.L.(1992). AircraftControland Simulation. JohnWiley&Sons,Inc. [Taoetal.,2004] Tao,G.,Chen,S.,Tang,X.,andJoshi,S.M.(2004). AdaptiveControl ofSystemswithActuatorFailures. Springer. [TaoandKokotovic,1996] Tao, G. and Kokotovic, P. V. (1996). Adaptive Control of SystemswithActuatorandSensorNonlinearities. JohnWiley&Sons,NewYork. [TarbouriechandGarcia,1997] Tarbouriech, S. and Garcia, G. (1997). Control of Uncertain Systems with Bounded Inputs. Lecture Notes in Control and Information Sciences.Springer-Verlag,NewYork. [Tarbouriechetal.,2006] Tarbouriech, S., Prieur, C., and da Silva, J. M. G. (2006). Stability analysis and stabilization of systems presenting nested saturations. IEEE TransactionsonAutomaticControl,51(8):1364–1371. [TeelandBuffington,1997] Teel, A. T. and Buffington, J. M. (1997). Anti-windup for an F-16’s daisy chain control allocator. In AIAA Guidance, Navigation, and Control Conference,NewOrleans,LA. [TothandVigo,2002] Toth,P.andVigo,D.,editors(2002). TheVehicleRoutingProb- lem. SIAM Monographs on Discrete Mathematics and Applications. Society for IndustrialandAppliedMathematics,Philadelphia. [TsakalisandIoannou,1993] Tsakalis, K. S. and Ioannou, P. A. (1993). Linear Time VaryingPlants: ControlandAdaptation. PrenticeHall,UpperSaddleRiver,NJ. [Turneretal.,2004] Turner, M. C., Herrmann, G., and Postlethwaite, I. (2004). Accounting for uncertainty in anti-windup synthesis. In Proceedings of the 2004 AmericanControlConference,volume6,pages5292–5297,Boston,MA. [TurnerandPostlethwaite,2004] Turner, M. C. and Postlethwaite, I. (2004). A new perspective on static and low order anti-windup synthesis. International Journal of Control,77(1):27–44. [Vidyasagar,1993] Vidyasagar,M.(1993). NonlinearSystemsAnalysis. PrenticeHall, EnglewoodCliffs,NJ. [Villotaetal.,2006] Villota, E., Kerr, M., and Jayasuriya, S. (2006). A study of con- figurations for anti-windup control of uncertain systems. In Proceedings of the 45th IEEEConferenceonDecisionandControl,pages6193–6198,SanDiego,CA. 186 [Wangetal.,2006] Wang, R., Paul, A., Stefanovic, M., and Safonov, M. G. (2006). Cost detectability and stability of adaptive control systems. International Journal of RobustandNonlinearControl,17(5-6):549–561. [Wangetal.,2002] Wang,Y.,Yeung,S.,andMurray,R.M.(2002). Bifurcationcontrol of rotating stall with actuator magnitude and rate limits: Part II - control synthesis andcomparisonwithexperiments. Automatica,38(4):611–625. [WarkandHolt,1994] Wark, P. and Holt, J. (1994). A repeated matching heuris- tic for the vehicle routing problem. Journal of the Operational Research Society, 45(10):1156–1167. [WestonandPostlethwaite,2000] Weston, P. F. and Postlethwaite, I. (2000). Linear conditioning for systems containing saturating actuators. Automatica, 36(9):1347– 1354. [Wharington,1998] Wharington, J. (1998). Autonomous Control of Soaring Aircraft by Reinforcement Learning. PhD thesis, Royal Melbourne Institute of Technology, Melbourne,Australia. [WharingtonandHerszberg,1998] Wharington, J. and Herszberg, I. (1998). Control of a high endurance unmanned air vehicle. In 21st Congress of the International CounciloftheAeronauticalSciences,Melbourne,Australia. [Wharington,2004] Wharington, J. M. (2004). Heuristic control of dynamic soaring. In5thAsianControlConference,volume2,pages714–722,Melbourne,Australia. [WuandGrigoriadis,1999] Wu, F. and Grigoriadis, K. M. (1999). LPV-based control of systems with amplitude and rate actuator saturation constraints. In Proceedings of the 1999 American Control Conference, volume 5, pages 3191–3195, San Diego, CA. [WuandSoto,2003] Wu, F. and Soto, M. (2003). Extended LTI anti-windup control withactuatormagnitudeandratesaturations. InProceedingsofthe42ndIEEECon- ferenceonDecisionandControl,pages2786–2791,Maui,HI. [WuandJayasuriya,2006] Wu, W. and Jayasuriya, S. (2006). An internal model con- trol based anti-windup scheme for stable uncertain plants with input saturation. In Proceedings of the 45th IEEE Conference on Decision and Control, pages 5424– 5428,SanDiego,CA. [Xiongetal.,2005] Xiong, Z., Hao, Y., Wei, J., and Li, L. (2005). Fuzzy adaptive kalman filter for marine INS/GPS navigation. In Proceedings of the IEEE Inter- national Conference on Mechatronics and Automation, volume 2, pages 747–751, NiagaraFalls,Canada. 187 [YangandZhao,2005] Yang, H. I. and Zhao, Y. J. (2005). A practical strategy of autonomous dynamic soaring for unmanned aerial vehicles in loiter missions. In AIAAInfotech@Aerospace,Arlington,VA. [Yazetal.,1999] Yaz, E., Gao, Y., and Olejniczak, K. (1999). Comparison of adaptive filterperformanceinestimatingthenoisestatisticsforharmonicmodels. InProceed- ings of the 38th Conference on Decision and Control, volume 5, pages 5070–5075, Phoenix,AZ. [ZaccarianandTeel,2002] Zaccarian,L.andTeel,A.R.(2002).Acommonframework for anti-windup, bumpless transfer and reliable designs. Automatica, 38(10):1735– 1744. [ZaccarianandTeel,2004] Zaccarian, L. and Teel, A. R. (2004). Nonlinear scheduled anti-windup design for linear systems. IEEE Transactions on Automatic Control, 49(11):2055–2061. [Zhao,2004] Zhao, Y. J. (2004). Optimal patterns of glider dynamic soaring. Optimal ControlApplicationsandMethods: WileyInterScience,25(2):67–89. [ZhaoandQi,2004] Zhao,Y.J.andQi,Y.C.(2004). Minimumfuelpowereddynamic soaringofunmannedaerialvehiclesutilizingwindgradients. OptimalControlAppli- cationsandMethods: WileyInterScience,25(5):211–233. 188 APPENDICES AppendixA TheApplicationofPontryagin’sMinimumPrinciple Given the horizontal distance, x and the altitude, h, [Metzger and Hedrick, 1974] discussesthesimplifiedequationsofmotionasfollows: ˙ x =V x , ˙ h =V z ⇒ dh dx = V z V x whereinourcasewehave h(t 0 ) def = 0 h(t b ) =h t 0 + x(t b )−x(t 0 ) V x 189 Ifthecostfunctionforeachglide-and-thermalintervalisdefinedasthetotaltime, t total = Z x(t b ) x(t 0 ) 1 V x dx+ −h(t b ) W , theHamiltonianfunction, H canbewrittenas: H = 1 V x +λ× V z V x , andtheEuler-Lagrangeequationsare dh dx = V z V x , h(t 0 ) = 0 dλ dx =− ∂H ∂h = 0, λ(x(t b )) = ∂ ∂h(t b ) −h(t b ) W = −1 W By using the Pontryagin’s Minimum Principle and the polar curve approximation, wecancalculatetheoptimalvalueofthehorizontalvelocity, V x : V opt =argmin Vx H(V x ,λ,V z ) =argmin Vx 1 V x − 1 W V z V x =argmin Vx 1 V x − 1 W aV 2 x +bV x +c V x which is in agreement with our speed-to-fly definition given in (2.20) where there is an extra multiplicative term, x(t b )−x(t 0 ) affecting the minimum value of t total but not thecorrespondingargument, V opt . 190 AppendixB Bellman-GronwallLemma Let λ(t), g(t), k(t) benonnegativepiecewisecontinuousfunctionsoftime t.If thefunction y(t) satisfiestheinequality y(t) ≤ λ(t)+g(t) Z t t 0 k(s)y(s)ds, ∀t≥t 0 ≥ 0, itcanbeshownthat y(t) ≤ λ(t)+g(t) Z t t 0 λ(s)k(s)exp Z t s k(τ)g(τ)dτ ds, ∀t≥t 0 ≥ 0 Inparticular,if λ(t) =λ isaconstantand g(t) = 1, then y(t) ≤ λexp Z t t 0 k(s)ds , ∀t≥t 0 ≥ 0 Forfurtherdetailsthereadermightreferto[Vidyasagar,1993]. 191 AppendixC Barbalat’sLemma If lim t→∞ R t 0 f(τ)dτ existsandisfinite,and f(t) isauniformlycontinuousfunction, then lim t→∞ f(t) = 0. [Popov,1973] A special case of this general result can be stated as another lemma [Ioannou and Sun,1996]asfollows: If f , ˙ f ∈ L ∞ and f ∈ L p forsome p ∈ [1,∞), then f(t)→ 0 as t→∞. 192
Abstract (if available)
Abstract
The objective of meeting higher endurance requirements remains a challenging task for any type and size of Unmanned Aerial Vehicles (UAVs). According to recent research studies significant energy savings can be realized through utilization of thermal currents. The navigation strategies followed across thermal regions, however, are based on rather intuitive assessments of remote pilots and lack any systematic path planning approaches. Various methods to enhance the autonomy of UAVs in soaring applications are investigated while seeking guarantees for flight performance improvements.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Adaptive control with aerospace applications
PDF
Multiple model adaptive control with mixing
PDF
Practical adaptive control for systems with flexible modes, disturbances, and time delays
PDF
LQ feedback formulation for H∞ output feedback control
PDF
Bumpless transfer and fading memory for adaptive switching control
PDF
Relaxing convergence assumptions for continuous adaptive control
PDF
Distributed adaptive control with application to heating, ventilation and air-conditioning systems
PDF
Synthesis and analysis of high-performance controllers for high-speed autonomous aerobatic flight
PDF
Exploitation of wide area motion imagery
Asset Metadata
Creator
Kahveci, Nazli Eylem
(author)
Core Title
Robust adaptive control for unmanned aerial vehicles
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
12/17/2007
Defense Date
11/01/2007
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
actuator saturation,constrained control,flight trajectory optimization,OAI-PMH Harvest,robust adaptive control,unmanned aerial vehicles,vehicle routing problem
Language
English
Advisor
Ioannou, Petros A. (
committee chair
), Jonckheere, Edmond A. (
committee member
), Redekopp, Larry G. (
committee member
), Safonov, Michael G. (
committee member
)
Creator Email
kahveci@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m958
Unique identifier
UC1293746
Identifier
etd-Kahveci-20071217 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-487151 (legacy record id),usctheses-m958 (legacy record id)
Legacy Identifier
etd-Kahveci-20071217.pdf
Dmrecord
487151
Document Type
Dissertation
Rights
Kahveci, Nazli Eylem
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
actuator saturation
constrained control
flight trajectory optimization
robust adaptive control
unmanned aerial vehicles
vehicle routing problem