<%BANNER%>

hp-Pseudospectral Method for Solving Continuous-Time Nonlinear Optimal Control Problems

Permanent Link: http://ufdc.ufl.edu/UFE0042778/00001

Material Information

Title: hp-Pseudospectral Method for Solving Continuous-Time Nonlinear Optimal Control Problems
Physical Description: 1 online resource (173 p.)
Language: english
Creator: DARBY,CHRISTOPHER LIU
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2011

Subjects

Subjects / Keywords: CONTROL -- MESH -- METHODS -- NUMERICAL -- OPTIMAL -- OPTIMIZATION -- ORBIT -- PSEUDOSPECTRAL -- TRANSFER
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
Genre: Mechanical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: In this dissertation, a direct hp-pseudospectral method for approximating the solution to nonlinear optimal control problems is proposed. The hp-pseudospectral method utilizes a variable number of approximating intervals and variable-degree polynomial approximations of the state within each interval. Using the hp-discretization, the continuous-time optimal control problem is transcribed to a finite-dimensional nonlinear programming problem (NLP). The differential-algebraic constraints of the optimal control problem are enforced at a finite set of collocation points, where the collocation points are either the Legendre-Gauss or Legendre-Gauss-Radau quadrature points. These sets of points are chosen because they correspond to high-accuracy Gaussian quadrature rules for approximating the integral of a function. Moreover, Runge phenomenon for high-degree Lagrange polynomial approximations to the state is avoided by using these points. The key features of the hp-method include computational sparsity associated with low-order polynomial approximations and rapid convergence rates associated with higher-degree polynomials approximations. Consequently, the hp-method is both highly accurate and computationally efficient. Two hp-adaptive algorithms are developed that demonstrate the utility of the hp-approach. The algorithms are shown to accurately approximate the solution to general continuous-time optimal control problems in a computationally efficient manner without a priori knowledge of the solution structure. The hp-algorithms are compared empirically against local (h) and global (p) collocation methods over a wide range of problems and are found to be more efficient and more accurate. The hp-pseudospectral approach developed in this research not only provides a high-accuracy approximation to the state and control of an optimal control problem, but also provides high-accuracy approximations to the costate of the optimal control problem. The costate is approximated by mapping the Karush-Kuhn-Tucker (KKT) multipliers of the NLP to the Lagrange multipliers of the continuous-time first-order necessary conditions. It is found that if the costate is continuous, the hp-pseudospectral method is a discrete representation of the continuous-time first-order necessary conditions of the optimal control problem. If the costate is discontinuous, however, and a mesh point is at the location of a discontinuity, the hp-method is an inexact discrete representation of the continuous-time first-order necessary conditions. The computational efficiency and accuracy of the proposed hp-method is demonstrated on several examples ranging from problems whose solutions are smooth to problems whose solutions are not smooth. Furthermore, a particular application of a multiple-pass aeroassisted orbital maneuver is studied. In this application, the minimum-fuel transfer of a small maneuverable spacecraft between two low-Earth orbits (LEO) with an inclination change is analyzed. In the aeroassisted maneuvers, the vehicle deorbits into the atmosphere and uses lift and drag to change its inclination. It is found that the aeroassisted maneuvers are more fuel efficient than all-propulsive maneuvers over a wide range of inclination changes. The examples studied in this dissertation demonstrate the effectiveness of the hp-method.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by CHRISTOPHER LIU DARBY.
Thesis: Thesis (Ph.D.)--University of Florida, 2011.
Local: Adviser: Rao, Anil.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2011
System ID: UFE0042778:00001

Permanent Link: http://ufdc.ufl.edu/UFE0042778/00001

Material Information

Title: hp-Pseudospectral Method for Solving Continuous-Time Nonlinear Optimal Control Problems
Physical Description: 1 online resource (173 p.)
Language: english
Creator: DARBY,CHRISTOPHER LIU
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2011

Subjects

Subjects / Keywords: CONTROL -- MESH -- METHODS -- NUMERICAL -- OPTIMAL -- OPTIMIZATION -- ORBIT -- PSEUDOSPECTRAL -- TRANSFER
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
Genre: Mechanical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: In this dissertation, a direct hp-pseudospectral method for approximating the solution to nonlinear optimal control problems is proposed. The hp-pseudospectral method utilizes a variable number of approximating intervals and variable-degree polynomial approximations of the state within each interval. Using the hp-discretization, the continuous-time optimal control problem is transcribed to a finite-dimensional nonlinear programming problem (NLP). The differential-algebraic constraints of the optimal control problem are enforced at a finite set of collocation points, where the collocation points are either the Legendre-Gauss or Legendre-Gauss-Radau quadrature points. These sets of points are chosen because they correspond to high-accuracy Gaussian quadrature rules for approximating the integral of a function. Moreover, Runge phenomenon for high-degree Lagrange polynomial approximations to the state is avoided by using these points. The key features of the hp-method include computational sparsity associated with low-order polynomial approximations and rapid convergence rates associated with higher-degree polynomials approximations. Consequently, the hp-method is both highly accurate and computationally efficient. Two hp-adaptive algorithms are developed that demonstrate the utility of the hp-approach. The algorithms are shown to accurately approximate the solution to general continuous-time optimal control problems in a computationally efficient manner without a priori knowledge of the solution structure. The hp-algorithms are compared empirically against local (h) and global (p) collocation methods over a wide range of problems and are found to be more efficient and more accurate. The hp-pseudospectral approach developed in this research not only provides a high-accuracy approximation to the state and control of an optimal control problem, but also provides high-accuracy approximations to the costate of the optimal control problem. The costate is approximated by mapping the Karush-Kuhn-Tucker (KKT) multipliers of the NLP to the Lagrange multipliers of the continuous-time first-order necessary conditions. It is found that if the costate is continuous, the hp-pseudospectral method is a discrete representation of the continuous-time first-order necessary conditions of the optimal control problem. If the costate is discontinuous, however, and a mesh point is at the location of a discontinuity, the hp-method is an inexact discrete representation of the continuous-time first-order necessary conditions. The computational efficiency and accuracy of the proposed hp-method is demonstrated on several examples ranging from problems whose solutions are smooth to problems whose solutions are not smooth. Furthermore, a particular application of a multiple-pass aeroassisted orbital maneuver is studied. In this application, the minimum-fuel transfer of a small maneuverable spacecraft between two low-Earth orbits (LEO) with an inclination change is analyzed. In the aeroassisted maneuvers, the vehicle deorbits into the atmosphere and uses lift and drag to change its inclination. It is found that the aeroassisted maneuvers are more fuel efficient than all-propulsive maneuvers over a wide range of inclination changes. The examples studied in this dissertation demonstrate the effectiveness of the hp-method.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by CHRISTOPHER LIU DARBY.
Thesis: Thesis (Ph.D.)--University of Florida, 2011.
Local: Adviser: Rao, Anil.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2011
System ID: UFE0042778:00001


This item has the following downloads:


Full Text

PAGE 1

HP –PSEUDOSPECTRALMETHODFORSOLVINGCONTINUOUS-TIMENONLINEAR OPTIMALCONTROLPROBLEMS By CHRISTOPHERL.DARBY ADISSERTATIONPRESENTEDTOTHEGRADUATESCHOOL OFTHEUNIVERSITYOFFLORIDAINPARTIALFULFILLMENT OFTHEREQUIREMENTSFORTHEDEGREEOF DOCTOROFPHILOSOPHY UNIVERSITYOFFLORIDA 2011

PAGE 2

c r 2011ChristopherL.Darby 2

PAGE 3

ACKNOWLEDGMENTS Iwouldliketothankmyadvisor,Dr.AnilV.Raoforguidingmet hroughthejourney ofobtainingaPhD.Hehastaughtmetopursueexcellenceinall thatIdo.Ithankhimfor hisguidance.Ithankhimforneverholdingmetothestandard ,butinstead,holdingme toahigherstandard. Iwouldalsoliketothankthemembersofmycommittee:Dr.Will iamHager, Dr.WarrenDixon,andDr.NormanFitz-Coy.Dr.Hager'sexper ienceandwisdom hasgivenmenewperspectivesonhowtocontinuetogrowasare searcherandmy conversationswithDr.DixonandDr.Fitz-Coyhavealwaysgi venmemotivationto continuemywork. IwouldliketoalsothankmyfellowmembersofVDOL.Iwillneve rforgetthetimes wehavesharedwhileIhavebeenamemberofthislab.Whileeach memberofour laboratoryhasenrichedmylife,Iwouldspecicallyliketo thankDivyaGargforher everwillingnesstolookintoandanswerorponderquestions Imayhaveregardingour research.Nomatterhoweasyorfarfetchedmyquestionsmayb e,sheisalwayswilling toponderthemwithme. Finally,Iwouldliketothankmyfriends,familyandChristi eforbeingthereto supportandencouragemenotonlyduringmypursuitforaPhD,b ut,inalwaysbeing thereformeforallofmypursuitsinlife.Iknowthatwhateve rIdoandwhereverIgo, theywillalwaysbemyconstant. 3

PAGE 4

TABLEOFCONTENTS page ACKNOWLEDGMENTS ..................................3 LISTOFTABLES ......................................7 LISTOFFIGURES .....................................8 ABSTRACT .........................................10 CHAPTER 1INTRODUCTION ...................................12 2MATHEMATICALBACKGROUND .........................21 2.1OptimalControl .................................21 2.2Finite-DimensionalOptimization .......................25 2.2.1UnconstrainedMinimizationofaFunction ..............25 2.2.2EqualityConstrainedMinimizationofaFunction ...........28 2.2.3InequalityConstrainedMinimizationofaFunction ..........30 2.2.4InequalityandEqualityConstrainedMinimizationofa Function ..31 2.2.5Finite-DimensionalOptimizationAsAppliedtoOptimal Control ..32 2.3NumericalQuadrature .............................33 2.3.1Low-OrderNumericalIntegrators ...................34 2.3.2GaussianQuadrature .........................35 2.3.3Multiple-IntervalGaussianQuadrature ................41 2.3.4Summary ................................42 2.4LagrangePolynomialApproximation .....................44 2.4.1Single-IntervalLagrangePolynomialApproximations ........44 2.4.2Multiple-IntervalLagrangePolynomialApproximatio ns .......48 2.4.3Summary ................................49 3MOTIVATIONFORAN HP –PSEUDOSPECTRALMETHOD ...........52 3.1MotivatingExample1:ProblemwithaSmoothSolution ..........52 3.2MotivatingExample2:ProblemwithaNonsmoothSolution ........53 3.3DiscussionofResults .............................56 4 HP –PSEUDOSPECTRALMETHOD ........................58 4.1Multiple-IntervalContinuous-TimeOptimalControlPro blemFormulation .58 4.2 hp -Legendre-GaussTranscription .......................61 4.2.1The hp –LGPseudospectralTransformedAdjointSystem ......65 4.2.2Sparsityofthe hp –LGPseudospectralMethod ............70 4.3 hp –Legendre-Gauss-RadauTranscription ..................73 4.3.1The hp -LGRPseudospectralTransformedAdjointSystem .....78 4

PAGE 5

4.3.2Sparsityofthe hp –LGRPseudospectralMethod ...........82 4.4Examples ....................................86 4.4.1Example1-ProblemwithaSmoothSolution ............87 4.4.2Example2-Bryson-DenhamProblem ................90 4.4.2.1Case1:LocationofDiscontinuityUnknown ........93 4.4.2.2Case2:LocationofDiscontinuityKnown .........93 4.5Summary ....................................97 5 HP –MESHREFINEMENTALGORITHMS .....................100 5.1ApproximationofErrorsinaMeshInterval ..................101 5.2 hp –MeshRenementAlgorithm1 .......................103 5.2.1IncreasingDegreeofApproximationorSubdividing .........103 5.2.2LocationofIntervalsorIncreaseofCollocationPoin ts .......105 5.2.3StoppingCriteria ............................106 5.2.4 hp –MeshRenementAlgorithm1 ...................106 5.3 hp –MeshRenementAlgorithm2 .......................107 5.3.1IncreasingDegreeofApproximationorSubdividing .........107 5.3.2IncreaseinPolynomialDegreeinInterval k .............108 5.3.3DeterminationofNumberandPlacementofNewMeshPoint s ...108 5.4 hp –MeshRenementAlgorithm2 .......................109 6EXAMPLES ......................................112 6.1Example1-MotivatingExample1 ......................113 6.2Example2-ReusableLaunchVehicleRe-entryTrajectory .........116 6.3Example3:Hyper-SensitiveProblem .....................120 6.3.1Solutionsfor t f =50 ..........................120 6.3.2Solutionsfor t f =1000 .........................121 6.3.3Solutionsfor t f =5000 .........................123 6.3.4SummaryofExample6.3 .......................124 6.4Example4:MotivatingExample2 ......................128 6.5Example5:MinimumTime-to-ClimbofaSupersonicAircraft .......130 6.6DiscussionofResults .............................134 7AEROASSISTEDORBITALINCLINATIONCHANGEMANEUVER .......138 7.1EquationsofMotion ..............................139 7.1.1ConstraintsOntheMotionoftheVehicle ...............140 7.1.1.1Initialandnalconditions ..................141 7.1.1.2Interiorpointconstraints ...................141 7.1.1.3Vehiclepathconstraints ...................142 7.1.2TrajectoryEventSequence ......................142 7.1.3PerformanceIndex ...........................143 7.2Re-FormulationofProblem ..........................145 7.3OptimalControlProblem ............................146 7.4All-PropulsiveInclinationChangeManeuvers ................147 5

PAGE 6

7.5Results .....................................147 7.5.1UnconstrainedHeatingRate ......................149 7.5.2ConstrainedHeatingRate .......................153 7.5.3OptimalTrajectoryStructure ......................153 7.6Discussion ...................................159 7.7Conclusions ...................................161 8CONCLUSIONS ...................................162 8.1DissertationSummary .............................162 8.2FutureWork ...................................164 8.2.1ContinuedWorkonMeshRenement ................164 8.2.2DiscontinuousCostate .........................165 REFERENCES .......................................166 BIOGRAPHICALSKETCH ................................173 6

PAGE 7

LISTOFTABLES Table page 3-1MatrixDensityandCPUTimefor p –and h –LGRMethodsforExample2. ...57 6-1InitialGridusedinExample6.1 ...........................114 6-2SummaryofAccuracyandSpeedforExample6.1 ................115 6-3InitialGridusedinExample6.2 ...........................119 6-4SummaryofAccuracyandSpeedforExample6.2 ................119 6-5InitialGridusedinExample6.3.1 ..........................121 6-6SummaryofAccuracyandSpeedforExample6.3.1 ...............121 6-7InitialGridusedinExample6.3.2 ..........................123 6-8SummaryofAccuracyandSpeedforExample6.3.2 ...............124 6-9InitialGridusedinExample6.3.3 ..........................124 6-10SummaryofAccuracyandSpeedforExample6.3.3 ...............127 6-11InitialGridusedinExample6.4 ...........................130 6-12SummaryofAccuracyandSpeedforExample6.4 ................130 6-13InitialGridusedinExample6.5 ...........................133 6-14SummaryofAccuracyandSpeedforExample6.5 ................134 7-1PhysicalandAerodynamicConstants. .......................140 7

PAGE 8

LISTOFFIGURES Figure page 2-1AnExtremalCurve x ( t ) andaComparisonCurve x ( t ) .............23 2-2 e vs. andaThreeIntervalTrapezoidRuleApproximation. ...........34 2-3Errorvs.NumberofIntervalsforTrapezoidRuleApproxima tionofEq.(2-66). .35 2-4Legendre-GaussPointsDistribution. ........................39 2-5Legendre-Gauss-RadauPointsDistribution. ....................40 2-6Errorvs.NumberofLGPointsforApproximationofEq.(2-66) .........41 2-7Errorvs.NumberofLGorLGRPointsforApproximationofEq.( 2-66) .....43 2-8Approximationof f ( )=1 = (1+25 2 ) UsingUniformlySpacedSupportPoints. 46 2-9Approximationof f ( )=1 = (1+25 2 ) UsingLGSupportPoints. ........47 2-10Errorvs.NumberofLinearApproximatingIntervalsofEq.( 2-86) ........49 2-11Errorvs.NumberofSupportPointsforApproximatingEq.(286) ........50 2-12Errorvs.NumberofLGorLGRIntervalsforApproximationo fEq.(2-86) ....51 3-1Errorfor p and h –LGRApproximationtoEqs.(3-1)–(3-3). ............54 3-2Errorvs.CollocationPtsfor p –and h –LGRApproximationtoEqs.(3-5)–(3-8) .56 4-1CollocationandDiscretizationPointsbyLGMethods. ..............65 4-2QualitativeRepresentationofconstraintJacobianfor LGpointswithLarge N k 74 4-3QualitativeRepresentationofconstraintJacobianfor LGpointswithSmall N k .75 4-4CollocationandDiscretizationPointsbyLGRMethods. .............79 4-5QualitativeRepresentationofconstraintJacobianfor LGRpointswithLarge N k 85 4-6QualitativeRepresentationofconstraintJacobianfor LGRpointswithSmall N k 86 4-7SolutiontoExample4.4.1 ..............................89 4-8ConvergenceRatesforLGMethodsforExample4.4.1. .............90 4-9ConvergenceRatesforLGRMethodsforExample4.4.1. ............91 4-10Case1:ConvergenceforLGmethodsforExample4.4.2. ............94 4-11Case1:ConvergenceforLGRmethodsforExample4.4.2. ...........95 8

PAGE 9

4-12Case1:SolutionbyaUniform h – 2 -Method. ....................96 4-13Case2:Convergencefor p and h – 3 –LGandLGRMethodsforExample4.4.2 .98 5-1LocationofIntervalsifSubdivisionisChosen. ...................105 5-2LocationofIntervalsbasedonCurvature. .....................110 6-1Statevs.TimeontheFinalMeshforExample6.2 ................118 6-2StateandControlvs.TimeontheFinalMeshforExample6.3. 1 ........122 6-3StateandControlvs.TimeontheFinalMeshforExample6.3. 3. ........125 6-4Statevs.TimeontheFinalMeshforExample6.3.3. ...............126 6-5Controlvs.TimeontheFinalMeshforExample6.4. ...............131 6-6Statevs.TimeontheFinalMeshforExample6.5for hp (2) – 3 and =10 4 ...133 7-1SchematicofTrajectoryDesignforAeroassistedOrbitalT ransferProblem. ...144 7-2SchematicofAll-PropulsiveSingleImpulseandBi-Parabolic Transfers. .....148 7-3 V min ,vs. i f ,and m ( t f ) = m 0 vs. n .........................150 7-4 V 1 V 2 V 3 ,and V D Impulsesvs. n ......................151 7-5Changein i byAtmosphere, i aero ,andtransfertime, T ,vs. n for Q = 1 ...152 7-6 V min ,vs. n forVariousMaxHeatingRatesandFinalInclinations. .......154 7-7 m ( t f ) = m 0 ,vs. n forVarious Q max ,and i f =(20,40,60,80) deg. .........155 7-8 S max ,vs. n ,forVarious Q max and i f =(20,40,60,80) deg. ............156 7-9 h vs. v Q vs. t ,and h vs. r for i f =40 -deg. ....................157 7-10 and vs. t for n =1 and i f =40 deg. ......................158 7-11 h vs. v for n =(1,2,3,4) and i f =40 degand Q max =400 (W/cm 2 ). .......159 7-12 Q ,vs. t for n =(1,2,3,4) i f =40 degand Q max =400 (W/cm 2 ). ........160 9

PAGE 10

AbstractofDissertationPresentedtotheGraduateSchool oftheUniversityofFloridainPartialFulllmentofthe RequirementsfortheDegreeofDoctorofPhilosophy HP –PSEUDOSPECTRALMETHODFORSOLVINGCONTINUOUS-TIMENONLINEAR OPTIMALCONTROLPROBLEMS By ChristopherL.Darby May2011 Chair:AnilV.RaoMajor:MechanicalEngineering Inthisdissertation,adirect hp –pseudospectralmethodforapproximatingthe solutiontononlinearoptimalcontrolproblemsisproposed .The hp –pseudospectral methodutilizesavariablenumberofapproximatinginterva lsandvariable-degree polynomialapproximationsofthestatewithineachinterva l.Usingthe hp –discretization, thecontinuous-timeoptimalcontrolproblemistranscribe dtoanite-dimensional nonlinearprogrammingproblem(NLP).Thedifferential-alg ebraicconstraintsofthe optimalcontrolproblemareenforcedatanitesetofcolloc ationpoints,wherethe collocationpointsareeithertheLegendre-GaussorLegend re-Gauss-Radauquadrature points.Thesesetsofpointsarechosenbecausetheycorresp ondtohigh-accuracy Gaussianquadraturerulesforapproximatingtheintegralo fafunction.Moreover, Rungephenomenonforhigh-degreeLagrangepolynomialappr oximationstothestateis avoidedbyusingthesepoints.Thekeyfeaturesofthe hp –methodincludecomputational sparsityassociatedwithlow-orderpolynomialapproximat ionsandrapidconvergence ratesassociatedwithhigher-degreepolynomialsapproxim ations.Consequently,the hp –methodisbothhighlyaccurateandcomputationallyefcie nt.Two hp –adaptive algorithmsaredevelopedthatdemonstratetheutilityofth e hp –approach.The algorithmsareshowntoaccuratelyapproximatethesolutio ntogeneralcontinuous-time optimalcontrolproblemsinacomputationallyefcientman nerwithoutaprioriknowledge ofthesolutionstructure.The hp –algorithmsarecomparedempiricallyagainstlocal( h ) 10

PAGE 11

andglobal( p )collocationmethodsoverawiderangeofproblemsandarefo undtobe moreefcientandmoreaccurate. The hp –pseudospectralapproachdevelopedinthisresearchnoton lyprovidesa high-accuracyapproximationtothestateandcontrolofano ptimalcontrolproblem, butalsoprovideshigh-accuracyapproximationstothecost ateoftheoptimalcontrol problem.ThecostateisapproximatedbymappingtheKarush-K uhn-Tucker(KKT) multipliersoftheNLPtotheLagrangemultipliersofthecon tinuous-timerst-order necessaryconditions.Itisfoundthatifthecostateiscont inuous,the hp –pseudospectral methodisadiscreterepresentationofthecontinuous-time rst-ordernecessary conditionsoftheoptimalcontrolproblem.Ifthecostateis discontinuous,however, andameshpointisatthelocationofadiscontinuity,the hp –methodisaninexact discreterepresentationofthecontinuous-timerst-orde rnecessaryconditions. Thecomputationalefciencyandaccuracyoftheproposed hp –methodis demonstratedonseveralexamplesrangingfromproblemswho sesolutionsaresmooth toproblemswhosesolutionsarenotsmooth.Furthermore,ap articularapplication ofamultiple-passaeroassistedorbitalmaneuverisstudie d.Inthisapplication,the minimum-fueltransferofasmallmaneuverablespacecraftb etweentwolow-Earth orbits(LEO)withaninclinationchangeisanalyzed.Intheae roassistedmaneuvers, thevehicledeorbitsintotheatmosphereandusesliftanddr agtochangeitsinclination. Itisfoundthattheaeroassistedmaneuversaremorefuelef cientthanall-propulsive maneuversoverawiderangeofinclinationchanges.Theexam plesstudiedinthis dissertationdemonstratetheeffectivenessofthe hp –method. 11

PAGE 12

CHAPTER1 INTRODUCTION Theobjectiveofanoptimalcontrolproblemistodeterminet hestateandcontrol thatoptimizeaperformanceindexsubjecttodynamicconstr aints,boundaryconditions, andpathconstraints.Theperformanceindexgenerallyrepr esentsadesiredmetricor combinationofmetrics(e.g.,fuelconsumption,time,ener gy).Solvingoptimalcontrol problemsfundamentallyhasitsrootsinthecalculusofvari ations.Variationalcalculus isusedtoformulateasetofnecessaryconditionsthatdene extremalsolutionsto anoptimalcontrolproblem[ 1 – 3 ].Usingthecalculusofvariations,analyticsolutions canbefoundtoanarrowrangeofsimpleoptimalcontrolprobl ems.Unfortunately,in practice,analyticsolutionstogeneraloptimalcontrolpr oblemsaredifcultorimpossible todetermine.Instead,suchproblemsmustbesolvednumeric ally. Numericalmethodsforoptimalcontrolfallintotwocategor ies:indirectanddirect methods[ 4 5 ].Inanindirectmethod,thenecessaryconditionsfromthec alculus ofvariationsarederivedandthesenecessaryconditionsre sultinaHamiltonian boundary-valueproblem(HBVP).Examplesofindirectmethodsin cludeshooting, multipleshooting[ 6 7 ],nitedifference[ 8 ],andcollocation[ 9 10 ].Indirectmethods arehighlyaccurateandoptimalityofthecontinuous-times ystemisreadilyveriable becausethesemethodscomputeextremalsolutions.Itisnot ed,however,thatindirect methodsareimpractical.Therstdisadvantageisthatthen ecessaryconditionsfrom thecalculusofvariationsmustbeformulatedwhichisdifc ultforproblemsofeven simpletomoderatecomplexity.Next,therst-orderoptima lityconditionsaremostoften extremelydifculttosolve.Furthermore,anindirectmeth odcanonlybeusedifgood initialguessestothestateandcostate(aquantityinheren tinthenecessaryconditions) areprovided.Inparticular,guessesonthecostatearevery difculttogeneratebecause thecostatehasnophysicalinterpretation.Finally,forpr oblemswhosesolutionshave 12

PAGE 13

activepathconstraints,aprioriknowledgeoftheswitchin gstructureofthesolutionmust beknown. Inadirectmethod,theoptimalcontrolproblemistranscrib edintoanite-dimensional nonlinearprogrammingproblem(NLP).Manywell-knownandro bustalgorithmsexist tosolveNLPs[ 11 – 13 ].Inadirectmethod,thenecessaryconditionsfromthecalc ulus ofvariationsarenotformulated.Moreover,directmethods donotrequireguesseson thecostate.Ingeneral,directmethodsarenotasaccuratea sindirectmethodsand ingeneraltherst-ordernecessaryconditionsfromthecal culusofvariationsarenot satised. Directshooting[ 14 – 16 ]isamethodwhichhasbeenextensivelyusedforlaunch vehicleandorbittransferapplications[ 4 ].Awell-knowncomputerimplementation ofdirectshootingisPOST[ 17 ].Successfuldirectshootingmethodsapproximate theoptimalcontrolproblemusingafewnumberofNLPvariabl es.Thesuccessof directshootingmethodsinoptimizingtrajectoriesforlau nchvehicleandorbittransfer applicationsislargelyduetothefactthatmanyofthesepro blemsaresimpleenough thattheirsolutionscanbeaccuratelyapproximatedusinga smallnumberofoptimization parameters.Asthenumberofvariablesusedinadirectshooti ngmethodgrows,the abilitytosuccessfullyuseadirectshootingmethodrapidl ydeclines.Furthermore, whenusingadirectshootingmethod,theswitchingstructur eofinactiveandactivepath constraintsmustbeknownapriori. Anotherclassofadirectmethodsaredirectcollocationmeth ods[ 5 14 18 – 23 ]. Inadirectcollocationmethod,thestateisapproximatedus ingasetoftrial(basis) functionsandasetofdifferential-algebraicconstraints areenforcedatanitenumber ofcollocationpoints.Incontrasttoindirectmethodsandd irectshooting,adirect collocationmethoddoesnotrequireanaprioriknowledgeof theactiveandinactivearcs forproblemswithinequalitypathconstraints.Furthermor e,directcollocationmethods aremuchlesssensitivetoinitialguessthantheaforementi onedindirectmethods 13

PAGE 14

anddirectshootingmethod.Someexamplesofcomputerimplem entationsofdirect collocationmethodsareSOCS[ 24 ],OTIS[ 18 ],DIRCOL[ 25 ],andGPOPS[ 26 ]. Inthepasttwodecades,directcollocationmethodshavebec omemorewidelyused inapproximatingthesolutionsofgeneralnonlinearoptima lcontrolproblems.Typically, directcollocationmethodsforoptimalcontrolaredevelop edas h (local)methods.Inan h –method,xedlow-degreeapproximationsareusedandthepr oblemisdividedinto manysubintervals.Mostoften,theclassofRunge-Kuttamet hods[ 4 19 21 27 – 30 ]are usedfor h –discretizationandconvergenceofthenumericaldiscreti zationisachievedby increasingthenumberofsubintervals[ 19 30 31 ]. h –methodsareutilizedextensively inapproximatingthesolutiontooptimalcontrolproblemsb ecausetheresultingNLP issparse.Foraspeciedproblemsize,sparsityintheNLPgr eatlyincreasesthe computationalefciency.Adrawbackto h –methods,however,isthatconvergenceisat apolynomialrate,andtherefore,anexcessivelylargenumb erofsubintervalsmaybe requiredtoaccuratelyapproximatethesolution.Forexamp le,inthedesignofanoptimal low-thrusttrajectorytothemoon,Ref.[ 32 ]requiredanapproximatingNLPof211031 variablesand146285constraints. Incontrastto h –methods,theclassofdirectcollocation p –methodsusesa smallxednumberofapproximatingintervals(oftenonlyas ingleintervalisused). Convergencetotheexactsolutionbya p –methodisachievedbyincreasingthe degreeofthepolynomialapproximationineachinterval.In optimalcontrol,single interval p –methodshavebeendevelopedsynonymouslyaspseudospectr almethods [ 22 26 33 – 52 ].Inapseudospectralmethod,thecollocationpointsarede nedas thepointsfromhighlyaccurateGaussianquadraturesrules .Furthermore,thebasis functionsaretypicallyChebyshevorLagrangepolynomials .Forproblemswhose solutionsaresmoothandwell-behaved,apseudospectralme thodconvergesatan exponentialrate[ 53 – 55 ]. 14

PAGE 15

Themostwelldevelopedpseudospectralmethodsarethe Gausspseudospectralmethod (GPM)[ 22 35 ],the Radaupseudospectralmethod (RPM)[ 38 – 40 ], andthe Lobattopseudospectralmethod (LPM)[ 33 ].Inthesemethods,thesetof collocationpointsaretheLegendre-Gauss(LG),LegendreGauss-Radau(LGR),and theLegendre-Gauss-Lobatto(LGL)points.Furthermore,th esethreesetsofpoints aredenedontheinterval 2 [ 1,1] anddifferprimarilyinhowtheendpointsare incorporated.TheLGpointsdonotincludeeitherendpoint, theLGRpointsincludeonly oneendpoint,andtheLGLpointsincludebothendpoints.Beca useoptimalcontrol problemsgenerallyhaveboundaryconditionsatboththeini tialandterminaltimes, theLGLpointsarethemostobviousofthethreesetsofGaussi anquadraturepoints touse.AlthoughLGLpointsarethemostobviouschoice,recen tresearchshows thattheymaynotbethemostappropriate.Ref.[ 42 ]showsthatthecostateestimate usingLGLpointstendstobenoisy.Inparticular,itwasexpl ainedinRef.[ 40 ]that theLGLdiscretecostateapproximationisrank-decientan dthenoisyLGLcostate approximationisduetotheoscillationsabouttheexactsol ution,andtheseoscillations havethesamebehaviorasthenullspaceoftheLGLapproximat ion.Furthermore,it wasshowninRef.[ 40 ]thatnotonlyistheLGLapproximationtothecostateinaccu rate, buttheapproximationtothecontrolmayalsobeinaccurate. Finally,itisnotedthat whenusingtheLGorLGRmethods[ 22 40 ],theresultingdiscreteapproximationsare high-accuracyapproximationstothesolutionoftheoptima lcontrolproblem. While p –methodsconvergeatanexponentialrateforproblemswhose solutionsare smoothandwell-behaved, p –methodshaveseverallimitations.First,ifthesolutiono f aproblemisnotsmooth,exponentialconvergencedoesnotho ldandtheconvergence ratestotheexactsolutionaresignicantlylower.Further more,ifthesolutiontoa problemissmoothbuthasveryhigh-orderbehaviors,anaccu ratesolutionisobtained onlyusinganunreasonablyhigh-degreepolynomialapproxi mation.Next,theNLP arisingfroma p –methodisdense.Becausethedensityofthe p –methodNLPis 15

PAGE 16

high,asthedegreeofpolynomialapproximationbecomeslar ge,theNLPbecomes computationallyintractable. Athirdclassofdirectcollocationmethodsthathasnotrece ivedasignicantamount ofattentioninapproximatingthesolutionstooptimalcont rolproblemsistheclassof hp –methods.An hp –methodisoneinwhichthenumberofintervals,placementof intervals,anddegreeofpolynomialapproximationineachi ntervalisallowedtovary. Previously, hp –methodshavebeenusedinmechanicsanduiddynamics.Inpa rticular, Refs.[ 56 – 60 ]describethemathematicalpropertiesof h p ,and hp –methodsfornite elements.Itshouldbenotedthatintheperiodfrom1975to19 95,morethan37,000 papersrelatedto h –methodsinniteelementshavebeenpublished[ 57 ].Inthatsame timeperiod,lessthan100publishedpapersrelatedtoeithe r p or hp –methodsinnite elementsexist[ 57 ].Therefore,historically,thevastmajorityofresearchh asbeenin h –methods.Somecommerciallyavailable h –niteelementsoftwareprogramsare MSC/NASTRAN,Cosmos/M,Abaqus,Aska,Adina,andAnsys[ 57 ].Withrespectto p – and hp –methodsinniteelements,someavailablesoftwareprogra msareMSC/PROBE, AppliedStructure,PHLEXandSTRIPE[ 57 ]. Inordertocombinethebenetsof h and p –methodsforapproximatingthesolution togeneraloptimalcontrolproblems,thatis,combinethesp arsityofan h –methodwith theexponentialconvergencerateofa p –method,inthisresearch,an hp –approach isdevelopedusingLGandLGRpoints.Furthermore,becauset hisresearchutilizes collocationatpointsdenedbyGaussianquadrature,wecal lthese hp –pseudospectral methods.Inadditiontodeveloping hp –pseudospectralmethodsthatprovidehigh accuracyapproximationstothestateandcontrol,wedevelo phighaccuracycostate estimationmethodsusingthe hp –LGandLGRmethods.Byprovidinganapproximation tothecostateandnotjustthestateandcontrol,vericatio nofoptimalityofthe originalcontinuous-timeproblemmaybedeterminedfromth edirect hp –method.The mappingoftheKarush-Kuhn-Tucker(KKT)multipliersoftheNLP tothecostateofthe 16

PAGE 17

continuous-timeoptimalcontrolproblemareshownbyrstd erivingtheKKTconditions oftheNLP,thenderivingthetransformedadjointsystem,an dthencomparingthe transformedadjointsystemtotherst-ordernecessarycon ditionsofthecontinuous-time optimalcontrolproblem.Becauseofthecontinuityconstrai ntonthestateatthemesh points,themappingfromtheKKTmultiplierstothecostatehav etwoforms.Ifthe costateiscontinuous,the hp –transformedadjointsystemisadiscreterepresentation ofthecontinuous-timerst-ordernecessaryconditionsof theoptimalcontrolproblem. If,however,thecostateisdiscontinuous,thetransformed adjointsystemisnolonger adiscreterepresentationofthecontinuous-timerst-ord ernecessaryconditions.It isfound,however,thathighlevelsofaccuracyarestillobt ainedinapproximatingthe solutiontothecontinuous-timeoptimalcontrolproblemif meshpointsareplacedatthe locationofthecostatediscontinuities. Todemonstratetheutilityofan hp –method,two hp –meshrenementalgorithms arepresented.Previousrenementalgorithmshavebeendeve lopedfor h –methods [ 30 61 ]whileanexampleofa p –meshrenementalgorithmisgiveninRef.[ 47 ]. Renementalgorithmsfor h –methodsaredesignedtoincreasethenumberoflow-order approximatingintervalswhereerrorsarelargestinapprox imatingthesolutionto thecontinuous-timeoptimalcontrolproblem.Byonlyincrea singthenumberof approximatingintervalswhereerrorsarelarge,theseren ementmethodsincrease thecomputationalefciencyovernaive h –methodapproximations(i.e.,auniformly distributedgridofapproximatingsubintervals).Next,in the p –renementalgorithmof Ref.[ 47 ],thealgorithmincreasedtheaccuracyoftheapproximatio n,butdidnotaddress computationalefciency.Large,densesingle-interval p –approximationsareusedon everyiterationofthealgorithmexceptforthenalone.Onl yonthenaliterationof thealgorithmistheapproximatingmeshdivided.Byusinglar ge,densesingle-interval p –approximations,theNLPissolvedinacomputationallyine fcientmanner. 17

PAGE 18

Inthe hp –algorithmsofthisresearch,thenumberofintervals,widt hofintervals,and polynomialdegreeineachintervalareallowedtovaryoneve ryiteration.Atwo-tiered strategyisemployedtodeterminethelocationsofinterval sanddegreeofpolynomial withineachintervaltoachieveaspeciedsolutionaccurac y.Ifthesolutionbeing approximatedbyanintervalissmooth,thenthedegreeofapp roximatingpolynomialis increased.Otherwise,ifthesolutionbeingapproximatedb yanintervalisnon-smooth, theintervalissubdividedintomorelow-ordersubinterval s.Thetwo hp –algorithmsare differentbecausetherstalgorithmismore p –biased,thatis,increasingthedegree ofanapproximatingpolynomialischosenmoreoften,andthe secondalgorithmis more h –biased,thatis,subdivisionofapproximatingsubinterva lsintomorelow-degree subintervalsismoreoftenchosen.Furthermore,the hp –algorithmsofthisresearchare signicantlydifferentfromtheaforementioned h and p –algorithms.The hp –algorithms ofthisresearcharedifferentfrom h –algorithmsbecausethedegreeofpolynomial isallowedtovaryineachapproximatinginterval.Finally, the hp –algorithmsofthis researcharedifferentfromthe p –algorithmofRef.[ 47 ]becausethenumberofintervals andlocationofintervalsareallowedtovaryoneachiterati on. Thecontributionsofthisdissertationareasfollows.Firs t, hp –directcollocation methodsaredescribedusingLegendre-GaussandLegendre-G auss-Radaupoints. Ifthecostateiscontinuous,thetransformedadjointsyste msofthe hp –LGandLGR methodsaredemonstratedtobediscreterepresentationsof thecontinuous-time rst-orderoptimalityconditionsofanoptimalcontrolpro blem.Asaresult,high-accuracy approximationsforthestate,controlandcostateofthesol utiontoanoptimalcontrol problemareobtained.Next,the hp –methodsofthisresearchcombinethecomputational sparsityof h –methodsandtheexponentialconvergenceratesof p –methodsonintervals whosesolutionissufcientlysmooth.Therefore,the hp –methodsresultinhighlyefcient andhighlyaccurateapproximations. 18

PAGE 19

Two hp –adaptivealgorithmsarepresented.Noaprioriknowledgeo fthesolution structureofanoptimalcontrolproblemisrequiredinorder forthesealgorithmsto efcientlyobtainhighlyaccuratesolutions.Thealgorith msiterativelyrenethe approximatingmeshusingatwo-tieredstrategytoeitherin creasethedegreeof polynomialapproximationinanintervalortosubdividethe intervalintomorelow-degree approximatingsubintervals.Therstalgorithmpresented isa p –biasedalgorithm, whilethesecondalgorithmpresentedisan h –biasedalgorithm.Thetwoalgorithms areusedtosolveseveralexamplesandtheutilityofeachalg orithmisdemonstrated. Furthermore,the hp –algorithmsofthisresearcharecomparedagainsteachothe rand with p and h –methods.Itisdemonstratedthattheexibilityofan hp –approachadds considerablerobustnesstoefcientlyapproximatingopti malcontrolproblemstoa highlevelofaccuracyregardlessofthestructureofthesol utionofanoptimalcontrol problem. Thisdissertationisdividedintothefollowingchapters.C hapter 2 denesthe mathematicalbackgroundnecessarytounderstandthe hp –methodspresented.A nonlinearoptimalcontrolproblemisstatedandthecontinu ous-timerst-ordernecessary conditionsforthisproblemarederived.Next,nite-dimen sionaloptimizationofa functionsubjecttoequalityandinequalityconstraintsis examinedandtheKKTsystem isderived.Finally,numericalintegrationandpolynomial interpolationareexamined. Chapter 3 providesmotivationforan hp –method.Itisdemonstratedthatthemethodthat bestapproximatesageneralnonlinearoptimalcontrolprob lemis problemdependent Therefore,strict h or p –methodsmaybecomputationallyexpensiveforaspecic problemtype.Chapter 4 denesthe hp –LGand hp –LGRdiscreteapproximationto thecontinuous-timeoptimalcontrolproblem.Forthe hp –LGand hp –LGRmethods, theKKTconditionsarederivedandarerelatedtotherst-orde rnecessaryconditions oftheoptimalcontrolproblem.Finally,thesparsitypatte rnoftheNLPdenedbythe hp –methodsofthisresearchareexamined.Itisshownthatthes parsityofthe hp –LG 19

PAGE 20

andLGRmethodsisdependentuponthedegreeofapproximatin gpolynomialineach interval.Chapter 5 describesthetwo hp –adaptivealgorithmsofthisdissertation.First,a descriptionisgivenonhowerrorsareassessedinanyapprox imatinginterval.Second, thetwo hp –algorithmsaredescribedindetail.Therstalgorithmre nesthemeshby examiningtherelativedistributionoferrorwithineachap proximatingintervalandthe secondalgorithmrenesthemeshbyanalyzingthecurvature ofthestate.Chapter 6 presentsseveralexamplesthataresolvedbythetwo hp –algorithmsofthisdissertation aswellas h and p –methods.Itisdemonstratedthatthe hp –algorithmsarerobustto problemtypeandresultincomputationallyefcient,higha ccuracysolutionstoawide rangeofproblems.InChapter 7 ,aparticularapplicationofamultiple-passaeroassisted orbitalmaneuverisstudied.Theoptimalminimum-fueltran sferofasmallmaneuverable spacecraftbetweentwolow-Earthorbitswithaninclination changeisanalyzed.Optimal aeroglidemaneuversarecomparedwithall-propulsivemane uversoverawiderange ofinclinationchange.Furthermore,avariableheatingrat econstraintisenforcedfor theaeroglidemaneuvers.Itisdemonstratedthatanoptimal aeroglidemaneuveris morefuelefcientthananall-propulsivemaneuver.Finall y,Chapter 8 summarizesthe contributionsofthisdissertationandsuggestsfutureres earchdirections. 20

PAGE 21

CHAPTER2 MATHEMATICALBACKGROUND Thischapterdescribesthemathematicalbackgroundnecess arytounderstandthe hp –methodpresentedinthisresearch.Theadvancementsofthi sresearcharefounded onoptimalcontroltheoryandgeneratinghighlyaccurate,c omputationallyefcientNLP approximationstoacontinuous-timeoptimalcontrolprobl em.Numericalapproximation, functionapproximation,andquadratureapproximationthe oryarealsofoundationsof thisresearchandarediscussedindetail. 2.1OptimalControl Theobjectiveofanoptimalcontrolproblemistodeterminet hestateandcontrol thatoptimizeaperformanceindexsubjecttodynamicconstr aints,boundaryconditions, andpathconstraints.Theoptimalcontrolproblemstudiedi nthisresearchisgivenin Bolzaformasfollows.Minimizethecostfunctional J =( x ( t 0 ), t 0 x ( t f ), t f )+ Z t f t 0 g ( x ( t ), u ( t )) dt (2–1) subjecttothedynamicconstraints x ( t ) d x dt = f ( x ( t ), u ( t )), (2–2) theinequalitypathconstraints C ( x ( t ), u ( t )) 0 (2–3) andtheboundaryconditions ( x ( t 0 ), t 0 x ( t f ), t f )= 0 (2–4) where x ( t ) 2 R n isthestate, u ( t ) 2 R m iscontrol, f : R n R m R n C : R n R m R s : R n R R n R R q and t istime.Furthermore, : R n R R n R R istheMayercostand g : R n R m R istheLagrangian.Todetermineanextremal solution,thecalculusofvariationsisused.Thecalculuso fvariationsisusedtogenerate 21

PAGE 22

asetofrst-ordernecessaryconditionsthatmustbesatis edforacandidateoptimal solution(alsoknownasanextremalsolution).Anextremalso lutionmaybeaminimum, maximum,orsaddle.Second-ordersufciencyconditionsmus tbecheckedtoverifythat anextremalsolutionisaminimum.Thecontinuous-timerst -ordernecessaryconditions fortheoptimalcontrolproblemdenedbyEqs.( 2–1 )–( 2–4 )arenowderived. First,theconstraintsareadjoinedtothecostfunctionalt oformanaugmentedcost functionalas J a =( x ( t 0 ), t 0 x ( t f ), t f ) T ( x ( t 0 ), t 0 x ( t f ), t f ) (2–5) + Z t f t 0 [ g ( x ( t ), u ( t )) T ( t )(_ x ( t ) f ( x ( t ), u ( t ))) r T ( t ) C ( x ( t ), u ( t ))] dt where ( t ) isthecostateand r ( t ) and aretheLagrangemultiplierscorresponding toEqs.( 2–3 ),and( 2–4 ),respectively.Therst-ordervariationwithrespecttoa llfree variablesisgivenas J a = @ @ x ( t 0 ) x 0 + @ @ t 0 t 0 + @ @ x ( t f ) x f + @ @ t f t f (2–6) T T @ @ x ( t 0 ) x 0 T @ @ t 0 t 0 T @ @ x ( t f ) x f T @ @ t f t f +(( g T (_ x f ) r T C ) j t = t f ) t f (( g T (_ x f ) r T C ) j t = t 0 ) t 0 + Z t f t 0 @ g @ x x + @ g @ u u T (_ x f )+ T @ f @ x x + T @ f @ u u T x r T C r T @ C @ x x r T @ C @ u u dt Next,integratingbypartsthetermcontaining x suchthatitcanbeexpressedinterms of x ( t 0 ), x ( t f ), and x gives Z t f t 0 T x dt = T ( t f ) x ( t f )+ T ( t 0 ) x ( t 0 )+ Z t f t 0 T x dt (2–7) 22

PAGE 23

Tworelationshipsrelating x f to x ( t f ) and x 0 to x ( t 0 ) aregivenas x 0 = x ( t 0 )+_ x ( t 0 ) t 0 (2–8) x f = x ( t f )+_ x ( t f ) t f (2–9) TherelationshipsgiveninEq.( 2–8 )and( 2–9 )areinterpretedbyexaminingFig. 2-1 InFig. 2-1 x ( t ) representsanextremalcurveand x ( t ) representsaneighboring comparisoncurvewhere t 0 and x 0 arexed. x ( t f ) isthedifferencebetween x ( t ) and x ( t ) evaluatedat t f x f isthedifferencebetween x ( t ) and x ( t ) evaluatedattheendof eachcurve,respectively.AscanbeseeninFig. 2-1 ,Eq.( 2–9 )representsarst-order approximationto x f .Thesameresultholdsforfreeinitialconditionson t 0 and x 0 ,giving therelationshipofEq.( 2–8 ).SubstitutingEqs.( 2–7 ),( 2–8 )and( 2–9 )intoEq.( 2–6 )and Figure2-1.AnExtremalCurve x ( t ) andaComparisonCurve x ( t ) 23

PAGE 24

groupingterms,thevariationof J a isgivenas J a = @ @ x ( t 0 ) T @ @ x ( t 0 ) + T ( t 0 ) x 0 (2–10) + @ @ x ( t f ) T @ @ x ( t f ) T ( t f ) x f T + @ @ t 0 T @ @ t 0 g ( t 0 ) T ( t 0 ) f ( t 0 )+ r T ( t 0 ) C ( t 0 ) t 0 + @ @ t f T @ @ t f + g ( t f )+ T ( t f ) f ( t f ) r T ( t f ) C ( t f ) t f + Z t f t 0 ( @ g @ x + T @ f @ x r T @ C @ x + ) x +( @ g @ u + T @ f @ u r T @ C @ u ) u T (_ x f ) r T C i dt Therst-orderoptimalityconditionsareformedbysetting thevariationof J a equalto zerowithrespecttoeachfreevariable.Inordertowritethe seconditionsinaconvenient form,rst,theaugmentedHamiltonianisdened.Deningth eaugmentedHamiltonian as H ( x ( t ), u ( t ), ( t ), r ( t ))= g ( x ( t ), u ( t ))+ T ( t ) f ( x ( t ), u ( t )) r T ( t ) C ( x ( t ), u ( t )), (2–11) therst-orderoptimalityconditionsarethenexpressedas x T ( t )= @ H @ (2–12) T ( t )= @ H @ x (2–13) 0 = @ H @ u (2–14) T ( t 0 )= @ @ x ( t 0 ) + T @ @ x ( t 0 ) (2–15) T ( t f )= @ @ x ( t f ) T @ @ x ( t f ) (2–16) H ( t 0 )= T @ @ t 0 + @ @ t 0 (2–17) H ( t f )= T @ @ t f @ @ t f (2–18) 24

PAGE 25

Furthermore,usingthecomplementaryslacknesscondition r takesthevalue r i ( t )=0, when C i ( x ( t ), u ( t )) < 0, i =1,..., s (2–19) r i ( t ) < 0, when C i ( x ( t ), u ( t ))=0, i =1,..., s (2–20) Thenegativityof r i when C i =0 isinterpretedsuchthatimprovingthecostmayonly comefromviolatingtheconstraint[ 2 ].Furthermore, r i ( t )=0 when C i < 0 statesthat thisconstraintisinactive,andtherefore,ignored.Theco ntinuous-timerst-order optimalityconditionsofEqs.( 2–12 )–( 2–20 )deneasetofnecessaryconditions thatmustbesatisedforanextremalsolutionofanoptimalc ontrolproblem.The second-ordersufciencyconditionscanbeinspectedtodet erminewhethertheextremal solutionisaminimizingsolution. 2.2Finite-DimensionalOptimization Inadirectcollocationapproximationofanoptimalcontrol problem,thecontinuous-time optimalcontrolproblemistranscribedintoanonlinearpro grammingproblem(NLP). Whilethecontinuous-timeoptimalcontrolproblemisaninn ite-dimensionalproblemof ndingcontinuous-timefunctionsthatminimizeacostfunc tionalsubjecttodifferential andalgebraicconstraints,anNLPisanite-dimensionalpr oblemofndingavectorof staticparametersthatminimizeacostfunctionsubjecttoa lgebraicconstraints.Inthis sectionthetopicsofunconstrainedminimization,equalit yconstrainedminimization,and inequalityconstrainedminimizationofafunctionaredisc ussed.Moreover,aniterative Newtonmethodisshownfordeterminingaminimumofafunctio n. 2.2.1UnconstrainedMinimizationofaFunction Considerthefollowingproblemofdeterminingtheminimumo fafunctionsubjectto multiplevariableswithoutanyconstraints.Minimizetheo bjectivefunction J ( x ), (2–21) 25

PAGE 26

where x T =( x 1 ,..., x n ) .For x tobealocallyminimizingpoint,theobjectivefunction mustbegreaterwhenevaluatedatanyneighboringpoint,i.e ., J ( x ) > J ( x ). (2–22) Inordertodevelopasetofsufcientconditionsdeningalo callyminimizingpoint x rst,athreetermTaylorseriesexpansionaboutsomepoint, x ,isusedtoapproximate theobjectivefunctionas J ( x )= J ( x )+ g T ( x )( x x )+ 1 2 ( x x ) T H ( x )( x x )+ higherorderterms (2–23) wherethegradientvector, g ( x ) ,is g ( x ) r x J @ J @ x T = 266666664 @ J @ x 1 @ J @ x 2 ... @ J @ x n 377777775 (2–24) andthesymmetric n n Hessianmatrixis H ( x ) r xx J @ 2 J @ x 2 = 266666664 @ 2 J @ x 2 1 @ 2 J @ x 1 @ x 2 ... @ 2 J @ x 1 @ x n @ 2 J @ x 2 @ x 1 @ 2 J @ x 2 2 ... @ 2 J @ x 2 @ x n ... @ 2 J @ x n @ x 1 @ 2 J @ x n @ x 2 ... @ 2 J @ x 2 n 377777775 (2–25) For x tobealocalminimizingpoint,twoconditionsmustbesatis ed.First,anecessary conditionisthat g ( x ) mustbezero,i.e., g ( x )= 0 (2–26) Thenecessaryconditionbyitselfonlydenesanextremalpo intwhichcanbealocal minimum,localmaximum,orsaddlepoint.Inordertoensure x isalocalminimum,then 26

PAGE 27

anadditionalconditionthatmustbesatisedis ( x x ) T H ( x )( x x ) > 0. (2–27) Eqs.( 2–26 )and( 2–27 )togetherdenesufcientconditionsforalocalminimum. Inordertoiterativelyndaminimumoftheobjectivefuncti on,rst,consideragain theapproximationofEq.( 2–23 )usingaTaylorseriesexpansionaboutsomepoint x ( i ) suchthatanapproximationatanewpoint x ( i +1) isgivenas J ( x ( i +1) ) J ( x ( i ) )+ g ( i ) T ( x ( i +1) x ( i ) )+ 1 2 ( x ( i +1) x ( i ) ) T H ( i ) ( x ( i +1) x ( i ) ) (2–28) where g ( i ) = g ( x ( i ) ) and H ( i ) = H ( x ( i ) ) .Next,deningthesearchdirection, p ( i +1) = x ( i +1) x ( i ) ,theapproximationto J ( x ( i +1) ) canthenbewrittenas J ( x ( i +1) ) J ( x ( i ) )+ g ( i ) T p ( i +1) + 1 2 p ( i +1) T H ( i ) p ( i +1) (2–29) Inordertondaminimumof J ( x ) ,theTaylorseriesexpansioninEq.( 2–29 )is differentiatedwithrespectto x ( i +1) andthisderivativeissettozerosuchthat g ( i +1) = 0 = g ( i ) + H ( i ) p ( i +1) (2–30) RearrangingEq.( 2–30 ),theNewtonsearchdirectionisgivenas p ( i +1) = H ( i ) 1 g ( i ) (2–31) Finally,toensurethatthesearchisforaminimumandnotast ationarypointora maximum,thedescentconditionisenforced, g ( i ) T p ( i +1) < 0. (2–32) Eq.( 2–31 )isthenappliediterativelyfromsomeinitialguess x ( i ) suchthatanewpoint x ( i +1) isfound.If, J ( x ( i +1) ) < J ( x ( i ) ), (2–33) 27

PAGE 28

foreachiteration, i ,thenEq.( 2–31 )isappliediterativelyuntilaminimumisfound. 2.2.2EqualityConstrainedMinimizationofaFunction Considerndingtheminimumoftheobjectivefunction J ( x ) (2–34) subjecttoasetof m n constraints f ( x )= 0 (2–35) Findingtheminimumofanobjectivefunctionsubjecttoequa lityconstraintsusesan approachsimilartothecalculusofvariationsapproachfor determiningtheextremalof functionals.DenetheLagrangianas L ( x )= J ( x ) T f ( x ) (2–36) where T =( 1 ,..., m ) aretheLagrangemultipliers.Then,anecessaryconditionf or theminimumoftheLagrangianisthatthepoint ( x ) satises r x L ( x )= 0 (2–37) r L ( x )= 0 (2–38) Furthermore,thegradientof L withrespectto x and is r x L = g ( x ) G T ( x ) (2–39) r L = f ( x ) (2–40) 28

PAGE 29

wheretheJacobianmatrix, G ( x ) ,isdenedas G ( x ) @ f @ x = 266666664 @ f 1 @ x 1 @ f 1 @ x 2 ... @ f 1 @ x n @ f 2 @ x 1 @ f 2 @ x 2 ... @ f 2 @ x n ... @ f m @ x 1 @ f m @ x 2 ... @ f m @ x n 377777775 (2–41) ItisnotedthatataminimumoftheLagrangian,theequalityc onstraintofEq.( 2–35 )is satised.Next,thisnecessaryconditionalonedoesnotspe cifyaminimum,maximumor saddlepoint.Inordertospecifyaminimum,rst,denetheH essianoftheLagrangian as H L = r xx L = r xx J m X i =1 i r xx f i (2–42) Then,asufcientconditionforaminimumisthat v T H L v > 0 (2–43) foranyvector v intheconstrainttangentspace. InordertoiterativelyndaminimumoftheLagrangian,atwo -termTaylorseries expansionofEqs.( 2–39 )and( 2–40 )aboutsomepoint ( x ( i ) ( i ) ) isequatedtozero giving 0 = g ( i ) G ( i ) T ( i ) + H ( i ) L ( x ( i +1) x ( i ) ) G ( i ) T ( ( i +1) ( i ) ) (2–44) 0 = f ( i ) G ( i ) ( x ( i +1) x ( i ) ) (2–45) where G ( i ) = G ( x ( i ) ) and H ( i ) L = H L ( x ( i ) ) .Aftersimplication,theseequationsleadtothe followinglinearsystemofequations: 264 H ( i ) L G ( i ) T G ( i ) 0 375 264 p ( i +1) ( i +1) 375 = 264 g ( i ) f ( i ) 375 (2–46) 29

PAGE 30

ItisnotedthatEq.( 2–46 )iscalledtheKarush-Kuhn-Tucker(KKT)system[ 19 ]which isimplementedinanNLPsuchasSNOPT[ 11 ].Eq.( 2–46 )issolvediterativelyuntila minimumtotheLagrangianisfound.ItisnotedthatEq.( 2–46 )onlydenesanextremal point.Inordertondalocallyminimizingpoint,Eq.( 2–43 )mustalsobesatised.In determiningaminimumtotheLagrangian,aminimumto J ( x ) subjecttotheconstraints f ( x )= 0 isfound. 2.2.3InequalityConstrainedMinimizationofaFunction Considertheproblemofminimizingtheobjectivefunction J ( x ) (2–47) subjecttotheinequalityconstraints c ( x ) 0 (2–48) Inequalityconstrainedproblemsaresolvedbydividingthe inequalityconstraintsintoa setofactiveconstraints,andasetofinactiveconstraints .Atthesolution x ,someofthe constraintsaresatisedasequalities,thatis c i ( x )=0, i 2 A (2–49) where A iscalledtheactiveset,andsomeconstraintsarestrictlys atised,thatis c i ( x ) < 0, i 2 A 0 (2–50) where A 0 iscalledtheinactiveset. Byseparatingtheinequalityconstraintsintoanactiveseta ndaninactiveset,the activesetcanbedealtwithasequalityconstraintsasdescr ibedintheprevioussection, andtheinactivesetcanbeignored.Theaddedcomplexityofi nequalityconstrained problemsisindeterminingwhichsetofconstraintsareacti ve,andwhichareinactive.If theactiveandinactivesetsareknown,aninequalityconstr ainedproblembecomesan 30

PAGE 31

equalityconstrainedproblemstatedasminimizetheobject ivefunction J ( x ) (2–51) subjecttotheconstraints c i ( x )=0, i 2 A (2–52) andthesamemethodologyisappliedasintheprevioussectio ntodeterminea minimum.2.2.4InequalityandEqualityConstrainedMinimizationofa Function Finally,considertheproblemofndingtheminimumoftheob jectivefunction J ( x ) (2–53) subjecttotheequalityconstraints f ( x )= 0 (2–54) andtheinequalityconstraints c ( x ) 0 (2–55) DenetheLagrangianas L ( x ( A ) )= J ( x ) T f ( x ) ( A ) T c ( A ) ( x ) (2–56) where aretheLagrangemultiplierswithrespecttotheequalityco nstraintsand ( A ) aretheLagrangemultipliersassociatedwiththeactiveset ofinequalityconstraints. Furthermore,notethat ( A 0 ) = 0 ,andtherefore,theinactivesetofofconstraintsare ignored.AnecessaryconditionfortheminimumoftheLagran gianisthatthepoint ( x ( A ) ) satises r x L ( x ( A ) )= 0 (2–57) r L ( x ( A ) )= 0 (2–58) r ( A ) L ( x ( A ) )= 0 (2–59) 31

PAGE 32

Next,applyingatwo-termTaylorseriesapproximationtoth enecessaryconditionsof Eqs.( 2–57 )–( 2–59 )gives 0 = g ( i ) G ( i ) T f ( i ) G ( i ) T c ( A ) ( i )( A ) + H ( i ) L ( x ( i +1) x ( i ) ) (2–60) G ( i ) T f ( ( i +1) ( i ) ) G ( i ) T c ( A ) ( ( i +1)( A ) ( i )( A ) ) 0 = f ( i ) G ( i ) f ( x ( i +1) x ( i ) ) (2–61) 0 = c ( i )( A ) G ( i ) c ( A ) ( x ( i +1) x ( i ) ) (2–62) where G f ( x ) istheJacobianwithrespecttotheequalityconstraintsand G c ( A ) ( x ) isthe Jacobianwithrespecttotheactiveinequalityconstraints .Furthermore, H L inthiscaseis denedas H L = r 2xx L = r 2xx J m X i =1 i r 2xx f i q X i =1 ( A ) i r 2xx c ( A ) i (2–63) where q isthenumberofactiveinequalityconstraints.Therefore, theresultingKKT systemisgivenas 266664 H ( i ) L G ( i ) T f G ( i ) T c ( A ) G ( i ) f 00 G ( i ) c ( A ) 00 377775 266664 p ( i +1) ( i +1) ( i +1)( A ) 377775 = 266664 g ( i ) f ( i ) c ( i )( A ) 377775 (2–64) Fortheproblemofminimizinganobjectivefunctionsubject toequalityandinequality constraints,theKKTsystemofEq.( 2–64 )issolvediterativelysuchthataminimizing points ( x ( A ) ) isfound.Againitisnotedthatfor ( x ( A ) ) tobeaminimizing point,Eq.( 2–43 )mustbesatisedaswell. 2.2.5Finite-DimensionalOptimizationAsAppliedtoOptim alControl NLPsareusedtoapproximateacontinuous-timeoptimalcontr olproblembyadirect collocationmethodsuchthatthenite-dimensionaloptimi zationtechniquesofSection 2.2 canbeusedandthedifcultproblemofsolvingthecontinuou s-timerst-order necessaryconditionsisavoided.TheKKTsystemofEq.( 2–46 )orEq.( 2–64 )issolved iterativelyuntilaminimumisfoundfortheobjectivefunct ion.Therearetwoprimary 32

PAGE 33

considerationsinconstructingNLPsthatapproximateconti nuous-timeoptimalcontrol problems.Therstandmostimportantconsiderationisifth ediscreteapproximationwill provideanaccurateapproximationtothesolutionofaconti nuous-timeoptimalcontrol problem.Thesecondconsiderationisthatifthediscreteap proximationisaccurately approximatingthesolutionofacontinuous-timeoptimalco ntrolproblem,isitalsoa computationallyefcientproblemtosolve.IftheKKTsystemo fEq.( 2–46 )orEq.( 2–64 ) growslarge,thatis,thenumberofvariablesand/ornumbero fconstraintsislarge, iterativelysolvingEq.( 2–46 )orEq.( 2–64 )willbecomeexcessivelycomputationally expensive.Furthermore,iftheHessianorJacobianmatrice sareverydense,that is,alargenumberoftheirentriesarenon-zero,thenagain, theproblemwillbecome excessivelycomputationallyexpensive.Inordertomainta incomputationaltractability oftheNLP,theKKTsystemshouldbeconstructedtobesimultane ouslyassmalland sparseaspossiblewhileprovidinganaccurateapproximati onofthecontinuous-time optimalcontrolproblem. 2.3NumericalQuadrature Twoconceptsareimportantinconstructingthediscretized nite-dimensional optimizationproblemofthisresearchfromthecontinuoustimeoptimalcontrolproblem. Thesetwoconceptsarenumericalquadrature(numericalint egration)andpolynomial interpolation.Numericalquadraturecangenerallybedesc ribedassamplingan integrandatanitenumberofpointsandusingaweightedsum ofthesepointsto approximatetheintegral.Inthissection,rst,low-order numericalintegratorsarebriey discussed.Next,Gaussianquadraturerulesarepresenteda nditisshownthatthe quadraturerulesassociatedwiththeLGpointsarethemosta ccuratequadraturerules. Moreover,thequadraturerulesassociatedwiththeLGRpoin tsarethesecondmost accuratequadraturerules. 33

PAGE 34

2.3.1Low-OrderNumericalIntegrators Acommontechniqueusedtoapproximatetheintegralofafunc tionistouse manylow-degreeapproximatingsubintervalsthatareeache asilyintegrated.The approximationtotheintegralofafunctionisthenthesumof theintegralofeach approximatingsubinterval.Awell-knownlow-ordermethod isthecompositetrapezoidal rule[ 62 ].Thecompositetrapezoidruledividesafunctionintomany uniformlydistributed subintervalsandeachsubintervalisapproximatedbyastra ightlineapproximationthat passesthroughthefunctionattheendsofthesubinterval.F ig. 2-2 shows exp( ), 2 [ 1,1], andathreeintervalcompositetrapezoidruleapproximatio nofthefunction.For f ( )=exp( )f ( )Approximation 0 0 1.5 2 2.5 3 1 0.5 0.5 0.5 1 1 Figure2-2. e vs. andaThreeIntervalTrapezoidRuleApproximation. N approximatingintervals,thecompositetrapezoidruleisg ivenas[ 62 ] Z t f t 0 f ( t ) dt t f t 0 2 N ( f ( t 0 )+2 f ( t 1 )+2 f ( t 2 )+ +2 f ( t N 1 )+ f ( t N )) (2–65) 34

PAGE 35

where t T =( t 0 ,..., t N ) arethegridpointssuchthatpoints ( t i 1 t i ) denethebeginning andendofthe i th interval.Considerusingthecompositetrapezoidruletoap proximate Z 1 1 exp( ) d (2–66) Fig. 2-3 showsthe log 10 errorfromthetrapezoidruleforapproximatingEq.( 2–66 )asa functionofthe log 10 numberofapproximatingintervals.AsseeninFig. 2-3 ,inorderto achievehighlevelsofaccuracyintheapproximation,manya pproximatingintervalsare required.For onehundredthousand approximatingintervals,theerroris O (10 10 ) 2 2.5 3 4 5 4 6 8 10 12 log 10 Errorlog 10 NumberofApproximatingIntervals n 3.54.5 Figure2-3.Errorvs.NumberofIntervalsforTrapezoidRuleAp proximationof Eq.(2-66). 2.3.2GaussianQuadrature Incontrasttolow-orderpiecewiseapproximatingmethodss uchasthecomposite trapezoidrule,anotherquadraturemethodistousehighlya ccurateGaussianquadratures. 35

PAGE 36

Inusingan n pointGaussianquadrature,twoquantitiesmustbespecied foreach point,namely,thelocationofthepoint,and,aweightingfa ctormultiplyingthevalue ofthefunctionevaluatedatthispoint.Thethreesetsofpoi ntsdenedbyGaussian quadraturearetheLegendre-Gauss(LG),Legendre-Gauss-R adau(LGR),andthe Legendre-Gauss-Lobatto(LGL)points.FortheLGpoints,th epointsarefreetobe placedanywhereonthedomain 2 [ 1,1] ,fortheLGR,oneendpointisxedat 1 = 1 andfortheLGL,bothendpointsarexedat ( 1 n )=( 1,1) .Thegoalofa Gaussianquadratureistoconstructan n pointapproximation n X i =1 w i f i (2–67) where w i isthequadratureweightassociatedwithpoint i and f i = f ( i ) ,suchthat theapproximationinEq.( 2–67 )isexactforpolynomials f ( ) ofaslargeadegree aspossible[ 62 ].Let E n ( f ) betheerrorbetweentheintegralof f ( ) andan n point quadratureapproximationtotheintegralof f ( ) ,thatis, E n ( f )= Z 1 1 f ( ) d n X i =1 w i f i (2–68) First,itisnotedthatfor E n ( a 0 + a 1 + + a m m )= a 0 E n (1)+ a 1 E n ( )+ + a m E n ( m ) (2–69) that E n ( f )=0 foreverypolynomialofdegree m ifandonlyif E n ( i )=0,( i =0,1,..., m ). (2–70) Next,considertheconstructionofaquadraturetobeexactf orapolynomialofas largeadegreeaspossiblewhilerequiring 1 = 1 (thatis,considertheconstruction oftheLGRpoints).Twoexamplesareprovidedtodemonstrate howthisquadratureis constructed. 36

PAGE 37

Example1: n =1 :For n =1 ,thereisonefreeparameter, w 1 ,and 1 = 1 isxed. Becausethereisonefreeparameter,thegoalistohavezeroer rorforazerothdegree polynomial,i.e., E n (1)=0. (2–71) FromEq.( 2–71 ),thereisoneequationtodeterminetheunknownweight w 1 givenas Z 1 1 1 d w 1 =0. (2–72) Inordertosatisfythisequation, w 1 =2 Example2: n =2 :For n =2 ,therearethreefreeparameters, w 1 w 2 and t 2 and again t 1 = 1 isxed.Becausetherearethreefreeparameters,itisdesire dtohave zeroerrorforaseconddegreepolynomial,i.e., E n (1)=0, (2–73) E n ( )=0, (2–74) E n ( 2 )=0. (2–75) Thisleadstothesetofequations w 1 + w 2 =2, (2–76) w 1 + w 2 t 2 =0, (2–77) w 1 + w 2 t 2 2 = 2 3 (2–78) InsolvingEqs.( 2–76 )and( 2–78 )for w 1 w 2 ,and t 2 itisfoundthat w 1 =0.5, w 2 =1.5, and 2 =1 = 3 .As n isincreased,asimilarapproachisused; 2 n 1 freeparameters aredeterminedsuchthatthequadratureisexactforpolynom ialsofdegree 2 n 2 .In constructingsetsofpointsinthisfashion,theLGRquadrat urepointsandweightsare constructed.Similarly,fortheLGmethod,because t 1 isnotxed,for n points,theLG quadraturerulehas 2 n freeparameterstodetermineinordertoconstructquadratu re 37

PAGE 38

approximationsaccurateforfunctionsofdegreeuptoandin cluding 2 n 1 .Asa result,foragivennumberofpoints,theLGmethodprovidest hehighestaccuracy quadratureapproximationofanyquadraturerules,whileth eLGRmethodprovidesthe nexthighestprecision.TheLGLquadraturerule,for n points,has 2 n 2 freeparameters andareaccurateforpolynomialsofdegreeuptoandincludin g 2 n 3 .Itshouldbe notedonceagainthattheLGandLGRpointsarestudiedinthis researchbecausethe LGLpointshavepreviouslybeendemonstratedtoresultinin accuratepseudospectral approximationstocontinuous-timeoptimalcontrolproble m[ 40 42 ].As n becomes large,determiningtheLG,LGR,orLGLpointsandweightsbec omesincreasingly moredifcultbythemethodexempliedintheaforementione dtwoexamples.The resultingsetofequationsarenonlinearandtheirsolution sarenotobvious.Fortunately, ithasalsobeenshownthattheLG,LGR,andLGLpointsarealso denedbyalinear combinationofaLegendrepolynomialand/oritsderivative s[ 63 ].The n LGpointsare denedastherootsofan n th –degreeLegendrepolynomial P n ( )= 1 2 n n d n d n ( 2 1) n (2–79) andthecorresponding n weightsaregivenas w i = 2 1 2 i P n ( i ) 2 ,( i =1,..., n ) (2–80) where P n isthederivativeofthe n th –degreeLegendrepolynomial.Furthermore,the n LGRpointsaredenedastherootsofthesumofthe n th and ( n 1) th –degreeLegendre polynomials,i.e., P n ( )+ P n 1 ( ) ,andthecorrespondingquadratureweightsforthe LGRpointsare w 1 = 2 n 2 (2–81) w i = 1 (1 i ) P n 1 ( i ) 2 ,( i =2,..., n ). (2–82) 38

PAGE 39

Figs. 2-4 and 2-5 showtheLGandLGRpointsforvariousvaluesof n .First,itis notedthatfortheLGpointstheendpoints, =( 1,+1) ,arenotincludedandfor theLGRpointstheinitialpoint, = 1 ,isincluded.Furthermore,bothsetsofpoints demonstrateasimilarproleofhavingmorepointsnearthee ndsoftheinterval.Finally, itisnotedthattheLGpointsaresymmetricabouttheorigina ndtheLGRpointsare asymmetric. 0 0 1 0.5 0.5 1 5 10 NumberofLGPoints15 20 25 30 Figure2-4.Legendre-GaussPointsDistribution. TheaccuracyoftheLGandLGRpointsarenowdemonstratedona pproximating theintegralgiveninEq.( 2–66 ).Fig. 2-6 showsthe log 10 errorfortheapproximation ofEq.( 2–66 )versusthenumberofLGandLGRpoints.Itisseenthatasthen umber ofLGorLGRpointsisincreased,theaccuracyinapproximati onrapidlyincreases. Furthermore,therateofconvergencebytheLGandLGRquadra turestotheexact solutionisexponential.ForeightLGorLGRpoints,theappr oximationisexactwith 39

PAGE 40

0 0 1 0.5 0.5 1 5 10 15 20 25 30 NumberofLGRPoints f Figure2-5.Legendre-Gauss-RadauPointsDistribution.respecttonumericalroundoff.Whereasthetrapezoidrulere quiredonehundred thousandapproximatingintervals(i.e.,onehundredthous andandoneevaluationsof exp( ) )toapproximatetheintegralof exp( ) to O (10 10 ) ,only six pointsarerequiredto achievethissamelevelofaccuracyusingLGorLGRpoints(i. e.,only six evaluations of exp( ) ).Therefore,forapproximatingtheintegralofEq.( 2–66 ),theLGorLGR quadraturesaresignicantlymoreaccuratethanalow-orde rmethodsuchasthe compositetrapezoidrule.Finally,itisnotedthattheLGpo intsresultinaquadrature approximationthatisslightlymoreaccuratethantheLGRpo intsforapproximating Eq.( 2–66 ).ThisisconsistentwiththefactthattheLGquadratureise xactfora polynomialonedegreehigherthantheLGRpointsforanygive nnumberofquadrature points. 40

PAGE 41

0 2 4 2 4 6 8 10 12 14 16 6 8 10 1214 16 NumberofLGPointslog 10 Error LGPts LGRPts Figure2-6.Errorvs.NumberofLGPointsforApproximationofEq .(2-66). 2.3.3Multiple-IntervalGaussianQuadrature Gaussianquadraturemayalsobeimplementedinamultiple-i ntervalsensesuch thatahighlyaccurateGaussianquadratureisappliedtoeac happroximatinginterval. WhilethepointsdenedbyGaussianquadratureareonthedoma in 2 [ 1,1] ,itis notedthattheintegralofafunctionoveranarbitrarytimed omain t 2 [ t a t b ] maybe approximatedas Z t b t a f ( t ) dt = t b t a 2 Z 1 1 f ( t b t a 2 + t a + t b 2 ) d t b t a 2 n X i =1 w i f ( t b t a 2 i + t a + t b 2 ). (2–83) Next,considerapproximatingEq.( 2–66 )againusinguniformlydistributedmultiple-interval LGorLGRquadratures.Theapproximationto exp( ), 2 [ 1,1] using K uniformly 41

PAGE 42

distributedintervalsis Z +1 1 exp( ) d K X k =1 t k t k 1 2 n X i =1 w i f ( t k t k 1 2 i + t k + t k 1 2 ) (2–84) where t k 1 and t k denethebeginningandendofthe k th approximatinginterval, respectively,and n arethenumberofLGorLGRquadraturepointsineachinterval Figs. 2-7A and 2-7B showthe log 10 errorversusthetotalnumberofLGorLGRpoints, Kn ,forLGandLGRquadratureapproximationsusing n =(2,3,4) .InFigs. 2-7A and 2-7B ,itisseenthathighlevelsofaccuracyareobtainedbyusing amultiple-interval Gaussianquadrature.Furthermore,for n =2 ,twohundredLGpointsresultinan approximationofaccuracy O (10 10 ) asopposedtotheonehundredthousandpoints necessarytoobtainthislevelofaccuracyusingthetrapezo idrule.Furthermore,as n is increased,fewertotalpointsarerequiredtoobtainhighle velsofaccuracy.For n =3 theapproximationisexactfor150LGpoints,andfor n =4 ,theapproximationisexact for40LGpoints.InFig. 2-7B itisseenthatthemultiple-intervalLGRquadratureis lessaccuratethanthemultiple-intervalLGquadrature.Bec ausetheLGquadratureis accurateforapolynomialofonedegreehigher,forlowandmi d-degreeapproximations, thedifferenceinaccuracyisquitenoticeable.For n =4 ,themultiple-intervalLG quadratureisexactfor40LGpoints,whereastheLGRquadrat urerequires100points beforeprovidinganexactapproximationwithrespecttonum ericalroundoff. 2.3.4Summary Numericalintegrationmethodshavebeendiscussed.First, thecompositetrapezoid ruleisgivenasanexampleofalow-orderapproximationmeth od.Furthermore,this ruleisusedtoapproximateEq.( 2–66 )anditisseenthatmanyapproximatingintervals arerequiredtoachievehighlevelsofaccuracy.Next,Gauss ianquadraturerules arediscussed.Gaussianquadraturerulesareconstructeds uchthatpolynomials ofashighadegreeaspossibleareapproximatedusingasetof n pointsand n weightingfactors.Itisshownthatveryhighlevelsofaccur acyareobtainedforvery 42

PAGE 43

0 0 10 100 50 150 200 5 15 NumberofMulti-IntervalLGPointslog 10 Error n =2 n =3 n =4 ALGPoints. 0 0 10 100 50 150 200 5 15 Maxim n =2 n =3 n =4 NumberofMulti-IntervalLGRPointslog 10 Error BLGRPoints. Figure2-7.Errorvs.NumberofMultiple-IntervalLGorLGRPo intsforApproximationof Eq.( 2–66 ). 43

PAGE 44

fewpointsbyusingLGorLGRpointstoapproximatetheintegr alofEq.( 2–66 ). Finally,multiple-intervalGaussianquadraturearediscu ssed.Highlevelsofaccuracy inquadratureapproximationareobtainedagainforthiscas eforafewnumberof approximatingpointswhencomparedwiththecompositetrap ezoidrule.Furthermore, theLGquadratureismoreaccuratethantheLGRquadrature.I tisnotedthatwhile morepointsarerequiredtoaccuratelyapproximateEq.( 2–66 )usingamultiple-interval Gaussianquadraturethanasingleinterval,athoroughmoti vationisgivenastowhya multiple-intervalGaussianquadratureisdesiredinconst ructingthe hp –pseudospectral methodofthisresearchinChapter 3 andChapter 4 2.4LagrangePolynomialApproximation Thenextimportantconceptusedinconstructingthediscret izednite-dimensional optimizationproblemofthisresearchispolynomialinterp olation.Inthisresearch, Lagrangepolynomialsareusedasacontinuous-timeapproxi mationtothestate equationsofthecontinuous-timeoptimalcontrolproblem. InaLagrangepolynomial approximation,asetof n supportpoints, ( t 1 ,..., t n ) ,areusedandapolynomialof degree n 1 isttedtothesepoints.ALagrangepolynomialapproximati ontoafunction, f ( t ) ,isgivenby ~ f ( t )= n X j =1 f j L j ( t ), L j ( t )= n Y l =1 l 6 = j t t l t j t l (2–85) where f j = f ( t j ) .Itisnotedthatatthesupportpoints t j for j =(1,..., n ) ~ f i = f i 2.4.1Single-IntervalLagrangePolynomialApproximations ItiswellknownthatRungephenomenonexistsforLagrangepo lynomialsforan arbitrarysetofsupportpoints.Rungephenomenonistheocc urrenceofoscillationsin theapproximatingfunction, ~ f ,betweensupportpoints.Asthenumberofsupportpoints grows,themagnitudeoftheoscillationsbetweensupportpo intsalsogrows.Runge phenomenoncanbedemonstratedonthefunction f ( )= 1 1+25 2 2 [ 1,+1]. (2–86) 44

PAGE 45

Fig. 2-8 showstheLagrangepolynomialapproximationtothisfuncti onutilizing n = (8,32) uniformlydistributedsupportpoints.Itisseenthatasthe numberofpointsare increased,theapproximationneartheendsoftheintervalb ecomeincreasingly worse Next,consideraLagrangepolynomialapproximationtothes amefunctionutilizing n =(8,32) LGsupportpoints.ThisapproximationisshowninFig. 2-9 anditisseen thatbyusingLGsupportpoints,Rungephenomenonisavoided ,andtheapproximations become better asthenumberofsupportpointsareincreased.Insightintot hebehavior oftheLagrangepolynomialapproximationsisexempliedby thefollowingtheorem[ 62 ]. Theorem :Let t 0 t 1 ,..., t n bedistinctrealnumbers,andlet f beagivenrealvalued functionwith n +1 continuousderivativesontheinterval I t = H ( t t 0 t 1 ,..., t n ) ,with t somegivenrealnumber.Thenthereexists 2 I t suchthat f ( t ) ~ f ( t )= ( t t 0 ) ( t t n ) ( n +1)! f n +1 ( ). (2–87) Eq.( 2–87 )denesalimitontheerrorbetweenthefunctionbeingappro ximatedanda Lagrangepolynomialapproximationofthefunction.Itshou ldbenotedatanysupport point,theerroriszero.Indesigningan n pointLagrangepolynomialapproximation, theonlydesignfreedomisinplacementofthesupportpoints ,i.e.,changingthe behaviorofthenumeratorterminEq.( 2–87 ).Foruniformlydistributedpoints,the numeratortermbecomesexcessivelylargeattheendsofthei ntervals.Thisresults intheRungephenomenonthatisseeninFig. 2-8 .Furthermore,pointsdenedby Gaussianquadraturehaveastructurethattheyaremoredens elyplacedneartheends oftheinterval(seeagainFigs. 2-4 and 2-5 ).Thisdistributionreducesthemagnitude ofthenumeratorterminEq.( 2–87 )attheendsoftheintervals,therefore,reducing theeffectsofRungephenomenon.Fig. 2-11 showsthe log 10 maximumerrorasa functionofdegreeofapproximatingLagrangepolynomialfo runiformlydistributed supportpointsandsupportpointsdenedbyLGandLGRpoints forapproximating Eq.( 2–86 ).Itisseenthatasthedegreeofpolynomialincreases,theL agrange 45

PAGE 46

0 0 1 0.5 0.5 1 1 0.2 0.4 0.6 0.8 f ( ) ~ f ( ) f ( ) A8UniformSupportPoints. 0 0 1 0.5 0.5 1 200 400 600 800 200 Maxim f ( ) ~ f ( )f ( ) B32UniformSupportPoints. Figure2-8.Approximationof f ( )=1 = (1+25 2 ) UsingUniformlySpacedSupport Points. 46

PAGE 47

1.2 0.2 0 0 1 0.5 0.5 1 1 0.2 0.4 0.6 0.8 f ( ) f ( ) ~ f ( ) A8LGSupportPoints. 0 0 1 0.5 0.5 1 1 0.2 0.4 0.6 0.8 Each f ( ) f ( ) ~ f ( ) B32LGSupportPoints. Figure2-9.Approximationof f ( )=1 = (1+25 2 ) UsingLGSupportPoints. 47

PAGE 48

polynomialapproximationusingauniformlydistributedse tofsupportpoints diverges Fora 100 th degreepolynomialdenedwithuniformlydistributedsuppo rtpoints,errors are O (10 30 ) .ForaLagrangepolynomialdenedbyeitherLGorLGRsupport points,the approximation converges tothefunctioninEq.( 2–86 ).Fora 100 th degreepolynomial, theerrorsinapproximatingthisfunctionare O (10 10 ) .Furthermore,unlikethecase ofusingLGandLGRpointsfortheapproximationtointegrals ,itisnotedthatwhen usedasthesupportpointstoapproximatethefunctioninEq.( 2–86 )byaLagrange polynomial,thereisvirtuallynodifferenceinapproximat ionerrorbyLGorLGRpoints. 2.4.2Multiple-IntervalLagrangePolynomialApproximati ons ConsiderapproximatingEq.( 2–86 )usingmultiple-intervalapproximationsoflow ormediumdegree.WithoutusingLGorLGRpoints(oranyothers etofpointsthat minimizesRungephenomenoninaLagrangepolynomialapprox imation),commonly, Rungephenomenonisavoidedbysimplylimitingthedegreeof approximationtobe small.Considerauniformlydistributedsetofpolynomiala pproximationssuchthatthe degreeofapproximationisone,i.e.,linear,ineachinterv al.Furthermore,denethe supportpointsforeachintervaltobethebeginningandendo ftheinterval.Fig. 2-11 showsthe log 10 maximumerrorinapproximatingEq.( 2–86 )asafunctionof log 10 numberofapproximatingintervals.Itisseenthatforeveno nehundredthousand approximatingintervals,errorsarestill O (10 8 ) Next,considerauniformlydistributedmultiple-interval approximationutilizingLG andLGRpointsasthesupportpointsineachintervaltoappro ximateEq.( 2–86 ).Itis notedthattheLGandLGRpointscanbetransformedfromtheti medomain 2 [ 1,1] toanarbitrarytimedomain t 2 [ t a t b ] bythetransformation t = t b t a 2 + t a + t b 2 (2–88) Using n =(4,6,8) LGandLGRsupportpointsineachinterval,Fig. 2-12B shows the log 10 maximumerrorinapproximatingEq.( 2–86 )versusnumberofapproximating 48

PAGE 49

2 2.5 3 4 5 2 4 6 8 10 f 3.54.5log 10 MaximumErrorlog 10 NumberofApproximatingInterval Figure2-10.Errorvs.NumberofLinearApproximatingInterva lsofEq.( 2–86 ). intervalsusingaLagrangepolynomialdenedby n =(4,6,8) LGorLGRsupportpoints ineachinterval.Itisseenthathighlevelsofaccuracyareo btainedforarelativelyfew numberofintervals.Furthermore,asthenumberofsupportp ointsareincreasedin eachinterval,theapproximationbecomesincreasinglyacc urateforagivennumberof intervals.Finally,itisagainnotedthatwhenusingeither LGorLGRpointsassupport pointsthereisvirtuallynodifferenceintheerrorofappro ximation. 2.4.3Summary AdiscussiononLagrangepolynomialapproximationtoafunc tionisgiven. Rungephenomenonisdemonstratedtoexistwhenusinghigh-d egreepolynomial approximationstoafunctionusingasetofuniformlydistri butedsupportpoints. Furthermore,itisdemonstratedthatRungephenomenondoes notexistwhenusingLG 49

PAGE 50

0 0 10 10 40 60 80100 20 20 30 f NumberofSupportPointslog 10 MaximumErrorUniformPts LGRPts LGPts Figure2-11.Errorvs.NumberofSupportPointsforApproximati ngEq.( 2–86 ). orLGRsupportpointsandhigh-accuracyapproximationsare obtainedforhigh-degree Lagrangepolynomialsapproximations.Next,inapproximat ingafunctionusingmultiple intervals,itisshownthataverylow-orderapproximationr equiresatremendousnumber ofapproximatingintervalstoachievehighlevelsofaccura cy.Byusingmoderately sizedmultiple-intervalLagrangepolynomialapproximati onsusingLGorLGRpoints, itisdemonstratedthathighlevelsofaccuracyinapproxima tingafunctioncanbe obtainedformanyfewerapproximatingintervalsthanwhenu singverylow-degree approximations. 50

PAGE 51

0 0 10 100 50150 200 5 15 NumberofIntervalswith n PointsInEachIntervallog 10 MaximumError n =4 n =6 n =8 ALGPoints. 0 0 10 100 50150 200 5 15 NumberofIntervalswith n PointsInEachIntervallog 10 MaximumError n =4 n =6 n =8 BLGRPoints. Figure2-12.Errorvs.NumberofLGorLGRIntervalsforApproxi mationofEq.( 2–86 ). 51

PAGE 52

CHAPTER3 MOTIVATIONFORAN HP –PSEUDOSPECTRALMETHOD Amotivationisgivenfordevelopingan hp –pseudospectralmethod.Themotivation istocombinethecomputationalsparsityof h –methodsandtherapidratesofconvergence of p –methods.Inordertomotivatethedesireforan hp –approach,twoexamples withdistinctlydifferentsolutionsareprovided.Inther stexample,aproblemwitha smoothsolutionisstudied.Inthesecondexample,aproblem whosesolutionhas adiscontinuouscontrolisanalyzed.Becauseofthedisconti nuityinthecontrol,the solutiontothissecondproblemisnon-smooth. 3.1MotivatingExample1:ProblemwithaSmoothSolution Considerthefollowingoptimalcontrolproblem[ 64 ].Minimizethecostfunctional J = t f (3–1) subjecttothedynamicconstraints x =cos( ),_ y =sin( ), = u (3–2) andtheboundaryconditions x (0)=0, y (0)=0, (0)= x ( t f )=0, y ( t f )=0, ( t f )= (3–3) ThesolutiontotheoptimalcontrolproblemofEqs.( 3–1 )–( 3–3 )is ( x ( t ), y ( t ), ( t ), u ( t ))=( sin( t ), 1+cos( t ), t ,1) (3–4) where t f =2 .Thisproblemisnowsolvedwithasingle-interval p –LGRmethod [ 40 ]andauniformlydistributedmultiple-interval h –LGRmethodthatisdescribedin Chapter 4 .The h –LGRapproximationusestwocollocationpointsineachinte rval, i.e.,aseconddegreepolynomialapproximationtothestate .Becausethisproblemhas asmoothsolution,itisexpectedthata p –methodwillperformbetterwhencompared 52

PAGE 53

toan h –method.Fig. 3-1A showsthe log 10 maximumerroratthecollocationpoints versusthetotalnumberofcollocationpointsforthesingle -interval p –LGRmethod.It isseenthatfor20collocationpoints,theapproximationis essentiallyexact(relativeto roundofferror).Next,Fig 3-1B showsthe log 10 maximumerroratthecollocationpoints versusthe log 2 numberofapproximatingintervalsforauniformlydistribu ted h –LGR approximation.Usingthe h –LGRmethodtoapproximatethesolution,theconvergence ratesaresignicantlyslowerthanusinga p –method.Foranaccuracyof O (10 8 ) ,the h –LGRmethodrequires1024approximatingintervalsand2048 collocationpoints whilethe p –methodusesasingleapproximatingintervalwithonly13co llocationpoints. Strictlylow-ordermethodsareunabletoobtainhighlevelso faccuracyforthisproblem inacomputationallyefcientmanner.Forthisproblem,hig hlevelsofaccuracyinthe approximationaremosteasilyachievedbyusinga p –method. 3.2MotivatingExample2:ProblemwithaNonsmoothSolution Considerthefollowingoptimalcontrolproblemofasoftlun arlandingasgivenin Ref.[ 65 ].Minimizethecostfunctional J = Z t f 0 udt (3–5) subjecttothedynamicconstraints h = v ,_ v = g + u (3–6) theboundaryconditions ( h (0), v (0), h ( t f ), v ( t f ))=( h 0 v 0 h f v f )=(10, 2,0,0), (3–7) andthecontrolinequalityconstraint u min u u max (3–8) 53

PAGE 54

0 0 5 5 10 10 5 15 15 20 25 log 10 MaximumError NumberofCollocationPoints A p –LGRApproximation. 0 0 2 4 2 4 6 8 10 6 810 Numberof log 10 MaximumError log 2 NumberofCollocationPoints B h –LGRApproximation. Figure3-1.Errorvs.NumberofCollocationPointsfor p and h –LGRApproximationto Eqs.( 3–1 )–( 3–3 ). 54

PAGE 55

where u min =0 u max =3 g =1.5 t f isfree.Theoptimalsolutiontothisexampleis ( h ( t ), v ( t ), u ( t ))= 8><>: ( 3 4 t 2 + v 0 t + h 0 3 2 t + v 0 ,0), t < t s ( 3 4 t 2 +( 3 t s + v 0 ) t + 3 2 t s 2 + h 0 3 2 t +( 3 t s + v 0 ),3), t > t s (3–9) where t s is t s = t f 2 + v 0 3 (3–10) with t f = 2 3 v 0 + 4 3 r 1 2 v 2 0 + 3 2 h 0 (3–11) FortheboundaryconditionsgiveninEq.( 3–7 ), t s =1.4154 and t f =4.1641 .First, itisnotedinEq.( 3–9 )thatthecontrolisbang-bang,i.e.,thecontrolisatitsmi nimum valuefor t 2 [0, t s ] andisatitsmaximumvaluefor t 2 [ t s t f ] .Becauseofthe discontinuityinthecontrol,theexactsolutionforthesta teisnon-smooth.Fig. 3-2 shows the log 10 maximumerroratthecollocationpointsversusthenumberof collocationpoints forapproximatingthisproblemby p anduniformlydistributed h –LGRmethods.Itis seenthatbecausethesolutiontotheproblemisnon-smooth, the p –methoddoesnot convergeatanexponentialrate.Furthermore,theratesofc onvergenceofthe p and h –methodsareessentiallythesame.Therefore,whencompari ng p and h –methods withrespecttonumberofcollocationpoints,nosignicant differenceexistsbetween a p or h –methodapproximation.ThesparsityoftheapproximatingN LPsusing p and h –methodsare,however,quitedifferent.A p –methodresultsinamuchmoredense NLPascomparedtoan h –method.Table 3-1 showsthedensityoftheapproximating NLPby p and h –methodsforagivennumberofcollocationpoints.Furtherm ore,the CPUtimerequiredtosolvethegivenproblemisalsoshown.Iti sseeninTable 3-1 that foranequivalentnumberofcollocationpoints,the p –methodresultsinaproblemthat issignicantlymorecomputationallyexpensivebecauseof thedensityoftheNLP.As thenumberofcollocationpointsgrowforthe p –methodapproximation,theCPUtime 55

PAGE 56

0 0 1 1 2 4 100 50150 200 Error 3 log 10 MaximumError NumberofCollocationPoints p h – 2 Figure3-2.Errorvs.CollocationPtsfor p –and h –LGRApproximationto Eqs.( 3–5 )–( 3–8 ). requiredtosolvetheNLPrapidlygrowsaswell.Forthe h –method,becauseofthe sparsityoftheNLP,thetimerequiredtosolvetheproblemdo esnotgrowasrapidly aswiththe p –method.Furthermore,itisnotedthatfor200collocationp oints,the h –methodsolutiontimeisnearlyequivalenttothesolutiont imeofthe40collocation point p –method.Therefore,iftheexponentialconvergenceofa p –methodcannotbe guaranteed,amuchmorecomputationallyefcientNLPisobt ainedusingan h –method approximation. 3.3DiscussionofResults Thetwomotivatingexamplesdemonstrateakeypointinappro ximatingthe solutiontogeneraloptimalcontrolproblemsbydirectcoll ocationmethods:thebest approximationmethodis problemdependent .Furthermore,insomecasesthatare 56

PAGE 57

CollocationPoints p –MatrixDensity h –MatrixDensity p –CPUTime(s) h –CPUTime(s) 20 0.37 0.10 0.07 0.06 40 0.35 0.05 0.22 0.11 60 0.35 0.037 0.38 0.12 80 0.34 0.028 0.54 0.14 100 0.34 0.023 0.88 0.16 200 0.34 0.012 2.97 0.19 Table3-1.MatrixDensityandCPUTimefor p –and h –LGRMethodsforExample2. simpleandsmooth,a p –methodapproximationwillprovidethebestapproximating NLP. Inthecaseofanon-smoothproblemwithouthigherorderbeha viorsinthesolution,an h –methodapproximationwillprovidethebestapproximating NLP.Inthegeneralcase, asolutiontoanoptimalcontrolproblemwillhavemixedbeha viorswheresomeregions ofthetrajectorywillbesmooth,andotherregionswillbeno n-smooth.Therefore,a methodthatallowsforthecombinationofthestrengthofeit her h or p –methodswill becapableofprovidinghighlyaccuratesolutionsatalowco mputationalcost.An hp –pseudospectralmethodhastheexibilityofusingthecomp utationalsparsityof low-order h –methodsandthehighconvergenceratesof p –methods.The hp –LGand LGRmethodsarenowdescribedinChapter 4 57

PAGE 58

CHAPTER4 HP –PSEUDOSPECTRALMETHOD The hp –LGandLGRpseudospectralmethodsofthisresearcharenowd escribed. First,thecontinuous-timeoptimalcontrolproblemofSecti on 2.1 isreformulatedina multiple-intervalformandthecontinuous-timerst-orde rnecessaryconditionsarestated forthismultiple-intervalformulation.The hp –LGandLGRpseudospectraldiscretizations ofthemultiple-intervaloptimalcontrolproblemarethenp osed.TheKKTconditions arederivedforthe hp –LGandLGRmethods.Then,thetransformedadjointsystemfo r eachmethodisderivedfromtheKKTconditionsandamappingfro mthetransformed adjointsystemofeachmethodtothecontinuous-timerst-o rdernecessaryconditionsis derived.Finally,thesparsityoftheNLPdenedbythe hp –LGandLGRpseudospectral methodsisstudied. 4.1Multiple-IntervalContinuous-TimeOptimalControlPro blemFormulation Consideragainthecontinuous-timeoptimalcontrolproble mofSection 2.1 .Thatis, minimizethecostfunctional J =( x ( t 0 ), t 0 x ( t f ), t f )+ Z t f t 0 g ( x ( t ), u ( t )) dt (4–1) subjecttothedynamicconstraints x ( t )= f ( x ( t ), u ( t )), (4–2) theinequalitypathconstraints C ( x ( t ), u ( t )) 0 (4–3) andtheboundaryconditions ( x ( t 0 ), t 0 x ( t f ), t f )= 0 (4–4) Next,considerdividingtheproblemdenedinEqs.( 4–1 )–( 4–4 )into K meshintervals.A meshinterval k isdenedtobeginatthemeshpoint t k 1 andendatthemeshpoint t k 58

PAGE 59

Furthermore,let t 0 < t 1 < t 2 < < t K where t K = t f .Thetimedomainineachinterval k t 2 [ t k 1 t k ] ,ischangedto 2 [ 1,1] bythetransformation = 2 t ( t k 1 + t k ) t k t k 1 (4–5) Furthermore,theinversetransformationisgivenas t = t k t k 1 2 + t k 1 + t k 2 (4–6) anditisnotedthat dt d = t k t k 1 2 (4–7) Next,leteachmeshpoint t k beexpressedasafunctionoftheinitialtime, t 0 ,naltime, t K ,andaconstant a k k =(0,..., K ) ,as t k = t 0 + a k ( t K t 0 ),( k =0,..., K ). (4–8) Itisnotedthatfor t 0 a 0 =0 andfor t K a K =1 .Furthermore,because t k 1 < t k a k 1 < a k for k =1,..., K .Byexpressingeachmeshpointasafunctionof t 0 and t K thenumberoffreevariablescorrespondingtotimeisreduce dfrom K +1 freevariables correspondingtoeachmeshpointtotwofreevariablescorre spondingtotheinitialand naltimes.Thisisimportantforconstructingaconciseset ofnecessaryconditionsfor theoptimalcontrolproblemandforreducingthenumberofva riablesinthediscretized NLP.Let x ( k ) ( ) and u ( k ) ( ) bethestateandcontrolinthe k th meshinterval.Finally, itisnotedthatthecontinuityinthestateisenforcedatthe interiormeshpointsand therefore, x ( k ) (+1)= x ( k +1) ( 1),( k =1,..., K 1) Thecontinuous-timeoptimalcontrolproblemdenedbyEqs.( 4–1 )–( 4–4 ) expressedasa K intervalproblemisthentominimizethecost J =( x (1) ( 1), t 0 x ( K ) (+1), t K )+ K X k =1 t k t k 1 2 Z +1 1 g ( x ( k ) ( ), u ( k ) ( )) d (4–9) 59

PAGE 60

subjecttothedynamicconstraints d x ( k ) ( ) d = t k t k 1 2 f ( x ( k ) ( ), u ( k ) ( )),( k =1,..., K ) (4–10) theinequalitypathconstraints t k t k 1 2 C ( x ( k ) ( ), u ( k ) ( )) 0 ,( k =1,..., K ) (4–11) theboundaryconditions ( x (1) ( 1), t 0 x ( K ) (+1), t K )= 0 (4–12) andtheinteriorpointconstraints x ( k ) (+1) x ( k +1) ( 1)= 0 ,( k =1,..., K 1). (4–13) Itisnotedthat t k t k 1 2 (4–14) isaddedtoEq.( 4–11 )withoutaffectingtheconstraintsuchthattheoptimality conditions canbeposedinasuccinctmanner. Therst-orderoptimalityconditionsofthemultiple-inte rvalcontinuous-timeoptimal controlproblemofEqs.( 4–9 )–( 4–13 )arenowstated.DenetheaugmentedHamiltonian ineachinterval k as H ( k ) ( x ( k ) ( k ) u ( k ) r ( k ) )= g ( k ) + h ( k ) f ( k ) ih r ( k ) C ( k ) i (4–15) where h i isusedtodenotethestandardinnerproductbetweentwovect ors.Usingthis denitionoftheaugmentedHamiltonianineachinterval,th econtinuous-timerst-order 60

PAGE 61

optimalityconditionsoftheproblemdenedbyEqs.( 4–9 )–( 4–13 )are d x ( k ) d = t k t k 1 2 r ( k ) H ( k ) ,( k =1,..., K ), (4–16) d ( k ) d = t k t k 1 2 r x ( k ) H ( k ) ,( k =1,..., K ), (4–17) 0 = r u ( k ) H ( k ) ,( k =1,..., K ), (4–18) (1) ( 1)= r x (1) ( 1) ( h i ) (4–19) ( k ) (1)= ( k +1) ( 1)+ ( k ) ,( k =1,..., K 1) (4–20) ( K ) (+1)= r x ( K ) (+1) ( h i ) (4–21) 0= K X i =1 a k 1 a k 2 Z 1 1 H ( k ) d + r t 0 ( h i ) (4–22) 0= K X i =1 a k a k 1 2 Z 1 1 H ( k ) d + r t f ( h i ) (4–23) where ( k ) isthedifferencein ( k ) (1) and ( k +1) ( 1) ifthecostateisdiscontinuousata meshpoint k [ 2 ]. 4.2 hp -Legendre-GaussTranscription The hp –LGpseudospectralmethodisnowdescribed.The hp –LGpseudospectral methodisconstructedbydiscretizingthemultiple-interv alcontinuous-timeoptimal controlproblemdescribedbyEqs.( 4–9 )–( 4–13 ).First,thestateineachinterval k of thecontinuous-timeoptimalcontrolproblemisapproximat edbyaLagrangepolynomial usingtheinitialpoint, ( k ) 0 = 1 ,andtheLGpoints, ( ( k ) 1 ,..., ( k ) N k ) ,where N k isthe numberofLGpointsusedininterval k .Thestateapproximationisthengivenas x ( k ) ( ) X ( k ) ( )= N k X j =0 X ( k ) j L ( k ) j ( ), (4–24) where X ( k ) j = X ( k ) ( ( k ) j ) .Furthermore,itisnotedthatinthisChapter,allvectorfu nctions oftimeare row vectors,i.e., X ( k ) ( )=[ X ( k ) 1 ( ),..., X ( k ) n ( )] .Next,anapproximationto thederivativeofthestateisgivenbydifferentiatingthea pproximationofEq.( 4–24 )with 61

PAGE 62

respectto .Therefore,thederivativeofthestateapproximationisgi venas d X ( k ) ( ) d = N k X j =0 X ( k ) j L ( k ) j ( ). (4–25) Thecollocationconditionsforinterval k areformedbycollocatingthederivativeofthe stateapproximationwiththerighthandsideofthedynamicc onstraintsofEq.( 4–10 )at the N k LGpointsas N k X j =0 X ( k ) j D ( k ) ij t k t k 1 2 f ( k ) i = 0 ,( i =1,..., N k ), (4–26) where f ( k ) i = f ( X ( k ) i U ( k ) i ) U ( k ) i = U ( k ) ( ( k ) i ) isthediscretizedcontrolapproximation,and D ( k ) ij = L ( k ) j ( ( k ) i ),( i =1,..., N k j =1,..., N k +1) (4–27) isthe N k ( N k +1) Gausspseudospectraldifferentiationmatrix [ 22 26 35 36 ]inthe k th meshinterval.Itshouldbenotedthatthestateapproximati onisacontinuous-time Lagrangepolynomialapproximationwith N k +1 supportpointswhilethecontrolis calculateddiscretelyatthe N k collocationpoints.Theinequalitypathconstraintsof Eq.( 4–11 )areenforcedatthe N k collocationpointsineachinterval k as t k t k 1 2 C ( k ) ( X ( k ) i U ( k ) i ) 0 ,( i =1,..., N k ). (4–28) Next,itisnotedthatthusfarthestateapproximationonlyu sestheinitialpoint, 0 = 1 andthe N k LGpointsineachinterval.Nodiscreteapproximationasofy ethasbeen madeforthestateatthenalpointineachinterval, ( k ) N k +1 =+1 .Inordertoapproximate thestateatthenalpointineachinterval,thefollowingLG quadratureisused X ( k ) N k +1 = X ( k ) 0 + t k t k 1 2 N k X i =1 w ( k ) i f ( k ) i (4–29) 62

PAGE 63

where w ( k ) j aretheLegendre-Gaussweightsininterval k .Next,thestateisenforcedto becontinuousbetweenmeshintervals,i.e., X ( k ) N k +1 = X ( k +1) 0 ,( k =1,..., K 1). (4–30) Inordertoincreasethecomputationalefciencyoftheappr oximatingNLP,asingle variableisusedfor X ( k ) N k +1 and X ( k +1) 0 for k =(1,..., K 1) .Therefore,Eq.( 4–29 )isthen expressedas X ( k +1) 0 = X ( k ) 0 + t k t k 1 2 N k X i =1 w ( k ) i f ( k ) i ,( k =1,..., K 1) (4–31) X ( K ) N K +1 = X ( K ) 0 + t K t K 1 2 N K X i =1 w ( K ) i f ( K ) i andEq.( 4–30 )canbeignoredasitisbeingexplicitlytakenintoaccount. Finally,the costfunctionalofEq.( 4–9 )isapproximatedusingamultiple-intervalLGquadraturea s J ( X (1)0 t 0 X ( K ) N K +1 t K )+ K X k =1 N k X j =1 t k t k 1 2 w ( k ) j g ( k ) j (4–32) where g ( k ) j = g ( X ( k ) j U ( k ) j ) Thediscretized hp –LGpseudospectralapproximationtothecontinuous-time optimalcontrolproblemofEqs.( 4–9 )–( 4–13 )isnowstated.The hp –LGpseudospectral approximationistominimizethecost J =( X (1)0 t 0 X ( K ) N K +1 t K )+ K X k =1 N k X j =1 t k t k 1 2 w ( k ) j g ( k ) j (4–33) subjecttothecollocationconditions N k X j =0 X ( k ) j D ( k ) ij t k t k 1 2 f ( k ) i = 0 ,( i =1,..., N k k =1,..., K 1), (4–34) theinequalitypathconstraints t k t k 1 2 C ( k ) ( X ( k ) i U ( k ) i ) 0 ,( i =1,..., N k k =1,..., K 1), (4–35) 63

PAGE 64

thequadratureconstraints X ( k +1) 0 = X ( k ) 0 + t k t k 1 2 N k X i =1 w ( k ) i f ( k ) i ,( k =1,..., K 1) (4–36) X ( K ) N K +1 = X ( K ) 0 + t K t K 1 2 N K X i =1 w ( K ) i f ( K ) i andtheboundaryconditions ( X (1)0 t 0 X ( K ) N K +1 t K )= 0 (4–37) TheproblemdenedbyEqs.( 4–33 )–( 4–37 )isthediscrete hp –LGpseudospectral approximationtothecontinuous-timeoptimalcontrolprob lemdenedbyEqs.( 4–9 )–( 4–13 ). Itisnotedthattheformulationpresentedallowsforanarbi trarynumberofintervals,with arbitraryplacement,andwithanarbitrarynumberofLGcoll ocationpointsineach interval.Thatis,the hp –LGmethodpresentedhastheexibilityofapproximatinga continuous-timeoptimalcontrolproblemin any mannerdesiredbyusingLGpoints. Afewkeypropertiesofthe hp –LGpseudospectralmethodarenowstated.First,the discretizationpointsaredenedasthepointsatwhichthes tateisapproximated.The discretizationpointsinaninterval k ,therefore,arethe N k +2 points, ( ( k ) 0 ,..., ( k ) N k +1 ) Itisnotedthatcontrolisapproximatedatthe N k collocationpoints, ( ( k ) 1 ,..., ( k ) N k ) Next,awell-dened N th k –degreepolynomialapproximationtothestateisgivenby Eq.( 4–24 ).However,nowell-denedapproximatingpolynomialisgiv enforthecontrol. Itisinterestingtonotehowthediscretizationpointsandc ollocationpointsdifferfor the hp –LGpseudospectralmethod.Whenusingasingleapproximatin gintervalas showninFig. 4-1A ,thediscretizationpointsdifferfromthecollocationpoi ntsinthatthe discretizationpointsarethecollocationpoints plus theendpoints.Furthermore,fora multiple-intervaldiscretizationasisseeninFig. 4-1B ,thediscretizationpointsdifferfrom thecollocationpointsinthatthediscretizationpointsin cludethecollocationpointsand 64

PAGE 65

7 9 0 3 1 0.5 0.5 1 4 5 6 8 10 PolynomialDegreeCollocationPoints DiscretizationPoints t ASingle-IntervalCollocationPointsandDiscretizationPoints. 7 0 2 3 1 0.5 0.5 1 4 5 6 8 CollocationPoints DiscretizationPoints tNumberofIntervals BMultiple-IntervalCollocationPointsandDiscretizatio n Points. Figure4-1.CollocationandDiscretizationPointsbyLGMet hods. allthemeshpoints.Itisinterestingtonotethatthecolloc ationpointsdonotinclude any meshpointsforthe hp –LGpseudospectralmethod. 4.2.1The hp –LGPseudospectralTransformedAdjointSystem The hp –LGtransformedadjointsystemisnowderived.Itisshownth atthe hp –LG transformedadjointsystemisadiscreterepresentationof thecontinuous-timerst-order 65

PAGE 66

necessaryconditionsofEqs.( 4–16 )–( 4–23 )ifthecostateoftheoptimalcontrolproblem iscontinuous.First,theLagrangianofthecontinuous-tim eoptimalcontrolproblem denedbythe hp –LGpseudospectralmethodisgivenas L =( X (1)0 X ( K ) N K +1 ) h ( X (1)0 X ( K ) N K +1 i (4–38) + K X k =1 N k X j =1 t k t k 1 2 ( w ( k ) j g ( k ) j h ( k ) j C ( k ) j i ) h ( k ) j D ( k ) j ,1: N k X ( k ) 1: N k + D ( k ) j ,0 X ( k ) 0 t k t k 1 2 f ( k ) j i K 1 X k =1 h ( k ) X ( k +1) 0 X ( k ) 0 t k t k 1 2 N k X j =1 w ( k ) j f ( k ) j i h ( K ) X ( K ) N k +1 X ( K ) 0 t K t K 1 2 N K X j =1 w ( K ) j f ( K ) j i where ( k ) ( k ) ( k ) ,and arethemultipliersassociatedwithEqs.( 4–34 ),( 4–35 ), ( 4–36 ),and( 4–37 ).Furthermore, D ( k ) j ,1: N k isthe j th rowandrstthrough N th k columns of D ( k ) and X ( k ) 1: N k aretherstthrough N th k rowsof X ( k ) .Next,theKKTconditionsare givenbydifferentiatingtheLagrangianwithrespectto X ( k ) j j =1,..., N K +1 and U ( k ) j j =1,..., N k ineachinterval k aswellas t 0 and t K andsettingthesederivatives 66

PAGE 67

equaltozero.TheKKTconditionsare 0 = r U w ( k ) j g ( k ) j + h ( k ) j + w ( k ) j ( k ) f ( k ) j i (4–39) h ( k ) j C ( k ) j i ,( k =1,..., K j =1,..., N k ) D ( k ) T j ( k ) = t k t k 1 2 r X w ( k ) j g ( k ) j + h ( k ) j + w ( k ) j ( k ) f ( k ) j i (4–40) h ( k ) j C ( k ) j i ,( k =1,..., K j =1,..., N k ) D (1) T 0 (1) (1) = r X (1)0 ( h i ), (4–41) D ( k ) T 0 ( k ) ( k ) = ( k 1) ,( k =2,..., K ), (4–42) ( K ) = r X ( K ) N K +1 ( h i ), (4–43) 0 = r t 0 ( h i ) (4–44) + K X k =1 N k X j =1 a k 1 a k 2 ( w ( k ) j g ( k ) j + h ( k ) j f ( k ) j ih ( k ) j C ( k ) j i ) K 1 X k =1 a k 1 a k 2 h ( k ) N k X j =1 w ( k ) j f ( k ) j i + a K 1 a K 2 h K N K X j =1 w ( K ) j f ( K ) j i ( k =1,..., K j =1,..., N k ) 0 = r t f ( < > ) (4–45) + K X k =1 N k X j =1 a k a k 1 2 ( w ( k ) j g ( k ) j + h ( k ) j f ( k ) j ih ( k ) j C ( k ) j i ) + K 1 X k =1 a k a k 1 2 h ( k ) N k X j =1 w ( k ) j f ( k ) j i + a K a K 1 2 h K N K X j =1 w ( K ) j f ( K ) j i ( k =1,..., K j =1,..., N k ) where D Tj isthe j th rowof D T .IncomparingEq.( 4–41 )withEq.( 4–19 ),thecostate approximationat t = t 0 isgivenas (1)0 = D (1) T 0 (1) + (1) (4–46) Next,incomparingEq.( 4–43 )withEq.( 4–21 ),thecostateapproximationat t = t f is givenas ( K ) N K +1 = ( K ) (4–47) 67

PAGE 68

Considerthechangeofvariables = (4–48) r ( k ) j = ( k ) j w ( k ) j ,( k =1,..., K j =1,..., N k ) (4–49) ( k ) j = ( k ) w ( k ) j + ( k ) ,( k =1,..., K j =1,..., N k ), (4–50) anddenethematrix D ( k ) y [ 40 ]as D ( k ) y ij = w ( k ) j w ( k ) i D ( k ) ji ,( k =1,..., K j =1,..., N k ) (4–51) D ( k ) y i N k +1 = N k X j =1 D ( k ) y ij ,( i =1,..., N k ). (4–52) D ( k ) y denesthedifferentiationmatrixassociatedwithaLagran gepolynomialusingthe LGpoints, ( ( k ) 1 ,..., ( k ) N k ) andthe nalpoint ( k ) N k +1 inaninterval k [ 40 ].Substitutingthe changeofvariablesgiveninEqs.( 4–46 )–( 4–52 )into( 4–39 )–( 4–43 ),theLGtransformed 68

PAGE 69

adjointsystemisgivenas ( D ( k ) y j ,1: N k ( k ) + D ( k ) y j N k +1 ( k ) )= t k t k 1 2 r X H ( X ( k ) j U ( k ) j r ( k ) j ( k ) j ), (4–53) ( k =1,..., K j =1,..., N k ) 0 = r U H ( X ( k ) j U ( k ) j r ( k ) j ( k ) j ), (4–54) ( k =1,..., K j =1,..., N k ) (1)0 = r X (1)0 ( h i ), (4–55) ( K ) N K +1 = r X ( K ) N K +1 ( h i ), (4–56) ( k ) ( k 1) = t k t k 1 2 N k X j =1 w ( k ) j r X H ( X ( k ) j U ( k ) j r ( k ) j ( k ) j ), (4–57) ( k =2,..., K ), 0 = K X k =1 a k 1 a k 2 N k X j =1 w ( k ) j H ( X ( k ) j U ( k ) j r ( k ) j ( k ) j ) (4–58) + r t 0 ( h r i ), 0 = K X k =1 a k a k 1 2 N k X j =1 w ( k ) j H ( X ( k ) j U ( k ) j r ( k ) j ( k ) j ) (4–59) + r t f ( h r i ), whereitisnotedthatEq.( 4–57 )isformedbycomparingEqs.( 4–40 )and( 4–42 )andthe identity[ 40 ] D ( k ) i ,0 = N k X j =1 D ( k ) ij ,( i =1,..., N k ). (4–60) FromEq.( 4–57 ),notingthat ( K ) N K +1 = ( K ) ( K ) 0 = ( K 1) (4–61) Supposethecostateiscontinuousatmeshpoint K 1 .Thisthenimplies ( K ) 0 = ( K 1) N k +1 andusingEq.( 4–57 )again,weobtain ( K 1) 0 = ( K 2) (4–62) 69

PAGE 70

Next,supposethecostateiscontinuousatallthemeshpoint s k =1,..., K 1 .Thenthe costateatthemeshpointsaregivenas ( k +1) 0 = ( k ) N k +1 = ( k ) ,( k =1,..., K 1). (4–63) Inthecasewherethecostateiscontinuousateverymeshpoin t,thetransformed adjointsystemisadiscreterepresentationofthecontinuo us-timerst-ordernecessary conditions. Next,supposethatthecostateisdiscontinuousatameshpoi nt M 2 [1,..., K 1] Then ( M ) = ( M +1) 0 = ( M ) N M +1 ( M ) (4–64) where ( M ) isthedifferencebetween ( M ) N M +1 and ( M +1) 0 .Ifthecostateisdiscontinuous atmeshpoint M ( M ) 6 = ( M ) N M +1 inEq.( 4–53 ).Therefore,Eq.( 4–53 )isnotadiscrete approximationtothecontinuous-timecostatedynamics.Asa result,thetransformed adjointsystemgiveninEqs.( 4–53 )–( 4–59 )isnotadiscreterepresentationofthe rst-orderoptimalityconditionsgiveninEqs.( 4–16 )–( 4–23 ). 4.2.2Sparsityofthe hp –LGPseudospectralMethod Thesparsityofthe hp –LGpseudospectralmethodisnowdescribed.Itisnotedthat theonlydenseblocksofnon-zeronumbersthatoccurintheNL Pdenedbythe hp –LG pseudospectralmethodarearesultoftheLagrangepolynomi alapproximationtothe stateineachintervalinEq.( 4–34 )andalsothequadratureconstraintsofEq.( 4–36 ). Thesizeofthesedenseblocksaredependentuponthedegreeo fapproximating polynomial, N k ,ineachinterval.TheconstraintJacobianforthe hp –LGpseudospectral methodisnowdescribedindetail.First,denethetotalnum berofcollocationpointsas N = K X i =1 N i (4–65) Next,thetotalvectorofNLPvariablescorrespondingtothe discretizedstate,control, andinitialandnaltimeisconstructed.Denethe ( N + K +1) n 1 vector z x asthe 70

PAGE 71

vectorcontainingeachcomponentofthestateatallthedisc retizationpoints,i.e., z x = 266664 1 ... n 377775 (4–66) where i isthe ( N + K +1) 1 vectorcorrespondingtothe i th componentofthestate evaluatedatthediscretizationpoints.Furthermore,den ethe Nm 1 vector z u asthe vectorcontainingeachcomponentofthecontrolevaluateda tthecollocationpoints,i.e., z u = 266664 1 ... m 377775 (4–67) where i isthe N 1 vectorcorrespondingtothe i th componentofthecontrolatallof thecollocationpoints.Therefore,thetotalvectorofNLPv ariablesis z = 266666664 z x z u t 0 t f 377777775 (4–68) andthetotalnumberofNLPvariablesis ( N + K +1) n + Nm +2 Next,thetotalvectorofNLPconstraintsareconstructed.D enethe Nn 1 vector y D asthevectorcontainingthecollocationconditionsforeac hcomponentofthestate, i.e., y D = 266664 1 ... n 377775 (4–69) where i isthe N 1 vectorcorrespondingtothecollocationconditionsforthe i th componentofthestate.Next,denethe Ns 1 vector y C asthevectorcontainingthe 71

PAGE 72

inequalitypathconstraints: y C = 266664 1 ... s 377775 (4–70) where i isthe N 1 vectorcorrespondingtothe i th inequalitypathconstraint.Denethe Kn 1 vector y Q tobethevectorofquadratureconstraintsinEq.( 4–36 ).Finally,dene the q 1 vector y E tobethevectorofboundaryconditionsinEq.( 4–37 ).Therefore,the totalvectorofconstraintsforthe hp –LGmethodis y = 266666664 y D y C y Q y E 377777775 (4–71) andthetotalnumberofNLPconstraintsis Nn + Ns + Kn + q .Itisinterestingtonote forthe hp –LGpseudospectralmethod,thesizeoftheNLPproblem(i.e. ,thenumber ofvariablesandconstraints)aredependentuponthenumber ofcollocationpoints, N andthenumberofapproximatingintervals, K .Asthenumberofcollocationpoints or numberofapproximatingintervalsgrowslargeforthe hp –LGpseudospectralmethod,so dothenumberofvariablesandconstraintsintheNLP.Finall y,theconstraintJacobian forthe hp –LGpseudospectralmethodisgivenas G ( z )= @ y @ z (4–72) ThedensityoftheconstraintJacobianisnowexamined.Fort heconstraints associatedwith y C ,eachrowcorrespondstoaninequalitypathconstrainteval uated atacollocationpoint.ByinspectionofEq.( 4–35 ),itisseenthatateachcollocation point,theinequalitypathconstraintsareonlyafunctiono fthestateandcontrolatthat collocationpointandtheinitialandnaltime.Therefore, therowsof y C aresparse(i.e.,, only n + m +2 elementsinanyrowarenon-zero;theremaining ( N + K ) n +( N 1) m 72

PAGE 73

elementsarezero).Furthermore, y E isonlyafunctionofthestateattheinitialandnal timeandtheinitialandnaltime.Therefore,therowsof y E aresparseaswell. y D and y Q introducethedensityofthe hp –LGpseudospectralmethod.Itisseenin Eq.( 4–34 )thatthecollocationconditionforanycomponentofthesta teevaluatedatany collocationpointisafunctionofthestateandcontrolatth atcollocationpoint,theinitial andnaltime,and N k additionalvariablescorrespondingtotheLagrangepolyno mial approximationtothestate.Therefore,as N k growslargeinanyapproximatinginterval, thenumberofnon-zeroelementsinanyrowof y D growslargeaswell.Finally,itis seeninEq.( 4–36 )thatforeachquadratureconstraint,theconstraintisafu nctionof N k n variablescorrespondingtothestate, N k m variablescorrespondingtothecontrol, theinitialandnaltime,andtwoadditionalvariablescorr espondingtotheinitialstate andnalstateininterval k .Therefore,again,as N k growslargeinanyapproximating intervals,sodoesthenumberofnon-zeroentriesinanyrowo f y Q Figs. 4-2 and 4-3 showsaqualitativerepresentationoftheconstraintJacob ian forthe hp –LGpseudospectralmethod.Fig. 4-2 representsan hp –LGpseudospectral approximationwhere N k islargeineachapproximatingintervalandFig. 4-3 represents an hp –LGpseudospectralapproximationwhere N k issmallineachapproximating interval.Itisseenthatforlarge N k ,large,denseblocksariseintheconstraintJacobian correspondingtotherowsof y D and y Q .Forsmallvaluesof N k ,thedenseblocksare muchsmaller.Therefore,thesparsityofan hp –LGpseudospectralmethodisdependent upon N k ,or,thedegreeofpolynomialapproximationwithinanyinte rval.Low-degree approximationsresultinasparseNLPwhilehigh-degreeapp roximationsresultina denseNLP. 4.3 hp –Legendre-Gauss-RadauTranscription The hp –LGRpseudospectralmethodisnowdescribed.The hp –LGRpseudospectral methodisconstructedbydiscretizingthecontinuous-time optimalcontrolproblem describedbyEqs.( 4–9 )–( 4–13 ).First,thestateineachinterval k ofthecontinuous-time 73

PAGE 74

Figure4-2.QualitativeRepresentationofconstraintJaco bianforLGpointswithLarge N k optimalcontrolproblemisapproximatedbyaLagrangepolyn omialusingthe N k LGR points, ( ( k ) 1 ,..., ( k ) N k ) ,andthenalpoint, ( k ) N k +1 =+1 .Itisnotedthat ( k ) 1 = 1 .The stateapproximationisthengivenas x ( k ) ( ) X ( k ) ( )= N k +1 X j =1 X ( k ) j L ( k ) j ( ). (4–73) 74

PAGE 75

Figure4-3.QualitativeRepresentationofconstraintJaco bianforLGpointswithSmall N k Next,thederivativeofthestateapproximationisgivenbyd ifferentiatingtheapproximation ofEq.( 4–73 )withrespectto .Thederivativeofthestateapproximationisthengivenas d X ( k ) ( ) d = N k +1 X j =1 X ( k ) j L ( k ) j ( ). (4–74) 75

PAGE 76

Thecollocationconditionsarethenenforcedatthe N k LGRpointsas N k +1 X j =1 X ( k ) j D ( k ) ij t k t k 1 2 f ( k ) i = 0 ,( i =1,..., N k k =1,..., K ), (4–75) where D ( k ) ij = L ( k ) j ( ( k ) i ),( i =1,..., N k j =1,..., N k +1), (4–76) isthe N k ( N k +1) Radaupseudospectraldifferentiationmatrix [ 39 ]inthe k th mesh interval.TheinequalitypathconstraintsofEq.( 4–11 )areenforcedatthe N k collocation pointsineachinterval k as t k t k 1 2 C ( k ) ( X ( k ) i U ( k ) i ) 0 ,( i =1,..., N k k =1,..., K ). (4–77) Next,thestateisenforcedtobecontinuousbetweenmeshint ervals,i.e., X ( k ) N k +1 = X ( k +1) 1 (4–78) Inordertoincreasethecomputationalefciencyoftheappr oximatingNLP,asingle variableisusedfor X ( k ) N k +1 and X ( k +1) 1 for k =(1,..., K 1) .Therefore,Eq.( 4–75 )isthen writtenas D ( k ) j ,1: N k X ( k ) 1: N k + D ( k ) j N k +1 X ( k +1) 1 t k t k 1 2 f ( k ) j = 0 ( k =1,..., K 1), (4–79) D ( K ) j ,1: N K X ( K ) 1: N K + D ( K ) j N K +1 X ( K ) N K +1 t K t K 1 2 f ( K ) j = 0 andEq.( 4–78 )canbeignoredasitisbeingexplicitlytakenintoaccount. Finally, thecostfunctionalofEq.( 4–9 )isthenapproximatedusingamultiple-intervalLGR quadratureas J ( X (1)1 t 0 X ( K ) N K +1 t K )+ K X k =1 N k X j =1 t k t k 1 2 w ( k ) j g ( k ) j (4–80) where w ( k ) j aretheLegendre-Gauss-Radauweights[ 63 ]ininterval k 76

PAGE 77

Thediscretized hp –LGRpseudospectralapproximationtothecontinuous-time optimalcontrolproblemisnowstated.The hp –LGRpseudospectralapproximationisto minimizethecost J =( X (1)1 t 0 X ( K ) N K +1 t K )+ K X k =1 N k X j =1 t k t k 1 2 w ( k ) j g ( k ) j (4–81) subjecttothecollocationconditions D ( k ) j ,1: N k X ( k ) 1: N k + D ( k ) j N k +1 X ( k +1) 1 t k t k 1 2 f ( k ) j = 0 ( k =1,..., K 1), (4–82) D ( K ) j ,1: N K X ( K ) 1: N K + D ( K ) j N K +1 X ( K ) N K +1 t K t K 1 2 f ( K ) j = 0 theinequalitypathconstraints t k t k 1 2 C ( k ) ( X ( k ) i U ( k ) i ) 0 ,( i =1,..., N k k =1,..., K ). (4–83) andtheboundaryconditions ( X (1)1 t 0 X ( K ) N K +1 t K )= 0 (4–84) TheproblemdenedbyEqs.( 4–81 )–( 4–84 )isthediscrete hp –LGRpseudospectral approximationtothecontinuous-timeoptimalcontrolprob lemdenedbyEqs.( 4–9 )–( 4–13 ). Itisnotedthattheformulationpresentedallowsforanarbi trarynumberofintervals,with arbitraryplacement,andanarbitrarynumberofLGRcolloca tionpointsineachinterval. Therefore,thismethodhastheexibilityofapproximating acontinuous-timeoptimal controlproblemin any mannerdesiredbyusingLGRpoints. Afewkeypropertiesofthe hp –LGRpseudospectralmethodarenowstated.First, thediscretizationpointsaredenedasthepointsatwhicht hestateisapproximated. Thediscretizationpointsinaninterval k ,therefore,are ( ( k ) 1 ,..., ( k ) N k +1 ) .Itisnotedthat thecontrolisapproximatedatthecollocationpoints, ( ( k ) 1 ,..., ( k ) N k ) .Next,an N th k degree polynomialapproximationtothestateisgivenbyEq.( 4–73 ).Thecontrolontheother hand,isdiscretelyapproximatedatthe N k collocationpointswithineachinterval.The 77

PAGE 78

discretizationpointsandcollocationpointsarenowcompa redforthe hp –LGRmethod. WhenusingasingleapproximatingintervalasisshowninFig. 4-4A ,thediscretization pointsdifferfromthecollocationpoints only atthenalmeshpoint, t K .Furthermore, foramultiple-intervaldiscretizationasseeninFig. 4-4B ,again,thediscretization pointsdifferfromthecollocationpointsonlyatthenalme shpoint.Whilethe hp –LG pseudospectralmethoddidnotapproximatethecontrolatan ymeshpoints,the hp –LGR pseudospectralmethodapproximatesthecontrolateveryme shpointexceptthenal meshpoint.4.3.1The hp -LGRPseudospectralTransformedAdjointSystem TheKKTconditionsoftheNLPdenedbythe hp –LGRpseudospectralmethodare nowderived.TheLagrangianofthe hp -LGRpseudospectralmethodisgivenas L =( X (1)1 X ( K ) N K +1 ) h ( X (1)1 X ( K ) N K +1 i (4–85) + K X k =1 N k X j =1 t k t k 1 2 ( w ( k ) j g ( k ) j h ( k ) j C ( k ) j i ) K 1 X k =1 N k X j =1 h ( k ) j D ( k ) j ,1: N k X ( k ) 1: N k + D ( k ) j N k +1 X ( k +1) 1 t k t k 1 2 f ( k ) j i N K X j =1 h ( K ) j D ( K ) j ,1: N K X ( K ) 1: N K + D ( K ) j N K +1 X ( K ) N K +1 t K t K 1 2 f ( K ) j i where ( k ) ( K ) ,and arethemultipliersassociatedwithEq.( 4–82 ),Eq.( 4–83 ),and Eq.( 4–84 )inmeshinterval k ,respectively.TheKKTconditionsofthemultiple-interval Legendre-Gauss-Radaudiscretizationarethengivenas D (1) T j (1) = t 1 t 0 2 r X w (1) j g (1) j + h (1)j f (1) j i (4–86) h (1)j C (1)j i 1 j r X (1)1 + r X (1)1 h i ,( j =1,..., N 1 ), D ( k ) T j ( k ) = t k t k 1 2 r X w ( k ) j g ( k ) j + h ( k ) j f ( k ) j i (4–87) h ( k ) j C ( k ) j i 1 j D ( k 1) T N k 1 +1 ( k 1) ,( k =2,..., K j =1,..., N k ). D ( K ) T N K +1 ( K ) = r X ( K ) N K +1 ( h i ), (4–88) 78

PAGE 79

7 9 0 2 3 1 0.5 0.5 1 4 5 6 8 10 DiscretizationPoints CollocationPointsPolynomialDegreet ASingle-IntervalCollocationPointsandDiscretizationP oints. 7 0 2 3 1 0.5 0.5 1 4 5 6 8 DiscretizationPoints CollocationPoints tNumberofIntervals BMultiple-IntervalCollocationPointsandDiscretizatio nPoints. Figure4-4.CollocationandDiscretizationPointsbyLGRMe thods. 0 = r U w ( k ) j g ( k ) j + h ( k ) j f ( k ) j ih ( k ) j C ( k ) j i (4–89) ( k =1,..., K j =1,..., N k ), 0 = r t 0 ( < > ) (4–90) + K X k =1 N k X j =1 a k a k 1 2 ( w ( k ) j g ( k ) j + h ( k ) j f ( k ) j ih ( k ) j C ( k ) j i ) 0 = r t f ( < > ) (4–91) + K X k =1 N k X j =1 a k a k 1 2 ( w ( k ) j g ( k ) j + h ( k ) j f ( k ) j ih ( k ) j C ( k ) j i ) 79

PAGE 80

Next,considerthechangeofvariables = (4–92) r ( k ) j = ( k ) j w ( k ) j ,( k =1,..., K j =1,..., N k ), (4–93) ( k ) j = ( k ) j w ( k ) j ,( k =1,..., K j =1,..., N k ), (4–94) ( k ) = D ( k ) T N k +1 ( k ) ,( k =1,..., K ). (4–95) Furthermore,denethematrix D ( k ) y [ 40 ]as D ( k ) y 11 = D ( k ) 11 1 w ( k ) 1 (4–96) D ( k ) y ij = w ( k ) j w ( k ) i D ( k ) ji otherwise (4–97) D ( k ) y isthecorrespondingdifferentiationmatrixforaLagrange polynomialapproximation usingthe N k LGRpointsassupportpointsininterval k .SubstitutingEqs.( 4–92 )–( 4–97 ) intoEqs.( 4–86 )–( 4–91 ),thetransformedadjointsystemisgivenas D (1) y j (1) = t 1 t 0 2 r X H ( X (1)j U (1)j (1)j r (1)j ) (4–98) + 1 j w (1) 1 ( r X (1)1 ( h i ) (1)1 ), D ( k ) y j ( k ) = t k t k 1 2 r X H ( X ( k ) j U ( k ) j ( k ) j r ( k ) j ) (4–99) + 1 j w ( k ) 1 ( ( k 1) ( k ) 1 ),( k =2,..., K j =1,..., N k ) 0 = r U H ( X ( k ) j U ( k ) j ( k ) j r ( k ) j ),( k =1,..., K j =1,..., N k ) (4–100) ( K ) = r X ( K ) N K +1 ( h i ) (4–101) 0 = K X k =1 a k 1 a k 2 N k X j =1 w ( k ) j H ( X (1)j U (1)j (1)j r (1)j ) (4–102) + r t 0 ( h r i ), 0 = K X k =1 a k a k 1 2 N k X j =1 w ( k ) j H ( X (1)j U (1)j (1)j r (1)j ) (4–103) + r t 0 ( h r i ), 80

PAGE 81

where D ( k ) y j isthe j th rowof D ( k ) y .FromtheresultsofRef.[ 39 ],weexpandEq.( 4–95 ) usingEqs.( 4–86 )and( 4–87 )toobtainthefollowingrelationships: r X (1)1 (+ h i )= (1) + t 1 t 0 2 N 1 X j =1 w (1) j r x H ( X (1)j U (1)j (1)j r (1)j ). (4–104) ( k 1) = ( k ) + t k t k 1 2 N k X j =1 w ( k ) j r X H ( X ( k ) j U ( k ) j ( k ) j r ( k ) j ).( k =2,..., K ) (4–105) BycomparingEq.( 4–21 )andEq.( 4–101 ),thecostateapproximationat t = t f isgiven as ( K ) N K +1 = ( K ) (4–106) For k = K ,theright-handsideofEq.( 4–105 )isthenanapproximationto ( K ) 1 Therefore, ( K ) 1 = ( K 1) (4–107) Furthermore,usingEq.( 4–107 ),thenalterminEq.( 4–99 )disappearsfor k = K .Next, letthecostatebecontinuousatmeshpoint k = K 1 .Then, ( K 1) N K 1 +1 = ( K ) 1 = ( K 1) (4–108) andfor k = K 1 ,therighthandsideofEq.( 4–105 )becomesanapproximationto ( K 1) 1 .Therefore, ( K 1) 1 = ( K 2) (4–109) andthenalterminEq.( 4–99 )disappearsfor k = K 1 .Ifthecostateiscontinuousat everymeshpoint,then ( k ) N k +1 = ( k +1) 1 = ( k ) ,1 k K 1 (4–110) andthenalterminEq.( 4–99 )disappearsfor k =2,..., K .Furthermore,thelefthand sideofEq.( 4–104 ),then,isanapproximationtothecostateat t 0 andtherefore, (1)1 = r X (1)1 (+ h i ) (4–111) 81

PAGE 82

andthenalterminEq.( 4–98 )disappears.Therefore,forthecasewhenthecostate iscontinuousateverymeshpoint,thenalterminEqs.( 4–98 )and( 4–99 )disappears andthetransformedadjointsystemisadiscreterepresenta tionofthecontinuous-time rst-orderoptimalityconditions. Next,supposethatthecostateisdiscontinuousataparticu larmeshpoint M 2 [1,..., K 1] .Then ( M ) = ( M +1) 1 = ( M ) N M +1 ( M ) (4–112) where ( M ) isthedifferencebetween ( M ) N M +1 and ( M +1) 1 .Therefore,ifthecostate isdiscontinuousatmeshpoint M ( M 1) 6 = ( M ) 1 inEq.( 4–99 ).Asaresult,the transformedadjointsystemgiveninEqs.( 4–98 )–( 4–103 )isnotadiscreterepresentation oftherst-orderoptimalityconditionsgiveninEqs.( 4–16 )–( 4–23 ). 4.3.2Sparsityofthe hp –LGRPseudospectralMethod Thesparsityofthe hp –LGRpseudospectralmethodisnowdescribed.Itisnoted thattheonlydenseblocksofnon-zeronumbersthatoccurint heNLPdenedbythe hp –LGRpseudospectralmethodarearesultoftheLagrangepoly nomialapproximation tothestateineachintervalinEq.( 4–82 ).Thesizeofthesedenseblocksaredependent uponthedegreeofapproximatingpolynomial, N k ,ineachinterval.Theconstraint Jacobianforthe hp –LGRpseudospectralmethodisnowdescribedindetail.Thet otal vectorofNLPvariablescorrespondingtothediscretizedst ate,controlandinitialand naltimeisconstructed.Denethe ( N +1) n 1 vector z x asthevectorcontainingeach componentofthestateatallthediscretizationpoints,i.e ., z x = 266664 1 ... n 377775 (4–113) where n isthe ( N +1) 1 vectorcorrespondingtothe i th componentofthestateatthe discretizationpoints.Furthermore,denethe Nm 1 vector z u asthevectorcontaining 82

PAGE 83

eachcomponentofthecontrolatallthecollocationpoints, i.e., z u = 266664 1 ... m 377775 (4–114) where i isthe N 1 vectorcorrespondingtothe i th componentofthecontrolatthe collocationpoints.Therefore,thetotalvectorofNLPvari ablesis z = 266666664 z x z u t 0 t f 377777775 (4–115) andthetotalnumberofNLPvariablesis ( N +1) n + Nm +2 .Next,thetotalvectorof NLPconstraintsareconstructed.Denethe Nn 1 vector y D asthevectorcontaining thecollocationconditionsforeachcomponentofthestate, i.e., y D = 266664 1 ... n 377775 (4–116) where i isthe N 1 vectorcorrespondingtothecollocationconditionsforthe i th componentofthestate.Next,denethe Ns 1 vector y C asthevectorcontainingthe theinequalitypathconstraints: y C = 266664 1 ... s 377775 (4–117) where i isthe N 1 vectorcorrespondingtothe i th inequalitypathconstraint.Finally, denethe q 1 vector y E tobethevectorofboundaryconditionsinEq.( 4–84 ). 83

PAGE 84

Therefore,thetotalvectorofconstraintsforthe hp –LGRmethodis y = 266664 y D y C y E 377775 (4–118) andthetotalnumberofNLPconstraintsis Nn + Ns + q .Itisnotedthatunlikethe hp –LGpseudospectralmethod,thereisnodependenceonthenum berofvariablesor numberofconstraintsonthetotalnumberofmeshintervals. Therefore,asthenumber ofapproximatingmeshintervalsgrowslarge,itdoesnotaff ectthesizeoftheNLP denedbythe hp –LGRpseudospectralmethod.Finally,theconstraintJacob ianforthe hp –LGRpseudospectralmethodisgivenas G ( z )= @ y @ z (4–119) ThedensityoftheconstraintJacobianisnowexamined.Fort heconstraints associatedwith y C ,eachrowcorrespondstoaninequalitypathconstrainteval uatedata collocationpoint.ExaminingEq.( 4–83 ),ateachcollocationpoint, y C isonlyafunction ofthestateandcontrolatthatcollocationpointandtheini tialandnaltime.Therefore, therowsof y C aresparse.Furthermore, y E isonlyafunctionofthestateattheinitial andnaltime,andtheinitialandnaltime.Therefore,ther owsof y E aresparseaswell. ThedensityofanLGRpseudospectralmethodisaresultofthe rowsofthe constraintJacobiancorrespondingto y D .ItisseeninEq.( 4–82 )thatforanycollocation point,Eq.( 4–82 )isafunctionofthestateandcontrolatthatcollocationpo int,theinitial andnaltime, and N k additionalvariablescorrespondingtotheLagrangepolyno mial approximationtothestate.Therefore,as N k growslargeinanyapproximatinginterval, thenumberofnon-zeroelementsinanyrowof y D growslargeaswell. Figs. 4-5 and 4-6 showsaqualititiverepresentationoftheconstraintJacob ianfor the hp –LGRpseudospectralmethod.Fig. 4-5 representsan hp –LGRpseudospectral approximationwhere N k islargeineachapproximatingintervalandFig. 4-6 represents 84

PAGE 85

an hp –LGRpseudospectralapproximationwhere N k issmallineachapproximating interval.Itisseenthatforlarge N k ,large,denseblocksariseintheconstraintJacobian. Forsmallvaluesof N k ,thesedenseblocksaremuchsmaller.Therefore,thesparsi tyof an hp –LGRpseudospectralmethodisdependentupon N k ,or,thedegreeofpolynomial approximationwithinanyinterval.Forlow-degreeapproxi mations,asparseNLPis obtained.Forhigh-degreeapproximations,theNLPismuchm oredense. Points Figure4-5.QualitativeRepresentationofconstraintJaco bianforLGRpointswithLarge N k 85

PAGE 86

Time Figure4-6.QualitativeRepresentationofconstraintJaco bianforLGRpointswithSmall N k 4.4Examples The hp –LGandLGRpseudospectralmethodsarenowappliedtotwoopt imal controlproblemswithanalyticsolutions.Intherstexamp le,theoptimalstate,control andcostatearesmooth,whileinthesecondexample,thestat e,controlandthesecond componentofthecostatearenon-smoothandcontinuous,and therstcomponentof thecostateisdiscontinuous.Thepurposeoftheexamplesis toassesstheaccuracyof themultiple-intervalLGandLGRmethodsforapproximating state,controlandcostate. 86

PAGE 87

Theproblemsaresolvedbyeither p or h –methods.For p –methodproblems,thenumber ofintervalsisxedandthedegreeofapproximationwithine achintervalisallowedto vary.For h –methodproblems,thedegreeofapproximationisxedwithi neachinterval andthenumberofintervalsisallowedtovary.Finally,thet erminology“ h – x ”denotesan h –methodwith x collocationpointsineverymeshinterval. 4.4.1Example1-ProblemwithaSmoothSolution Considerthefollowingoptimalcontrolproblem[ 39 ].Minimizethecostfunctional J = 1 2 Z t f 0 ( y + u 2 ) dt (4–120) subjecttothedynamicconstraint y =2 y +2 u p y (4–121) andtheboundaryconditions y (0)=2 (4–122) y ( t f )=1 (4–123) where t f =5 .Theexactsolutiontotheoptimalcontrolproblemisgivena s y ( t )= x 2 ( t ) (4–124) y ( t )= x 2 p y (4–125) u ( t )= 2 ( t ) p y ( t ) (4–126) where x ( t ) and x ( t ) aregivenas 264 x ( t ) x ( t ) 375 = exp ( A t ) 264 x 0 x 0 375 (4–127) 87

PAGE 88

where A = 264 1 1 1 1 375 (4–128) and x 0 = p 2 (4–129) x f =1 (4–130) x 0 = x f B 11 x 0 B 12 (4–131) B = 264 B 11 B 12 B 21 B 22 375 = exp ( A t f ). (4–132) Theconvergenceratesfor p anduniformlydistributed h – 2 h – 3 ,and h – 4 –LG andLGRmethodsareanalyzed.Furthermore,the p –methodusesonlyasingle approximatinginterval.Becausetheproblemhasasmoothsol ution,itisexpectedthata p –methodwilldemonstratethefastestconvergenceratesfor thisproblem.Furthermore, becausetheproblemhasacontinuouscostate,itisexpected thatthemultiple-interval LGandLGRmethodswillprovideaccurateapproximationstot hecontinuous-time optimalcontrolproblem. Fig. 4-7 showsthestate,costate,andcontrolsolutionforthisprob lem.Asisseen inFig. 4-7 ,thestate,costate,andcontrolareallsmooth.Figs. 4-8A 4-8B ,and 4-8C showthemaximumstate,control,andcostateerroratthecol locationpointsasa functionofthenumberofcollocationpoints,respectively ,forboththe p –LGmethod and h – x –LGmethods.Furthermore,themaximumstate,control,andc ostateerrorat thecollocationpointsasafunctionofthenumberofcolloca tionpointsfor p –LGRand h – x –LGRmethodsareshowninFigs 4-9A 4-9B ,and 4-9C .Becausethesolutionto thisproblemissmooth,forbothLGandLGRmethods,the p –methoddemonstratesthe mostrapidrateofconvergencetotheexactsolution.Furthe rmore,itisnotedthatthe p –methodsconvergeexponentially.Next,itisseenthatthe h –methodsconverge 88

PAGE 89

0 0 2 2 3 1 1 1 4 5 2 4 3 SolutionTime y ( t ) y ( t ) u ( t ) Error Figure4-7.SolutiontoExample 4.4.1 atasignicantlyslowerrateascomparedwiththe p –methods.Inthisexample, byeitherLGorLGRmethod,convergenceisachievedmuchmore rapidlyusinga p –method.Moreover,asthedegreeofpolynomialapproximati onisincreasedutilizing the h -methods,therateofconvergenceisincreasedaswell.Itis alsonotedthatthere isvirtuallynodifferencebetweentheconvergenceratesof theLGandLGRmethods forthestateandcontrol.Furthermore,becausethemultipl e-intervalLGRmethodis amorecomputationallyefcientNLPforagivennumberofcol locationpoints,ifonly anaccuratestateandcontrolaredesired,thisresultimpli esthattheLGRmethod shouldbeused.Interestingly,theerrorinthecostatefort he h –LGmethodsisnoticeably smallerthanthe h –LGRmethodsforthisproblem.Thisresultimpliesthatifan accurate costateapproximationisrequiredbyamultiple-intervalL GorLGRmethod,theLG methodshouldbeused. 89

PAGE 90

0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 AStateError. 0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 BControlError. 0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 CCostateError. Figure4-8.ConvergenceRatesUtilizing p h – 2 h – 3 h – 4 –LGMethodsforExample 4.4.1 4.4.2Example2-Bryson-DenhamProblem Considertheproblemtominimizethecostfunctional[ 2 ] J = 1 2 Z 1 0 a 2 dt (4–133) subjecttothedynamicequations x = v (4–134) v = a (4–135) 90

PAGE 91

0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 AStateError. 0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 BControlError. 0 0 2 2 4 6 8 100 50150 200 log 10 MaximumError NumberofCollocationPoints p h – 2 h – 3 h – 4 CCostateError. Figure4-9.ConvergenceRatesUtilizing p h – 2 h – 3 h – 4 –LGRMethodsforExample 4.4.1 theboundaryconditions x (0)= x (1)=0 (4–136) v (0)= v (1)=1 (4–137) andtheconstraint x ( t ) l 91

PAGE 92

Thesolutiontothisproblemisgivenas x = 8>>>><>>>>: l [1 (1 t 3 l ) 3 ]:0 t 3 l l :3 l t 1 3 l l [1 (1 1 t 3 l ) 3 ]:1 3 l t 1 (4–138) v = 8>>>><>>>>: (1 t 3 l ) 2 :0 t 3 l 0:3 l t 1 3 l (1 1 t 3 l ) 2 :1 3 l t 1 (4–139) a = 8>>>><>>>>: 2 3 l (1 t 3 l ):0 t 3 l 0:3 l t 1 3 l 2 3 l (1 1 t 3 l ):1 3 l t 1 (4–140) x = 8>>>><>>>>: 2 9 l 2 :0 t 3 l 0:3 l t 1 3 l 2 9 l 2 :1 3 l t 1 (4–141) v = 8>>>><>>>>: 2 3 l (1 t 3 l ):0 t 3 l 0:3 l t 1 3 l 2 3 l (1 1 t 3 l ):1 3 l t 1 (4–142) Unlikethepreviousexamplethathadasmoothsolutionandac ontinuouscostate, thesolutiontothisproblemisonlypiecewisesmoothandhas adiscontinuouscostate solutionfor x .Withinanypiecewiseinterval,thehighestorderbehaviori ntheexact solutionisorderthree.Therefore,theonlydifcultyinac curatelyapproximatingthe solutiontothisproblemisinapproximatingthenonsmoothf eaturesat t =3 l and 92

PAGE 93

t =1 3 l andthefactthat x isdiscontinuousat t =3 l and t =1 3 l .Inthisexample, l =1 = 12 4.4.2.1Case1:LocationofDiscontinuityUnknown p –LGandLGRcollocationmethodsanduniformlydistributed h – 2 –LGandLGR collocationmethodsareusedtoapproximatethesolutionto thisprobleminthecase thatthelocationofdiscontinuityin x isassumed not tobeknown.Inthiscase,mesh pointswill,ingeneral,notbeatthelocationofthediscont inuityin x .Figs. 4-10A – 4-10D andFigs. 4-11A – 4-11D showtheconvergenceratesforstate,control,andcostatef or thisproblemasafunctionofcollocationpointsfor p and h – 2 –LGmethodsand p and h – 2 –LGRmethods,respectively.First,itisseenthatthereise ssentiallynodifference intheconvergenceratesofthe p and h – 2 –methods.Becausethisproblemhasa non-smoothsolution,exponentialconvergenceratesfroma p –methoddonotexist. Furthermore,itisseenthatthereisnosubstantialdiffere nceinaccuracybetween LGandLGRmethods.Itisalsonotedthatthestate,control,a nd v doconverge totheexactsolutionasthenumberofcollocationpointsare increased.Itisalso notedthat x doesnotconvergeasthenumberofcollocationpointsareinc reased. Figs. 4-12A – 4-12E showtheapproximationforthe h – 2 –methodutilizing150collocation points.Visually,theproblemisreasonablyapproximatedex ceptfor x .Thedifculty inapproximating x isisolatedtotheregionsnearthediscontinuities(i.e.,t hepointat approximately t =0.25 ). 4.4.2.2Case2:LocationofDiscontinuityKnown Next,thesolutiontotheproblemisapproximatedknowingth eexactlocationsof discontinuitiesin x p and h – 3 –LGandLGRmethodsareusedtoapproximatethe solutiontothisproblem.First,denethetimedomains T 1 2 [0,3 l ] T 2 2 [3 l ,1 3 l ] and T 3 2 [1 3 l ,1] .Forthe p –methodapproximations,asingleapproximating intervalisusedforeachof T 1 T 2 ,and T 3 withavariabledegreepolynomialineach interval.Furthermore,forthe h – 3 –method,avariablenumberofuniformlydistributed 93

PAGE 94

0 0 1 2 4 100 50150 3 log 10 MaximumErrorNumberofCollocationPoints p h – 2 t ) AMaximumErrorintheState. -1.5 0 0 1 0.5 2 100 50150 log 10 MaximumErrorNumberofCollocationPoints p h – 2 BMaximumErrorintheControl. 0 1.5 0.5 1 100 50150 p log 10 MaximumErrorNumberofCollocationPoints p h – 2 CMaximumErrorin x -1.5 0 0 1 0.5 0.5 1 2 100 50150 log 10 MaximumErrorNumberofCollocationPoints p h – 2 DMaximumErrorin v Figure4-10.Case1:Convergencefor p –andUniformlySpaced h – 2 –LGmethodsfor Example 4.4.2 approximatingintervalsofthirddegreeareusedin T 1 and T 3 ,while T 2 isapproximated usingasingleintervalofthirddegree.Figs. 4-13A – 4-13F showthemaximumerrorsat thecollocationpointswhentheapproximationtothesoluti onoftheproblemisdivided inthemannerdescribed.InFigs. 4-13A 4-13C ,and 4-13E ,thex-axisdenotesthe numberofintervalsin T 1 and T 3 .InFig. 4-13B 4-13D ,and 4-13F ,thex-axisdenotes thedegreeofapproximatingpolynomialin T 1 T 2 ,and T 3 .First,itisnotedthatthe maximumdegreeoftheexactsolutionisthree.Therefore,th e p and h – 3 –methods shouldexactlyapproximatethesolutionofthisproblem.Any deviationfromanexact 94

PAGE 95

0 0 1 2 4 100 50150 3 log 10 MaximumErrorNumberofCollocationPoints p h – 2 t ) AMaximumErrorintheState. -1.5 0 0 1 0.5 0.5 1 2 100 50150 log 10 MaximumErrorNumberofCollocationPoints p h – 2 BMaximumErrorintheControl. 1.4 1.6 1.2 0 1 100 50150 0.8 p log 10 MaximumErrorNumberofCollocationPoints p h – 2 CMaximumErrorin x -1.5 0 0 1 0.5 0.5 1 2 100 50150 log 10 MaximumErrorNumberofCollocationPoints p h – 2 DMaximumErrorin v Figure4-11.Case1:Convergencefor p –andUniformlySpaced h – 2 –LGRMethodsfor Example 4.4.2 approximationisduetothefactthatthetransformedadjoin tsystemsareaninexact discreterepresentationofthecontinuous-timerst-orde roptimalityconditionsifmesh pointsareatthelocationsofdiscontinuityinthecostate. Itisseenthatbyeither p or h – 3 –LGmethodapproximations,theerrorforthe state,controlandcostateareallnearmachineprecisionfo reithera low numberof intervalsordegreeofapproximatingpolynomial.Asthenumb erofintervalsordegreeof polynomialisincreased,theapproximationaccuracybecom es worse .Furthermore,by p or h – 3 –LGRmethods,theerrorasafunctionofnumberofapproximat ingintervalsor 95

PAGE 96

0 0 1 0.20.4 0.6 0.8 SolutionTime x ( t ) x ( t ) 0.02 0.04 0.06 0.08 0.1 A x ( t ) and x ( t ) 0 0 1 0.5 0.5 1 1 0.20.4 0.6 0.8 SolutionTime v ( t ) v ( t ) B v ( t ) and v ( t ) 0 0 2 1 2 4 6 8 10 0.2 0.4 0.6 0.8 SolutionTime t ) a ( t ) a ( t ) C a ( t ) and a ( t ) Time 0 0 1 40 60 0.2 0.4 0.6 0.8 20 x ( t ) x ( t ) x ( t ) 20 40 D x ( t ) and x ( t ) 0 0 2 1 4 2 6 8 10 0.20.4 0.6 0.8 SolutionTime v ( t ) v ( t ) E v ( t ) and v ( t ) Figure4-12.Case1:SolutionbyaUniform h – 2 -MethodUsing150CollocationPoints ComparedwiththeExactSolution. 96

PAGE 97

degreeofapproximatingpolynomialneitherconvergesnord iverges.Rather,errorsare bandedbetween O (10 4 ) and O (10 16 ) .Whilehighlevelsofaccuracyaredemonstrated by p and h – 3 –LGandLGRmethods,thebehaviorsofthesemethodswithresp ectto thiscaseareerraticbecausethetransformedadjointsyste misaninexactdiscrete representationofthecontinuous-timerst-orderoptimal ityconditions. 4.5Summary The hp –LGandLGRpseudospectralmethodsforapproximatinganonl inear continuous-timeoptimalcontrolproblemhavebeenposed.L GandLGRpointsare usedbecausetheyprovidehighlyaccuratequadratureappro ximationsandavoidRunge phenomenonintheLagrangepolynomialapproximationofthe state.Furthermore,the hp –LGandLGRmethodshavebeenposedsuchthatanarbitrarynum berofintervals, witharbitrarylocation,andanarbitrarynumberofcolloca tionpointsperintervalmaybe used.Overregionswherethesolutionissmooth,increasing thenumberofcollocation pointsinapproximatingintervalswillresultinefcienth igh-accuracyapproximations. Forregionsofthetrajectorywherethesolutionisnon-smoo th,increasingthenumberof low-orderapproximatingintervalswillresultincomputat ionallyefcientapproximations ofthesenon-smoothregions.Next,theKKTconditionsofthe hp –LGandLGRmethods arederived.Thetransformedadjointsystemsarederivedfr omtheKKTconditions andareshowntobediscreterepresentationsofthecontinuo us-timerst-order necessaryconditionsifthecostateiscontinuous.Iftheco stateisdiscontinuous,then thetransformedadjointsystemsarenotdiscreterepresent ationsofthecontinuous-time rst-ordernecessaryconditions.Twoexamplesareusedtod emonstratetheaccuracy ofapproximationofthe hp –LGandLGRmethods.Thesolutiontotherstexample issmoothandhasacontinuouscostatewhilethesolutiontot hesecondexampleis nonsmoothandhasadiscontinuouscostate.Itisseeninthe rstexamplethatas thenumberofapproximatingintervalsordegreeofapproxim ationinanyintervalis increased,theaccuracyinapproximatingthestate,contro l,andcostateisincreasedas 97

PAGE 98

-18 0 6 8 10 12 14 16 10 40 50 20 30 LGPts log 10 MaximumError t ) NumberofIntervals LGRPts AStateErrorbyLGandLGR h – 3 –Methods. -18 0 6 8 10 12 14 16 10 40 50 20 30 LGPts log 10 MaximumErrorNumberofCollocationPoints LGRPts BStateErrorbyLGandLGR p –Methods. 0 4 6 8 10 12 14 16 10 40 50 20 30 LGPts p log 10 MaximumError NumberofIntervals LGRPts CControlErrorbyLGandLGR h – 3 –Methods. 0 4 6 8 10 12 14 16 10 40 50 20 30 LGPts log 10 MaximumErrorNumberofCollocationPoints LGRPts DControlErrorbyLGandLGR p –Methods. 0 4 6 8 10 12 14 10 40 50 20 30 LGPts oints log 10 MaximumError NumberofIntervals LGRPts ECostateErrorbyLGandLGR h – 3 –Methods. 0 2 4 6 8 10 12 14 10 40 50 20 30 LGPts log 10 MaximumErrorNumberofCollocationPoints LGRPts FCostateErrorbyLGandLGR p –Methods. Figure4-13.Case2:Convergencefor p and h – 3 –LGandLGRMethodsforExample 4.4.2 98

PAGE 99

well.Furthermore,thereislittledifferenceintheapprox imationaccuracyofthestate andcontrolusingLGorLGRpoints,but,thecostateapproxim ationissignicantlymore accurateusingLGpoints.Next,thesecondexamplestudiesa problemwhosesolution hasadiscontinuouscostate.Itisdemonstratedthatifmesh pointsarenotlocatedat thelocationofdiscontinuityinthecostate,signicanter rorsexistinapproximatingthe discontinuity.Furthermore,ifmeshpointsarelocatedatt helocationsofdiscontinuity inthecostate,highlevelsofaccuracyareobtainedintheap proximationalthoughthe trendsarequiteerratic.FortheLGmethods,afewnumberofa pproximatingintervals orlow-degreeapproximationsbestapproximatedthesoluti ontotheproblem.Forthe LGRmethods,theerrorsarebandedwithinalargerangefrom O (10 4 ) to O (10 16 ) Therefore,ifthecostateisdiscontinuous,itisbesttoput meshpointsattheexact locationofdiscontinuity,but,guaranteeinghigh-accura cybytheproposed hp –LGor LGRmethodsisdifcultforthecaseofdiscontinuouscostat esolutions.Finally,the sparsityofthe hp –LGandLGRmethodsisstudied.Itisseenthatasthenumberof collocationpointswithinanintervalgrows,sodoestheden sityoftheapproximatingNLP. Therefore,thereisatradeoffthatexistsbetweenincreasi ngthenumberofcollocation pointsinanintervalinordertoobtainhigheraccuracysolu tionsthatalsoincreasesthe densityoftheNLP.Ifanintervalisapproximatingasmoothr egionofasolution,the rapidratesofconvergencebyincreasingthenumberofLGorL GRpointsfaroutweigh theincreaseinNLPdensity.Ifanintervalisapproximating anon-smoothregionofa solution,thenitismoreefcienttoincreasetheaccuracyo ftheapproximatingNLP bydividingtheproblemintomorelow-order,sparse,approx imatingintervals.Finally, itisseenthatthesizeofthe hp –LGmethodisdependentuponboththenumberof collocationpointsandthenumberofmeshintervalswhereas thesizeofthe hp –LGR methodisonlydependentuponthenumberofcollocationpoin ts.Therefore,ifalarge numberofmeshintervalsaretobeused,the hp –LGRmethodresultsinamoreefcient NLP. 99

PAGE 100

CHAPTER5 HP –MESHREFINEMENTALGORITHMS Inthischapter,two hp -meshrenementalgorithmsarepresentedusingthe hp –LGandLGRmethodsofChapter 4 .Ineachalgorithm,anaiveinitialmeshis selected.Becauseitisassumedthatnoaprioriknowledgeoft hesolutionstructure ofacontinuous-timeoptimalcontrolproblembeingapproxi matedisknown,thenaive initialmeshesareeitherasingle-intervalofmoderatedeg reeorrelativelyfewuniformly distributedintervalsoflow-degree.Fromthisnaiveappro ximation,bothalgorithms assesstheapproximationerrorsineachinterval.Iftheacc uracyoftheapproximation inanintervalneedstobeimproved,eitherthedegreeofappr oximationintheinterval isincreasedortheintervalissubdividedintomoreapproxi matingsubintervals.These strategiesarethenrepeatediterativelyuntilanaccurate approximationisobtained. Becauseaknownexactsolutionisnotgenerallyavailable,so lutionerrorsinan intervalareapproximatedbyexamininghowcloselythediff erential-algebraicconstraint equationsaresatisedbetweencollocationpoints. Therstmeshrenementalgorithmismoreheavilyweighteda sa p –biased hp –meshrenementalgorithm.Anintervalissubdividedonlyif itisbelievedthat exponentialconvergencedoesnotholdinthatinterval.The secondalgorithmisan h –biased hp –meshrenementalgorithm.Inthissecondalgorithm,subdi visionofan intervalischosenmuchmorefrequentlyandanincreaseinde greeofpolynomial approximationischosenonlyifthesolutioninanintervali sconsideredsmooth.First, adiscussiononassessmentoftheerrorsinanapproximating gridarediscussed. Next, hp –meshrenementalgorithm1and2aredescribed.Itisnotedt hat hp –mesh renementalgorithm1and2areadaptationsofthealgorithm spresentedinRefs.[ 66 ] and[ 67 ]. 100

PAGE 101

5.1ApproximationofErrorsinaMeshInterval Consideran hp –LGorLGRapproximationtothesolutionofacontinuous-tim e optimalcontrolproblemasdenedinChapter 4 with K meshintervals.Furthermore, let t 0 < t 1 < ... < t K bethemeshpointssuchthat t k 1 denesthebeginningof the k th meshand t k denestheendofthe k th mesh.Next,let N k bethenumberof collocationpointsinmeshinterval k .Indesigningthe hp –meshrenementalgorithms anassessmentoftheapproximationaccuracyofeachapproxi matinginterval k must bemade without knowledgeoftheexactsolutiontoacontinuous-timeoptima lcontrol problem. Inordertoprovideanassessmentontheaccuracyoftheappro ximationwithout knowledgeoftheexactsolution,thesatisfactionofthedif ferential-algebraicconstraints betweencollocationpointsareexamined.Atthecollocation points,thedifferential-algebraic equationsareenforced.However,ifanapproximating hp –pseudospectralapproximation isaccuratelyapproximatingthe continuous-time optimalcontrolproblem,thenthe differential-algebraicconstraintsshouldbenearlysati sedatanarbitrarypoint.In ordertoevaluatethedifferential-algebraicconstraints betweencollocationpoints,the stateandcontrolmustbeinterpolatedtoanarbitrarypoint .Thestateisdenedatan arbitrarypointinaninterval k bytheLagrangepolynomialofEqs.( 4–24 )or( 4–73 ) forthe hp –LGorLGRmethods,respectively.Thecontrolontheotherha nd,isonly calculatedatthecollocationpointsanddoesnothaveauniq ueunderlyingpolynomial approximation.Therefore,anycontrolapproximationthat matchesthevalueofthe controlatthecollocationpointsisavalidapproximation( e.g.,linear,cubic,spline, Lagrangepolynomial).Inthisresearch,aLagrangepolynom ialapproximationof thecontrolisused.Ifthestateisbeingapproximatedbyahi gh-accuracyLagrange polynomial,thenitiscounter-productivetousealow-orde r,lessaccurateapproximation tothecontrol.Therefore,ineitherthe hp –LGorLGRmethod,theLagrangepolynomial approximationofthecontrolwithsupportpointsatthe N k collocationpointsininterval k 101

PAGE 102

isgivenas U ( k ) ( )= N k X i =1 U ( k ) i L ( k ) i ( ). (5–1) Usingtheaforementionedstateandcontrolapproximations ,thedifferential-algebraic constraintscannowbeevaluatedatanarbitrarypointiname shinterval k .Inorderto samplethedifferential-algebraicconstraintsbetweenth ecollocationpoints,dene L pointsininterval k as t ( k ) =( t ( k ) 1 ,..., t ( k ) L ) 2 [ t k 1 t k ] .Thedifferential-algebraic constraintscanthenbeevaluatedatthese L pointsas X ( k ) i ( t ( k ) l ) f ( k ) i ( X ( k ) l U ( k ) l ) = a ( k ) li (5–2) C ( k ) j ( X ( k ) l U ( k ) l )= b ( k ) lj (5–3) where i =1,..., n j =1,..., s ,and l =1,..., L .Next,dene e ( k ) max asthemaximum elementof a ( k ) li or b ( k ) lj .If e ( k ) max islessthanaspeciederrortolerance, ,thenthesolution intheintervalisconsideredaccurate.ItisnotedthatEq.( 5–2 )isthemagnitudeof theviolationinthecollocationequationswhileEq.( 5–3 )issimplyanevaluationofthe inequalitypathconstraintsat L points.Iftheinequalitypathconstraintsareinactive, thelefthandsideofEq.( 5–3 )willbenegative,andtherefore,belessthan .Whenthe pathconstraintsareactive,betweencollocationpointsth elefthandsideofEq.( 5–3 ) maybecomepositivethusviolatingthepathconstraint.Ifa nyelement a ( k ) li or b ( k ) lj is greaterthan ,thentheintervalisconsideredinaccurateandmustberen edsuch thatamoreaccurateapproximationismade.Thetwo hp –meshrenementalgorithms arenowpresentedthatdeterminehowtoreneinaccurateapp roximatingintervals.In hp –meshrenementalgorithm1,theerrorsestimatedinEqs.( 5–2 )–( 5–3 )areused directlytodetermineifaninaccurateintervalshouldbesu bdividedorifthedegreeof approximationshouldbeincreased.In hp –meshrenementalgorithm2,an h –biased strategyisused.Thecurvatureofthestateinanintervalis usedtodeterminethe placementofthemeshpointsandwhethertosubdivideaninte rvalorincreasethe degreeoftheapproximatingpolynomialwithinaninterval. Evenwhenitischosento 102

PAGE 103

increasethedegreeoftheapproximatingpolynomial,thede greeoftheapproximation remainssmall. 5.2 hp –MeshRenementAlgorithm1 Adescriptionof hp –meshrenementalgorithm1isnowgiven.Itisnotedthat hp –meshrenementalgorithm1issimilartothealgorithmpres entedinRef.[ 66 ]except forthefollowingdifferences.InRef.[ 66 ],acubicinterpolationtothecontrolwasused. Inthisdissertation,aLagrangepolynomialapproximation tothecontrolisusedwithin eachinterval.Furthermore,inRef.[ 66 ],errorsareassessedatthemidpointsbetween collocationpoints,whereasinthisdissertation,errorsa reassessedatthearbitrary points t ( k ) withineachinterval.Furthermore,Ref.[ 66 ]onlymodiedthemeshbased uponanalysisofthecollocationconditions.Inthisdisser tation,thecollocationconditions andpathconstraintsareusedtomodifythemesh.5.2.1IncreasingDegreeofApproximationorSubdividing Thefundamentalanddistinctpropertyofan hp –meshrenementalgorithmisthe stepthatdetermineswhethertoincreasethedegreeofpolyn omialapproximationinan inaccurateintervalortosubdividetheinterval.Consider againthematricesdeningthe violationofthedifferential-algebraicconstraintequat ionsinEqs.( 5–2 )and( 5–3 ) A ( k ) = a ( k ) li (5–4) B ( k ) = b ( k ) lj (5–5) wherethecolumn i of A ( k ) andcolumn j of B ( k ) denesthemagnitudeofviolation inthecollocationequationsforthe i th componentofthestateandviolationofthe j th pathconstraintequation,respectively,at t ( k ) .Letthestate X ( k ) m ( t ) m 2 [1, n ] bethe componentofthestateapproximationwiththemaximumvalue of a ( k ) li .Furthermore, lettheconstraint C ( k ) p p 2 [1, s ] bethe p th inequalitypathconstraintwiththemaximum valueof b ( k ) lj .Thecomponentofthestateorinequalitypathconstraintwi thmaximum errorisalwaysusedtodeterminehowtoproceedwiththemesh renement.If e ( k ) max is 103

PAGE 104

associatedwiththeinequalitypathconstraints,theinter valissubdivided.Itisassumed that,ifapathconstraintbecomesactiveinaninterval,the trajectoryisnon-smooth and,therefore,itisnecessarytosubdividetheinterval.I fthemaximumerroroccurs inthecollocationequations,thefollowingprocedureisus edtodeterminewhetherto increasethedegreeofpolynomialapproximationorsubdivi detheinterval.First,dene thecolumnof A ( k ) thatcontains e ( k ) max as e ( k ) = 266664 e ( k ) 1 ... e ( k ) L 377775 = 266664 a ( k ) 1 i ... a ( k ) Li 377775 (5–6) Next,denethearithmeticmeanof e ( k ) as e ( k ) = P Li =1 e i L (5–7) Finally,dene as = 266664 ( k ) 1 ... ( k ) L 377775 = 266664 e ( k ) 1 = e ( k ) ... e ( k ) L = e ( k ) 377775 (5–8) thendenestheratiooftheerroratapoint e ( k ) 1 comparedtotheaverageerror e ( k ) .Next,deneaparameter ,where isusedtodeterminewhethertoincrease thedegreeoftheapproximationintheintervalortosubdivi detheinterval.Ifevery elementof islessthan ,thennoparticularpointdominatestheerrorsininterval k Inthiscase,thedegreeofthepolynomialapproximationint heintervalisincreased. Furthermore,ifvaluesof aregreaterthan ,thenerrorsintheintervalaredominated atparticularpoints.Inthiscase,theproblemissubdivide datthesepointsofmaximum error. 104

PAGE 105

5.2.2LocationofIntervalsorIncreaseofCollocationPoin ts Nowthatadeterminationhasbeenmadeastowhetheraninaccu rateinterval k is tobemodiedbyincreasingthedegreeofthepolynomialorsu bdivisionoftheinterval, thedetailsofeachchoicearenowdiscussed.Considertheca sewhere e ( k ) max occursin thecollocationequationsandeveryelementof islessthan .Inthiscase,thenumber ofcollocationpointsinanintervalisincreasedbyaspeci edamount N (+) ,thatis, N ( i +1) k = N ( i ) k + N (+) (5–9) where N ( i ) k isthenumberofcollocationpointsininterval k ongriditeration i .Next, considerthecasewheresomeelementsof aregreaterthan or e ( k ) max occursin thepathconstraintequations.If e ( k ) max isinthecollocationequations,theintervalis subdividedateverylocalmaximumof thatisalsogreaterthan .Fig. 5-1 showshow todivideanintervalforthiscase.Furthermore,thenumber ofcollocationpointsineach newintervalissetto N (0) .If e ( k ) max occursinthepathconstraintequations,theinterval issubdividedatthelocationofmaximumviolationofthepat hconstraintandeachnew intervalhas N (0) collocationpoints. Figure5-1.LocationofIntervalsifSubdivisionisChosen. 105

PAGE 106

5.2.3StoppingCriteria Thealgorithmstopswhenthedynamicconstraintsandinequa litypathconstraints aresatisedtothetolerance ineveryintervalat L pointsbetweenthecollocation points.5.2.4 hp –MeshRenementAlgorithm1 hp –MeshRenementAlgorithm1isnowstated. 1.Ifthisistherstiteration,initializetheproblemchoo sing N (0) collocationpointsand K =1 numberofintervals,where N (0) ischosenbytheuser.Otherwise,solvethe NLPusingthecurrentmesh.Begin:For k =1,..., K 2.If e ( k ) max ,thendonotmodifyinterval k .Otherwise,continuetostep3. 3.Ifthemaximumerroroccursinthecollocationequationsa ndeveryelement of ( k ) < ,thenincreasethedegreeofapproximationofmeshinterval k by N (+) 4.Otherwise,iferrorsarelargestinthecollocationequat ions,renetheinterval atthepointsdeningthemaximumerrorsin .Iferrorsarelargestinthe pathconstraintequations,renetheintervalatthepointo fmaximumerror. Furthermore,setthenumberofcollocationpointsto N (0) =5 ineachnew interval. End:For k =1,..., K 5.ReturntoStep1. Severalimportantaspectsof hp –meshrenementalgorithm1arenowclaried.It isnotedthat servesasatuningparameterthatvariestheamountofweight ingfrom h to p –biases.Forverysmallvaluesof ,the h –biasisincreased.Forlargevaluesof ,thealgorithmismore p –biased. shouldbeselectedsuchthatdetectionoflarge isolatederrorsismadeandtheapproximationcanbesubdivi dedaboutthesepoints. Next,fortheproblemsanalyzedinthisdissertation, N (+) =10 becauseincreasing bymorethantencollocationpointsmaymaketheNLPunnecess arilydense.Next,all newlycreatedintervalsareinitializedwithvecollocati onpointsfortworeasons.First, ifexponentialconvergenceinanewlycreatedintervalhold s,vecollocationpointsmay 106

PAGE 107

providehigh-accuracysolutions.Second,initializingane wlycreatedintervalwithve collocationpointswillkeepthetotalnumberofcollocatio npointssmall. 5.3 hp –MeshRenementAlgorithm2 hp –meshrenementalgorithm2isnowdened.Differentfrom hp –meshrenement algorithm1, hp –meshrenementalgorithm2isdesignedtobeamuchmore h –biased algorithm.Itisnotedthat hp –meshrenementalgorithm2issimilartothealgorithm presentedinRef.[ 67 ]withtheonlymodicationinthisdissertationbeinghowth econtrol isinterpolated.Inthisdissertation,thecontrolininter val k 2 1,..., K isinterpolatedby aLagrangepolynomialusingsupportpointsatthe N k collocationpoints.InRef.[ 67 ], thecontrolisinterpolatedwithaLagrangepolynomialwith supportpointsatthe N k collocationpointsandoneadditionalnon-collocatedpoin tfor k 2 1,..., K 1 andbya Lagrangepolynomialwithsupportpointsatthe N K collocationpointsinthenalinterval K 5.3.1IncreasingDegreeofApproximationorSubdividing Ifthesolutioninameshinterval k isinaccurate,therstdeterminationiswhether thedegreeofapproximationintheintervalshouldbeincrea sedorwhethertheinterval shouldbesubdivided.Supposethat X ( k ) m ( ) isthecomponentofthestateapproximation inthe k th meshintervalthatcorrespondstothemaximumvalueof a ( k ) li .Furthermore,let thecurvatureofthe m th componentofthestateinmeshinterval k begivenby ( k ) ( )= j X ( k ) m ( ) j h 1+ X ( k ) m ( ) 2 i 3 = 2 (5–10) Next,dene ( k ) max and ( k ) bethemaximumandmeanvalueof ( k ) ( ) ,respectively. Furthermore,let r k betheratioofthemaximumtothemeancurvature: r k = ( k ) max ( k ) (5–11) 107

PAGE 108

If r k < r max ,where r max > 0 isauser-denedparameter,thenthecurvatureofmesh interval k isconsidereduniformandthedegreeofapproximationisinc reasedwithinthe interval.If r k > r max ,thenalargecurvaturerelativetotheremainderoftheinte rvalexists andtheintervalissubdivided.5.3.2IncreaseinPolynomialDegreeinInterval k If r k < r max ,thenthepolynomialdegreeinmeshinterval k needstobeincreased. Let M bethecurrentdegreeoftheapproximatingpolynomialinthe interval.Thedegree N k ofthepolynomialforthenextgriditerationinmeshinterva l k isgivenbytheformula N k = M + ceil (log 10 ( e ( k ) max = ))+ A (5–12) where A > 0 isanarbitraryintegerconstantand ceil istheoperatorthatroundstothe nexthighestinteger.5.3.3DeterminationofNumberandPlacementofNewMeshPoint s If r k > r max ,thenmeshinterval k issubdivided.Twoquantitiesarenowdened:The numberofsubdivisionsandthelocationofeachsubdivision .First,thenewnumberof meshintervals,denoted n k ,iscomputedas n k = ceil ( B log 10 ( e ( k ) max = )), (5–13) where B isanarbitraryintegerconstant.Thelocationsofthenewme shpointsare determinedusingtheintegralofacurvaturedensityfuncti oninamannersimilartothat giveninRef.[ 30 ].Specically,let ( ) bethedensityfunction[ 30 ]givenby ( )= c ( ) 1 = 3 (5–14) where c isaconstantchosensothat Z +1 1 ( ) d =1. 108

PAGE 109

Let F ( ) bethecumulativedistributionfunctiongivenby F ( )= Z 1 ( ) d (5–15) The n k newmeshpointsarechosensothat F ( i )= i 1 n k ,1 i n k +1. ItisshowninRef.[ 30 ]thatplacingthemeshpointsinthismannerisresultsinthe smallest L 1 normerrorwhenapiecewiselinearapproximationofacurvei sused.Finally, itisnotedthatif n k =1 ,thennosubintervalsarecreated.Therefore,theminimumv alue for n k shouldbeatleasttwo.Figs. 5-2A and 5-2B showqualitativelyhowtodividean intervalintomoreintervalsforaconstantcurvatureinani ntervalandanintervalwitha largecurvaturerelativetotheremainderoftheinterval,r espectively. 5.4 hp –MeshRenementAlgorithm2 hp –meshrenementalgorithm2isnowdened.Thealgorithmisi nitializedfora coarseuniformlydistributedmeshwithaxed,low-degreep olynomialineachinterval. Thisapproximationissolved,andtheneachmeshintervalis analyzed.Iferrorsare lessthan ,thentheintervalisnotmodied.Otherwise,iferrorsareg reaterthan thentheintervalismodiedaccordingtothediscussionofSe ctions 5.1 through 5.3.3 Therenementprocesscontinuesiterativelyuntilasuitab leapproximatingmeshis constructed. hp –meshrenementalgorithm2isnowgiven: 1.SolvetheNLPusingthecurrentmesh. Begin:For k =1,..., K 2.If e ( k ) max ,donotmodifyinterval k 3.Ifeither r k r max or N k > M ,thenrenethe k th meshintervalinto n k subintervalswhere n k isgivenby( 5–13 ).Setthedegreeofthepolynomials oneachofthesubintervalstobe M andproceedtothenext k 4.Otherwise,setthedegreeofthepolynomialsonthe k th subintervaltobe N k whichisgivenin( 5–12 ). End:For k =1,..., K 109

PAGE 110

AConstantCurvature. BVariableCurvature. Figure5-2.LocationofIntervalsbasedonCurvature.5.ReturntoStep1. Severalimportantaspectsof hp –meshrenementalgorithm2arenowclaried. hp –meshrenementalgorithm2isa h –biasedalgorithm.Itisnotedthatunlike hp –mesh renementalgorithm1,anincreaseinthedegreeofpolynomi alapproximationis onlyallowedoncewithinaninterval.Iftheerrorsinapprox imationarestilllargeafter increasingthedegreeofapproximation,theintervalisfor cedtosubdivideasopposed toincreasingthedegreeofapproximationagain.Furthermo re,Eq.( 5–12 )denesthe numberofcollocationpointsthatareaddedtoanintervalba sedona log 10 differencein 110

PAGE 111

thecurrenterroranddesirederrorplusanarbitraryconsta nt A .Ifthecontinuous-time optimalcontrolproblemisknowntohaveasmoothsolutionor ifcomputationalefciency oneachmeshiterationisnotalargeconcern, A canbechosentobelarge.Next, Eq.( 5–13 )denesthenumberofintervalscreatedifsubdivisionisch osen.Eq.( 5–13 ) isafunctionofthe log 10 differenceinthecurrenterrorandthedesirederrorandan arbitrarymultiplier B .Atrade-offexistsforlargeandsmallvaluesof B .Forsmallvalues of B ,fewerintervalswillbecreatedoneachmeshiteration,but thenumberofmesh iterationsmaygrowlarge.Forlargevaluesof B ,moreintervalswillbecreatedoneach meshiteration,andthetotalnumberofmeshiterationsmayb ereduced,butthenumber ofmeshintervalsmaygrowlarge. 111

PAGE 112

CHAPTER6 EXAMPLES The hp –meshrenementalgorithmsposedinChapter 5 arenowexaminedon avarietyofoptimalcontrolproblems.Therstexampleisth emotivatingexample ofSection 3.1 whosesolution,aswerecall,issmooth.Thesecondexamplei sa reusablelaunchvehicleentryproblemwhosesolutionappea rstobesmooth.Thethird exampleisaproblemwherethetimescaleofthedynamicsbeco mefastasthenaltime increases.ThefourthexampleisthemotivatingexampleofSe ction 3.2 whosesolution, aswerecall,isnotsmooth.Finally,thefthexampleisamin imumtime-to-climbof asupersonicaircraft.Thisnalexamplehastwoactivepath constraintsalongthe trajectory.TheobjectiveofthisChapteristoassesshowwe llthe hp –meshrenement algorithmsperformrelativeto p and h –methods. Theexamplesweresolvedwithmodicationsoftheopen-sour ceoptimalcontrol softwareGPOPS[ 37 ]usingtheNLPsolverSNOPT[ 68 ].Allcomputationswere performedusinga2.5GHZCore2DuoMacbookProrunningMacOS-X 10.5.8with MATLABR2009b.Furthermore,analyticderivativeswereused intheNLPsolver.The p –methoddiscretizationsweresolvedutilizingasingleapp roximatingintervalwitha variablenumberofcollocationpointsintheinterval.The p –methodsimplyassesses theviolationsofEqs.( 5–2 )–( 5–3 )at L uniformlydistributedpointsbetweencollocation pointswhere L =1000 .Iferrorsaregreaterthan ,thenumberofcollocationpoints isincreasedbyten.Next, hp –meshrenementalgorithm1( hp (1) )uses =3.5 and evaluatesineachintervaltheviolationsofEqs.( 5–2 )–( 5–3 )at L =100 uniformly distributedpointsbetweencollocationpoints. hp –meshrenementalgorithm2( hp (2) ) uses r max =2 A =1 B =2 ,and L =30 uniformlydistributedpointsbetweencollocation points.Finally,the h –methodisanadaptationof hp –meshrenementalgorithm2.For the h –method, hp –meshrenementalgorithm2isusedwith r max < 1 ,permittingonly segmentdivisionandpreventinggrowthinthedegreeofthea pproximatingpolynomial. 112

PAGE 113

Theterminology“ h – x –”denotesan h –methodwith x numberofcollocationpoints ineachintervalwhile“ hp (2) – x –”denotesthe hp (2) –methodwithaminimumof x collocationpointsineachinterval.Allexamplesweresolve dutilizingLGRpoints. Finally,amaximumlimitof200collocationpointsisgivenf orthe p –method.Ifasuitable approximationisnotobtainedfor200collocationpoints,t henitisdecidedthatan accurateapproximationcannotbeobtained.Furthermore,f orthe hp (1) hp (2) ,and h –methods,ameshiterationlimitof40isgiven.Ifasuitable approximationisnotfound within40meshiterations,thenitisdecidedthatanaccurat eapproximationcannotbe obtained. 6.1Example1-MotivatingExample1 ConsideragainthemotivatingexamplegivenbyEqs.( 3–1 )–( 3–3 ).Theproblemis nowrestated.Minimizethecostfunctional J = t f (6–1) subjecttothedynamicconstraints x =cos( ),_ y =sin( ), = u (6–2) andtheboundaryconditions x (0)=0, y (0)=0, (0)= x ( t f )=0, y ( t f )=0, ( t f )= (6–3) ThesolutiontotheoptimalcontrolproblemofEqs.( 6–1 )–( 6–3 )is ( x ( t ), y ( t ), ( t ), u ( t ))=( sin( t ), 1+cos( t ), t ,1) (6–4) where t f =2 .AswasdemonstratedinSection 3.1 ,a p –methodresultsinmost efcientapproximationtothisproblem.Therefore,itisex pectedthatthe p and hp (1) strategieswillperformthebest.Table 6-1 showstheinitialuniformlydistributedgridfor eachstrategystudiedandthenumberofcollocationpointsi neachinterval.Itisnoted 113

PAGE 114

Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 2 2 10 h – 2 2 10 Table6-1.InitialGridusedinExample 6.1 forEachStrategyUsingUniformlyDistributed Intervalswith N k CollocationPoints. thattheinitialgridforallofthemethodsresultsinsmallN LPsthataresolvedefciently. Next,Table 6-2 showsthemaximumabsoluteerrorevaluatedatthecollocati onpoints andthetotalCPUtimerequiredbytheNLPtosolvetheproblemf oreachvalueof Furthermore,thenumberofcollocationpoints,meshinterv als,andNLPdensityonthe nalgridareshown.First,itisnotedforanyvalueof ,themaximumabsoluteerror usinganystrategyislessthan .Next,itisnotedthatas isdecreased,themaximum absoluteerrorbyanymethoddecreases.Examiningthe p and hp (1) –methods,itis seenthatthe p and hp (1) –methodsproducetheexactsameapproximationgridsforthi s problem.Becausethesolutiontothisproblemissmooth,the hp (1) –methodchooses toincreasethedegreeoftheapproximatingpolynomialasop posedtosubdividingthe initialinterval.Moreover,whenexaminingtheresultsfor hp (2) – 2 ,again,itisnotedthat anincreaseinthedegreeoftheapproximatingpolynomialwa schosenineachinterval asopposedtosubdividingtheintervals.Thisisconsistent withthefactthatthesolution tothisproblemissmooth.Next,incomparingtheCPUtimesfor the hp (2) – 2 –method withthe hp (1) and p –methods,itisseenthatthesolutionstimesaresimilarfor anyvalue of .Eventhoughthe hp (2) – 2 –methodresultsinmorecollocationpointsthanthe hp (1) or p –methodsforanyaccuracylevel,theincreasedNLPsparsity ofthe hp (2) – 2 –method makesupfortheincreaseintotalnumberofcollocationpoin ts. Resultsaredrasticallydifferentforthe h – 2 –method.For =(10 2 ,10 4 ) ,theCPU timeforthe h – 2 –methodiscomparabletotheCPUtimeoftheotherthreemethod s.As theaccuracytoleranceistightened,however,the h – 2 –methodisseentobesignicantly lesscomputationallyefcientwhencomparedwiththeother threemethods.For = 114

PAGE 115

MaxAbs. Strategy CPU Colloc Mesh Mesh NLP Error Time(s) Pts Interval Iterations Density 10 2 1.51 10 5 p 0.0084 10 1 1 33.47 1.51 10 5 hp (1) 0.0084 10 1 1 33.47 5.96 10 3 hp (2) – 2 0.0093 20 10 1 9.1 5.96 10 3 h – 2 0.0093 20 10 1 9.1 10 4 1.51 10 5 p 0.0084 10 1 1 33.47 1.51 10 5 hp (1) 0.0084 10 1 1 33.47 1.21 10 7 hp (2) – 2 0.025 50 10 2 5.3 1.39 10 5 h – 2 0.087 160 80 3 1.2 10 6 4.89 10 15 p 0.027 20 1 2 29.64 4.89 10 15 hp (1) 0.027 20 1 2 29.64 4.88 10 11 hp (2) – 2 0.034 70 10 2 4.5 1.98 10 8 h – 2 5.36 1568 784 5 0.13 10 7 4.89 10 15 p 0.027 20 1 2 29.64 4.89 10 15 hp (1) 0.027 20 1 2 29.64 8.51 10 13 hp (2) – 2 0.037 80 10 2 4.26 3.09 10 10 h – 2 104.60 5760 2880 5 0.035 Table6-2.SummaryofAccuracyandSpeedforExample 6.1 UsingVariousCollocation StrategiesandAccuracyTolerances. 10 7 ,the p hp (1) ,and hp (2) – 2 –methodsallsolvetheprobleminapproximately 0.03 seconds.Bycomparison,the h – 2 –methodrequires 104.60 secondstosolvetheproblem and5760collocationpointsareused..EventhoughtheNLPiso nly0.035%dense,the increasedsizemakesthe h – 2 –discretizedNLPcomputationallyinefcient. Itisdemonstratedonthisproblemwithasimple,smoothsolu tionthata p –method providesthemostcomputationallyefcientapproximation .Furthermore, hp (1) behaves exactlylikethe p –method.Next,itisseenthattheperformanceofthe h –biased hp (2) – 2 –methodiscomparabletothe p –method.Bymodestlyincreasingthedegree ofpolynomialwithineachapproximatinginterval,efcien t,sparseapproximationsare obtainedusingrelativelyfewcollocationpointsfor hp (2) .Finally,fortightaccuracy tolerances,the h – 2 –methodbecomesmuchlesscomputationallyefcientwhen comparedwiththeothermethods.The h – 2 –methodfor =10 7 requiresmorethan 3000timestheCPUtimeoftheothermethodspresentedinorder toobtainanaccurate approximation. 115

PAGE 116

6.2Example2-ReusableLaunchVehicleRe-entryTrajectory Considerthefollowingoptimalcontrolproblem[ 69 ]ofmaximizingthecrossrange duringtheatmosphericentryofareusablelaunchvehicle.M inimizethecostfunctional J = ( t f ) (6–5) subjecttothedynamicconstraints r = v sin r = v cos r sin r cos = v cos r cos r v = D m g sin r r = L cos mv g v v r cos r = L sin mv cos r + v cos r sin tan r (6–6) andtheboundaryconditions r (0)=79248+ R e m r ( t f )=24384+ R e m (0)=0 deg ( t f )= Free (0)=0 deg ( t f )= Free v (0)=7803 m/s v ( t f )=762 m/s r (0)= 1 deg r ( t f )= 5 deg (0)=90 deg ( t f )= Free (6–7) where r isthegeocentricradius, isthelongitude, isthelatitude, v isthespeed, r istheightpathangle, istheazimuthangle, istheangleofattack, isthebank angle,and R e =6371203.92 mistheradiusoftheEarth.Moredetailsforthisproblem, includingtheaerodynamicmodel,canbefoundinRef.[ 69 ].SimilartoExample 6.1 thisproblemhasasmoothsolution.Figs. 6-1A – 6-1F showsthestateapproximationfor =10 6 forthe hp (2) – 3 –method.Itisseenthatallsixcomponentsofthestateappea r 116

PAGE 117

tobesmooth.Therefore,itisexpectedthata p –biasedmethodwillprovidethemost computationallyefcientapproximationstothisproblem. Table 6-3 showstheinitialuniformgridandnumberofcollocationpoi ntsineach intervalforeachofthemethodsstudied.Next,Table 6-4 showstheresultsforeach methodforvariousvaluesof .For =10 4 ,allmethodsperformcomparablywitha slightadvantagetothe hp (2) – x –and h – 2 –methods.Itisnotedforthisproblem,unlike thepreviousproblem,the hp (1) –methoddoesnotincreasethedegreeofapproximation fromtheinitialinterval,butinstead,actuallysubdivide stheinitialintervalintothree subintervalsofvecollocationpointseachthatsatisfyth eaccuracycriterion =10 4 For =10 6 ,the hp (1) method,interestingly,performstheworstwithrespecttoC PU time.Itisnoted,however,that hp (1) alsorequiresthemostmeshiterationstoobtain asolution.Furthermore,itisnotedthatalthoughthe p –methodresultsinthefewest numberofcollocationpoints,itisnotthemostcomputation allyefcientmethodfor =10 6 .Thereductionincomputationalefciencyofthe p –methodisduetothefact thatthe p –methodNLPdensityforthisproblemissignicantlygreate rthaneitherthe hp (2) – x –or h – 2 –methodNLPdensities.Interestingly,the hp (2) – x and h – 2 –methodsare themostcomputationallyefcientfor =10 6 For =(10 7 ,10 8 ) ,the h – 2 –methodissignicantlyslowerthantheothermethods. Becauseofthelargeproblemsize,theresultingNLPbecomesc omputationally intractable.Incomparingthe p hp (1) and hp (2) – x –methods,itisseenagainthatthe hp (1) methodrequiresthegreatestcomputationaltimebecauseit requiresthemostnumber ofmeshiterations.For =10 8 ,the p –methodprovidesthefastestcomputational speeds.Forthisproblem,forveryhighlevelsofaccuracy,t he p –methodisthemost computationallyefcient.Next,incomparingthe hp (2) – x –methods,itisseenthatthe hp (2) – 3 and hp (2) – 4 –methodsperformtwiceasfastasthe hp (2) – 2 –method.Anincrease inoneortwocollocationpointsforhighlevelsofaccuracy, therefore,stronglyaffectthe CPUtimesbyreducingtheoverallsizeoftheproblemfor =10 8 .Finally,itisnotedfor 117

PAGE 118

0 40 60 80 50 20 30 10002000 500 70 15002500 t (sec)h (km) A h (km)vs.Time. 0 0 40 60 80 20 10002000 500 15002500 t (sec) (deg) B (deg)vs.Time. 0 0 10 40 20 30 1000 2000 4985 500 15002500 t (sec) (deg) C (deg)vs.Time. 0 0 2 4 6 8 1000 2000 500 15002500 t (sec) v (km/s) D v (km/s)vs.Time. 0 0 2 2 4 6 1000 2000 500 15002500 t (sec) r (deg) E r (deg)vs.Time. 0 0 40 60 80 100 20 1000 2000 500 15002500 t (sec) (deg)F (deg)vs.Time. Figure6-1.Statevs.TimeontheFinalMeshforExample 6.2 118

PAGE 119

Strategy N k NumberofIntervals p 20 1 hp (1) 20 1 hp (2) – 4 4 20 hp (2) – 3 3 15 hp (2) – 2 2 20 h – 2 2 20 Table6-3.InitialGridusedinExample 6.2 forEachStrategyUsingUniformlyDistributed Intervalswith N k CollocationPoints. Strategy CPUTime(s) CollocPts MeshInterval MeshIterations NLPDensity 10 4 p 7.03 40 1 3 15.12 hp (1) 6.04 15 3 2 11.48 hp (2) – 4 1.58 44 11 2 3.86 hp (2) – 3 1.93 30 10 1 3.51 hp (2) – 2 1.16 42 20 2 3.52 h – 2 1.10 42 21 2 3.46 10 6 p 13.90 70 1 6 14.02 hp (1) 32.05 145 11 10 2.31 hp (2) – 4 4.22 136 32 4 1.31 hp (2) – 3 3.91 119 30 3 1.49 hp (2) – 2 4.18 132 34 4 1.32 h – 2 5.34 256 128 4 0.58 10 7 p 45.37 80 1 7 13.82 hp (1) 102.10 165 15 13 1.92 hp (2) – 4 30.33 254 61 3 0.70 hp (2) – 3 36.41 224 53 4 0.81 hp (2) – 2 56.01 272 64 4 0.66 h – 2 424.05 774 387 5 0.19 10 8 p 64.97 90 1 8 13.69 hp (1) 201.78 210 18 13 1.56 hp (2) – 4 105.3 366 84 5 0.48 hp (2) – 3 111.2 479 112 5 0.38 hp (2) – 2 217.2 578 114 4 0.33 h – 2 1815.6 2480 1240 6 0.06 Table6-4.SummaryofAccuracyandSpeedforExample 6.2 UsingVariousCollocation StrategiesandAccuracyTolerances. =10 8 thattheCPUtimesusing hp (2) – 3 and hp (2) – 4 arecomparabletotheCPUtimes usingthe p –method. 119

PAGE 120

6.3Example3:Hyper-SensitiveProblem Considerthefollowing hyper-sensitive [ 70 – 73 ]optimalcontrolproblemadaptedfrom Ref.[ 71 ].Minimizethecostfunctional J = 1 2 Z t f 0 ( x 2 + u 2 ) dt (6–8) subjecttothedynamicconstraint x = x 3 + u (6–9) andtheboundaryconditions x (0)=1, x ( t f )=1 (6–10) where t f isxed.Thisproblemisthenalexampleproblemconsidered inthisChapter whosesolutionissmooth.Thisproblemhasa“take-off',”“c ruise,”and“landing”structure wherealloftheinterestingbehavioroccursinthe“take-of f”and“landing”segments (i.e.,thersttenandlasttentimeunitsofthesolution).F orfurtherdetails,Ref.[ 71 ] containsamorein-depthdiscussionofthisproblem.Thesol utiontothisproblemisnow approximatedforthreedifferentvaluesof t f 6.3.1Solutionsfor t f =50 TheproblemdenedbyEqs.( 6–8 )–( 6–10 )isnowstudiedfor t f =50 .Theinitial approximatinggridsaredenedinTable 6-5 andresultsforthevariousapproximation methodsaregiveninTable 6-6 .Theapproximationtothestateandcontrolofthis problemfor t f =50 isshownforthe hp (2) – 2 –approximationwith =10 4 inFigs. 6-2A and 6-2B .ItisseenthatinFigs. 6-2A and 6-2B ,theinitialmeshduringthe“cruise” segmentwasnotmodied.Theonlyrenementthatoccurredwa stomoreaccurately approximatethe“take-off”and“landing”segments.Next,i tisinterestingtonote,that whilethesolutiontothisproblemissmooth,the p –methodwasnotabletoobtaina solutionfor =10 6 .Furthermore,eachmethodhasasimilarCPUtimefor = (10 2 ,10 4 ) .Aswiththeprevioustwoexamples,forthetightestaccuracy tolerance, 120

PAGE 121

Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 4 10 4 hp (2) – 3 10 3 hp (2) – 2 10 2 h – 2 10 2 Table6-5.InitialGridusedinExample 6.3.1 forEachStrategyUsingUniformly DistributedIntervalswith N k CollocationPoints. Strategy CPUTime(s) CollocPts MeshInterval MeshIterations NLPDensity 10 2 p 0.11 30 1 3 50.69 hp (1) 0.25 15 3 3 24.42 hp (2) – 4 0.12 56 14 3 6.75 hp (2) – 3 0.14 35 11 4 9.52 hp (2) – 2 0.15 30 14 4 9.43 h – 2 0.24 34 17 6 8.06 10 4 p 0.33 60 1 6 50.38 hp (1) 1.01 55 7 7 12.37 hp (2) – 4 0.25 86 20 3 4.75 hp (2) – 3 0.31 74 20 4 5.19 hp (2) – 2 0.24 65 21 4 5.71 h – 2 0.69 162 81 6 1.82 10 6 p 7.74 200 1 20 50.12 hp (1) 2.15 100 10 9 7.98 hp (2) – 4 0.96 168 36 6 2.61 hp (2) – 3 1.26 147 31 6 3.08 hp (2) – 2 0.80 129 28 5 3.58 h – 2 74.15 1328 664 7 0.22 Table6-6.SummaryofAccuracyandSpeedforExample 6.3.1 UsingVarious CollocationStrategiesandAccuracyTolerances. the h – 2 –methodrequiressignicantlymorecomputationaltimeino rdertoobtainan accurateapproximationwhencomparedwiththe hp –methods.Furthermore, unlike theprevioustwoexamples,the p –methodisunabletoobtainasolutionforthetightest accuracytolerance.Forthisexample,inordertoobtainane fcienthigh-accuracy solution,an hp –methodis required 6.3.2Solutionsfor t f =1000 TheproblemdenedbyEqs.( 6–8 )–( 6–10 )isnowexaminedfor t f =1000 withthe initialapproximatinggridsdenedinTable 6-7 andresultsforthevariousapproximation methodsgiveninTable 6-8 .First,the p –methodisthemostcomputationallyinefcient 121

PAGE 122

0 0 1 10 40 50 20 30 0.2 0.2 0.4 0.6 0.8 1.2 xt A x vs. t 0 0 1.5 2 2.5 0.5 0.5 1 10 40 50 20 30 tu B u vs. t Figure6-2.StateandControlvs.TimeontheFinalMeshforExam ple 6.3.1 for hp (2) – 2 and =10 4 122

PAGE 123

Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 4 10 4 hp (2) – 3 10 3 hp (2) – 2 10 2 h – 2 10 2 Table6-7.InitialGridusedinExample 6.3.2 forEachStrategyUsingUniformly DistributedIntervalswith N k CollocationPoints. approximationforthelowestaccuracytolerance.Thenumbe rofcollocationpoints forthe p –methodapproximationislargerthanthe hp (1) hp (2) – x ,and h – 2 –methods. Furthermore,theNLPdensityofthe p –methodismuchgreaterthananyoftheother methods.Therefore,the p –methodisthemostcomputationallyinefcientmethod. Next,itisseenthatthe p –methodisunabletoobtainanapproximationfortheaccurac y tolerancesof =(10 4 ,10 6 ) Finally,instudyingthe h – 2 approximations,for =(10 2 ,10 4 ) ,the h – 2 –method performssimilarlywhencomparedwiththe hp –methods.Againitisnotedforthe tightestaccuracytolerance,the h – 2 –methodissignicantlymoreinefcient,requiring approximately170secondstosolvetheproblemusing9meshi terationsand664 approximatingintervalsonthenalmesh.The hp –methods,however,required approximatelyvesecondstosolvetheproblemfor =10 6 .For t f =1000 ,again, high-accuracyefcientsolutionsareonlyobtainedbythe hp –methods. 6.3.3Solutionsfor t f =5000 Finally,theproblemisanalyzedfor t f =5000 withtheinitialapproximatinggrids givenbyTable 6-9 andresultsgiveninTable 6-10 .Itisnotedthatforthe p –method approximation, L =10000 inordertoincreasetheresolutionoftheerrorassessment. Forthisproblem,the p –methodisunabletoobtainasolutionfor any valueof .Itis againseenthatfor =(10 2 ,10 4 ) ,the hp and h –methodsperformcomparably,and forthetightestaccuracytolerance,the h – 2 –methodagainismuchlesscomputationally efcient. 123

PAGE 124

Strategy CPUTime(s) CollocPts MeshInterval MeshIterations NLPDensity 10 2 p 8.56 90 1 9 50.26 hp (1) 0.71 30 6 5 13.49 hp (2) – 4 0.12 56 14 3 6.75 hp (2) – 3 0.27 56 18 6 6.06 hp (2) – 2 0.27 38 18 6 7.53 h – 2 0.37 46 23 8 6.10 10 4 p 9.53 200 1 20 50.12 hp (1) 1.97 70 10 8 9.05 hp (2) – 4 0.25 86 20 3 4.75 hp (2) – 3 1.17 141 43 7 2.57 hp (2) – 2 0.66 104 37 8 3.42 h – 2 1.03 184 92 8 1.60 10 6 p – – – – – hp (1) 4.05 125 15 9 5.84 hp (2) – 4 4.71 214 46 6 2.06 hp (2) – 3 5.93 246 59 7 1.73 hp (2) – 2 2.97 209 53 8 2.08 h – 2 169.55 1328 664 9 0.22 Table6-8.SummaryofAccuracyandSpeedforExample 6.3.2 UsingVarious CollocationStrategiesandAccuracyTolerances. Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 4 10 4 hp (2) – 3 10 3 hp (2) – 2 10 2 h – 2 10 2 Table6-9.InitialGridusedinExample 6.3.3 forEachStrategyUsingUniformly DistributedIntervalswith N k CollocationPoints. Figs. 6-3A and 6-3B showthestateandcontrolapproximationfor =10 4 using the hp (2) – 2 –methodandFigs. 6-4A and 6-4B showthestateapproximationfortherst 15andlast15timeunits.Asexpected,itisseenthatrenemen toccursonlyatthestart andterminusofthesolution.6.3.4SummaryofExample 6.3 InExample 6.3 ,severalinterestingfeaturesaredemonstratedforapprox imating thesolutionofthisproblemby p hp ,and h –methods.Asthenaltimeincreases, thetimescaleofthisproblemdeningthe“take-off”and“la nding”segmentsbecome increasinglyfastwithrespecttothetotaltrajectorytime .Consequently,as t f is 124

PAGE 125

0 0 1 xt 0.2 0.2 0.4 0.6 0.8 1.2 10002000 30004000 5000 A x vs. t 0 0 1.5 2 2.5 0.5 0.5 1 tu 100020003000 4000 5000 B u vs. t Figure6-3.StateandControlvs.TimeontheFinalMeshforExam ple 6.3.3 for hp (2) – 2 and =10 4 125

PAGE 126

0 0 1 5 10 15 xt 0.2 0.4 0.6 0.8 A x vs. t 0 1 xt 0.2 0.4 0.6 0.8 5000 500550105015 4985 4990 4995 B x vs. t Figure6-4.Statevs.TimeontheFinalMeshforExample 6.3.3 for hp (2) – 2 and =10 4 126

PAGE 127

Strategy CPUTime(s) CollocPts MeshInterval MeshIterations NLPDensity 10 2 p 13.11 200 1 20 50.12 hp (1) 0.84 45 7 6 12.74 hp (2) – 4 0.70 76 19 7 5.05 hp (2) – 3 0.42 68 22 8 5.02 hp (2) – 2 0.81 46 22 8 6.27 h – 2 0.38 54 27 8 5.25 10 4 p – – – – – hp (1) 1.41 70 10 9 9.05 hp (2) – 4 1.14 161 39 6 2.50 hp (2) – 3 1.21 151 46 8 2.42 hp (2) – 2 2.74 126 45 10 2.82 h – 2 1.12 200 100 9 1.48 10 6 p – – – – – hp (1) 3.73 110 12 9 7.78 hp (2) – 4 8.42 258 59 7 1.64 hp (2) – 3 5.69 217 55 8 1.93 hp (2) – 2 4.13 188 53 8 2.24 h – 2 107.91 1368 684 10 0.22 Table6-10.SummaryofAccuracyandSpeedforExample 6.3.3 UsingVarious CollocationStrategiesandAccuracyTolerances. increased,the p –methodbecomesincreasinglylesscapableofaccuratelyap proximating thesolutiontothisproblem.For t f =50 ,the p –methodwasabletoapproximatethe solutionfor =(10 2 ,10 4 ) andfor t f =1000 the p –methodwasabletoapproximatethe solutionfor =10 2 .For t f =5000 ,the p –methodwasunabletoapproximatesolutionto theproblemforanyofthevaluesof The h – 2 –methodismuchmoreinsensitivethanthe p –methodtothevalueof t f for approximatingthesolutionofthisproblem.Usinga p –method,itisnotpossibletorene themeshspecicallyinregionswheretheerrormaybelarge. Becauserenement isnecessaryatthestartandterminusofthesolution,the h – 2 –methodismuchmore successfulthanthe p –methodinobtaininghigh-accuracyapproximationstothes olution ofthisproblem.However,becausethe h – 2 –methodonlyuseslow-orderapproximating polynomials,whenhighaccuracyisrequired,thenumberofa pproximatingpolynomials islarge.EventhoughtheNLPisverysparseforthe h – 2 –method,thesizeoftheNLP problemforhighlevelsofaccuracyresultsinacomputation allyinefcientNLPthat needstobesolved. 127

PAGE 128

Finally,both hp –methodswereabletoobtainhighaccuracysolutionsforthi s problemforallthevaluesof t f andallthevaluesof studied.First,whilebeingamore p –biasedmethod,the hp (1) –algorithmallowsforsubdivisionofapproximatinginterv als. Becauseoftheabilitytosubdividetheinitialsingle-inter valapproximation, hp (1) is abletodenselycollectintervalsonlyattheendsofthetraj ectoryregardlessofnal time.Furthermore,byallowinghigherdegreecollocationw hencomparedwiththe h – 2 –method,theproblemsizeismuchsmallerthanthe h – 2 –approximationforhigh levelsofaccuracy.Therefore, hp (1) resultsinanefcientapproximationforanyofthe accuracyrequirementspresented. hp (2) alsodemonstratesasimilarcapability.First, themeshisrenedonlyattheendsofthetrajectory.Second,b ecauseincreasesin thedegreeofapproximatingpolynomialareallowedwithini ntervals,highaccuracy approximationsaremadewithamuchsmallerNLPthanthe h – 2 –method.Forthis problem,itisdemonstratedthatbyallowingthenumberandl ocationofmeshintervals tovary and allowingthedegreeofapproximatingpolynomialtovary,co mputationally efcientapproximationsareobtainedforawiderangeofacc uracyrequirementswithan insensitivitytothedurationofthetrajectory. 6.4Example4:MotivatingExample2 NowwewillrevisittheoptimalcontrolproblemofSection 3.2 .Torestatethe problem,thisproblemistominimize J = Z t f t 0 udt (6–11) subjectto h = v v = g + u (6–12) theboundaryconditions h (0)=10, v (0)= 2 h ( t f )=0, v ( t f )=0, (6–13) 128

PAGE 129

andthecontrolpathconstraint 0 u 3 (6–14) where g =1.5 ,and t f isfree.Theoptimalsolutiontotheoptimalcontrolproblem given inEqs.( 6–11 )–( 6–14 )is ( h ( t ), v ( t ), u ( t ))= 8><>: ( 3 4 t 2 + v 0 t + h 0 3 2 t + v 0 ,0), t s ( 3 4 t 2 +( 3 s + v 0 ) t + 3 2 ( s ) 2 + h 0 3 2 t +( 3 s + v 0 ),3), t s (6–15) where s isgivenas s = t f 2 + v 0 3 (6–16) with t f = 2 3 v 0 + 4 3 r 1 2 v 2 0 + 3 2 h 0 (6–17) FortheboundaryconditionsgiveninEq.( 6–13 ),wehave ( s t f )=(1.4154,4.1641) Itisseenthattheoptimalcontrolforthisproblemis“bangbang”inthatitisatits minimum-valuefor t < s andatitsmaximumvaluefor t > s .Forthisproblem,the solutionisrepresentedbytwopiecewiselow-degreepolyno mials.Theonlydifculty, therefore,inaccuratelyapproximatingthesolutionofthi sproblemisthediscontinuity inthecontrolandthenon-smoothnessofthestateat s .Table 6-11 showstheinitial gridforeachapproximationstrategyandTable 6-12 summarizestheresultsfromeach approximationmethodforavarietyofvaluesof .Becausethesolutiontothisproblemis notsmooth,the p –methodisunabletoobtainanadequateapproximationforan yvalue of .Furthermore,the hp (1) methodperformssignicantlyworsethaneitherthe hp (2) – 2 – or h – 2 –methods.Moreover,for =10 6 ,nosolutionwasfoundbythe hp (1) method. Therefore,forthisproblem,the p –biasedmethodsunderperform. Next,foreachvalueof ,the hp (2) – 2 –and h – 2 –methodsperformedcomparably. Furthermore,foreachaccuracyrequirement,thesemethods providedhighlyefcient approximationstothesolution.Foranyaccuracylevel,the h –biasedstrategiesprovide 129

PAGE 130

Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 2 10 2 h – 2 10 2 Table6-11.InitialGridusedinExample 6.4 forEachStrategyUsingUniformly DistributedIntervalswith N k CollocationPoints. MaxAbs. Strategy CPU Colloc Mesh Mesh NLP Error Time(s) Pts Interval Iterations Density 10 2 1.76 10 2 p 4.68 200 1 20 33.78 3.02 10 3 hp (1) 0.88 110 22 22 2.97 2.74 10 2 hp (2) – 2 0.037 20 10 1 10.50 2.74 10 2 h – 2 0.037 20 10 1 10.50 10 4 p – – – – – 2.64 10 3 hp (1) 1.18 155 31 25 2.12 5.20 10 3 hp (2) – 2 0.15 44 21 8 5.18 4.01 10 3 h – 2 0.094 32 16 5 6.82 10 6 p – – – – – 3.57 10 4 hp (1) 26.57 580 80 40 0.92 3.88 10 5 hp (2) – 2 0.72 60 30 19 3.75 8.15 10 5 h – 2 0.30 54 27 8 4.15 Table6-12.SummaryofAccuracyandSpeedforExample 6.4 UsingVariousCollocation StrategiesandAccuracyTolerances. themostefcientapproximations.Fig. 6-5 showsthecontrolapproximationandthe exactsolutionfor =10 4 forthe h – 2 –approximation.Becausetheonlydifcult featureinthisproblemisinapproximatingthediscontinui tyinthecontrol,itisseenthat signicantrenementoccursnearthecontroldiscontinuit y. 6.5Example5:MinimumTime-to-ClimbofaSupersonicAircraft Considerthefollowingoptimalcontrolproblemthatisavar iationofthe minimum time-to-climbofasupersonicaircraft [ 72 74 75 ].Minimizethecostfunctional J = t f (6–18) subjecttothedynamicconstraints h = v sin r ,_ v = T D m g sin r ,_ r = g v ( n cos r ), (6–19) 130

PAGE 131

0 0 1.5 2 2 2.5 3 3 0.5 1 1 4 5 ut u ( t ) u ( t ) Figure6-5.Controlvs.TimeontheFinalMeshforExample 6.4 AlongwiththeExact Solution. theboundaryconditions h (0)=0, h ( t f )=19995 m v (0)=129.31 m/s v ( t f )=295.09 m/s r (0)=0 deg r ( t f )=0 deg (6–20) andtheinequalitypathconstraints r 45 deg (6–21) h 0 m where h isthealtitude, v isthespeed, r istheightpathangle, g isthelocalacceleration duetogravity,and n istheloadfactor(andisthecontrolforthisexample).Furt her 131

PAGE 132

detailsofthevehiclemodelandthenumericalvaluesofthec onstantsforthismodel canbefoundinRef.[ 72 ]and[ 75 ].Inthisproblem,duringsomesegmentsofthe trajectory,non-smoothfeaturesareencounteredbytheact ivepathconstraintsof Eq.( 6–21 ).Furthermore,inothersegmentsofthetrajectory,thetra jectoryissmooth. Figs. 6-6A – 6-6C showthestateapproximationtothesolutionofthisproblem for =10 4 usingthe hp (2) – 3 –method.First,itisseenbyexaminationofFig. 6-6A and 6-6C thatsignicantrenementoccursnearthelocationoftheac tivepathconstraints inaltitudeandightpathangle.Furthermore,itisnotedin thesmoothregionsofthe trajectorythatamuchmoresparsegridofpointsisused.Mos tofthecomputational effortisinapproximatingthenon-smoothfeatures(i.e.,t heswitchbetweenactivityand inactivityofthepathconstraints). Table 6-13 showstheinitialgridsforeachmethodandTable 6-14 showstheresults forapproximatingthesolutionofthisproblembyeachstrat egy.For =10 2 ,each methodperformscomparably.For =10 4 ,the p –methodwasunabletoobtaina solution.Furthermore,itisnotedthatfor200collocation points,the p –methodrequired approximately12minutesto not obtainanaccurateapproximation. Next,incomparingthe h – 2 hp (1) ,and hp (2) – x –methods,itisdemonstratedthat forthehighestaccuracyrequirementof =10 6 ,the hp (2) – x –methodsarethemost successfulinefcientlyapproximatingthesolutionofthi sproblem.The hp (2) – x –methods highlyrenethemesharoundthenon-smoothregionsofthetr ajectory.Furthermore, inthesmoothregionsofthetrajectory,higherdegreecollo cationisused.Therefore, the hp (2) – x –methodsareabletocombinethebenetsofsubdividinginte rvalsand increasingthedegreeofapproximatingpolynomial. hp (1) isslowerthan hp (2) becauseof thedensityoftheNLPandthenumberofmeshiterationsaregr eaterwhenusing hp (1) asopposedtousing hp (2) .Finally,forhighlevelsofaccuracy,the h – 2 –methodisslower than hp (2) becauseoftheexcessivesizeoftheNLP. 132

PAGE 133

0 0 5 10 100 50150 200 15 20 h (km)t (sec) 400 A h vs. t 0 100 100 50150 200 200 t (sec)v (m/s) 300 400 500 B v vs. t 0 0 10 10 40 100 50 50150 200 20 30 t (sec) r (deg) C r vs. t Figure6-6.Statevs.TimeontheFinalMeshforExample 6.5 for hp (2) – 3 and =10 4 Strategy N k NumberofIntervals p 10 1 hp (1) 10 1 hp (2) – 4 10 4 hp (2) – 3 10 3 hp (2) – 2 10 2 h – 2 10 2 Table6-13.InitialGridusedinExample 6.5 forEachStrategyUsingUniformly DistributedIntervalswith N k CollocationPoints. 133

PAGE 134

Strategy CPUTime(s) CollocPts MeshInterval MeshIterations NLPDensity 10 2 p 0.63 40 1 4 27.43 hp (1) 0.60 35 7 4 7.46 hp (2) – 4 0.37 48 12 2 5.01 hp (2) – 3 0.44 36 12 4 5.94 hp (2) – 2 0.25 24 12 3 7.73 h – 2 0.21 24 12 2 7.73 10 4 p 707.67 200 1 20 25.51 hp (1) 4.20 125 15 9 4.30 hp (2) – 4 1.74 86 21 6 2.89 hp (2) – 3 0.28 69 20 6 3.42 hp (2) – 2 0.95 78 29 6 2.85 h – 2 0.87 110 55 4 1.79 10 6 p – – – – – hp (1) 69.70 245 17 18 7.56 hp (2) – 4 14.75 217 50 8 1.20 hp (2) – 3 19.84 241 62 10 1.06 hp (2) – 2 19.09 359 86 7 0.73 h – 2 115.01 934 467 6 0.21 Table6-14.SummaryofAccuracyandSpeedforExample 6.5 UsingVariousCollocation StrategiesandAccuracyTolerances. 6.6DiscussionofResults Theresultsoftheveexamplesdemonstratesomeofthekeych aracteristicsof approximatingproblemsby hp –meshrenementalgorithm1and hp –meshrenement algorithm2.InExample 6.1 ,thesolutionissmoothandthe p and hp (1) methodsperform exactlythesame.Asexpected,themostefcientapproximati ontothisproblemisby a p –method.Furthermore,asexpected,fortightaccuracytole rances,the h – 2 –method iscomputationallyinefcient. hp –meshrenementalgorithm2hadcomparableCPU timestothe p –methodforallvaluesof .Byarelativelymodestincreaseindegreeof approximationcoupledwiththeincreasedsparsityoftheap proximationin hp –mesh renementalgorithm2,efcientapproximationstothesolu tionoftherstexample problemwereobtained. InExample2,thesolutiontotheproblemappearstobesmooth, butismore complicatedthanExample 6.1 .Forthisexample,astheaccuracytolerancewas tightened,the p –methodprovidedthebestapproximation.Furthermore,ast he accuracytolerancewastightened,the h – 2 –methodbecamecomputationallyinefcient. 134

PAGE 135

Surprisinglyforthisexample, hp –meshrenementalgorithm1waslessefcientthan hp –meshrenementalgorithm2.Theinefciencyof hp –meshrenementalgorithm1 islargelyduetotheexcessivenumberofmeshiterationsreq uiredtondasolution. Furthermore,whilethe p –methodprovidesthemostefcientsolutionapproximation tothisproblemforahighaccuracyrequirement,itisintere stingtonotethat hp –mesh renementalgorithm2doesnotsignicantlyunderperform. Whileasinglehigh-degree approximatingintervalprovidesthemostefcientapproxi mationtothesolutionofthis problem,itisinterestingtonotethatmanyintervalsofmod eratedegreealsoaccurately andefcientlyapproximatethesolutionofthisproblem.Th ecombinedaffectsofusing moderatedegreeintervalsandasparseNLPprovidescompara bleefciencytothe p –method. Example 6.3 studiesaproblemforthreedifferentvaluesof t f wherethechosen valuesof t f changetheinherenttimescaleoftheproblem.Inthisproble m,a p –method hasanincreasinglydifculttimendingasuitableapproxi mationasthenaltimeis increased.Furthermore,forveryhighaccuracyrequiremen ts,the h – 2 -methodisagain seentoprovideinefcientresults.Finally,both hp –algorithmsperformedwellforany valueof t f .Thecombinedabilitytorenethelocationofmeshinterval sandtheuseof avariabledegreepolynomialineachintervalmakethe hp –algorithmssuccessfulfor approximatingthesolutionofthisproblem. ThesolutionstructuretoExample 6.4 issignicantlydifferentfromtheprevious threeexamples.Thesolutiontothisproblemisrepresented bytwopiecewiselow-degree polynomials.Becausethecontrolisdiscontinuoustheconne ctionbetweenthesetwo piecewisepolynomialsisnon-smooth.Inthisproblem,a p –methodwasunabletond anaccurateapproximationforanyvalueof becauseofthenon-smoothnatureofthe solutiontothisproblem. hp –meshrenementalgorithm1performedsignicantlyworse whencomparedwith hp –meshrenementalgorithm2orthe h – 2 –methodandwasnot 135

PAGE 136

abletondasolutionforthetightestaccuracytolerance.F oranyvalueof hp –mesh renementalgorithm2orthe h – 2 –methodperformedwell. Thenalexampleisamodiedversionoftheminimumtime-toclimbofa supersonicaircraft.Thisexamplehasnon-smoothregionso fthetrajectoryandsmooth regionsofthetrajectory.Becauseofthenon-smoothnesstha texists,the p –method wasnotabletondapproximationsfortightaccuracytolera nces.Furthermore, hp –meshrenementalgorithm2outperformed hp –meshrenementalgorithm1and the h – 2 –method. hp –meshrenementalgorithm2moresuccessfullysimultaneou sly handlednon-smoothandsmoothfeaturesthan hp –meshrenementalgorithm1orthe h –method. Insummary,themethodthatwillbestapproximatesthesolut ionofanoptimal controlproblemis problemdependent .Ifaproblemissimpleandsmooth,a p –method performsthebest.Inanyothercase,a p –methoddoesnotperformwell. hp –mesh renementalgorithm1improvesupontherobustnessofastri ctly p –methodfor approximatingthesolutiontoproblemsthatmaybenon-smoo th.Forproblemswith nonsmoothfeatures, hp –meshrenementalgorithm1requiresmorecomputationalti me than hp –meshrenementalgorithm2orthe h – 2 –method. Byan h –method,low-accuracysolutionscanbeobtainedtoawideva rietyof problemsinacomputationallyefcientmanner.Ontheother hand,high-accuracy solutionsrequireatremendouscomputationalburdenbecau seofthelargeproblem sizeofthe h –methods. hp –meshrenementalgorithm2performedmuchbetter. Furthermore,for everyexampleproblemofthischapter hp –meshrenementalgorithm 2performedeitherbetterthantheothermethodsorwascompa rabletothebest methodforeachproblem(i.e., p forproblemswithsimplesmoothsolutionsor h for problemswithnonsmoothsolutionswithnohigh-orderbehav iors).Itisdemonstratedby hp –meshrenementalgorithm2thathighaccuracysolutionsto problemsdonotrequire high-degreepolynomials.Thehighestdegreepolynomialus edby hp –meshrenement 136

PAGE 137

algorithm2inanyoftheexamplesiseighthorder.Ifexponen tialconvergenceis expectedinaparticularinterval,ahigh-degreepolynomia lapproximationshould notberequired.Furthermore,becauselowandmid-degreepo lynomialsareused, hp –meshrenementalgorithm2resultsinanNLPthatissparse. Moreover,by specicallyreningonlypartsofanapproximationwiththe largesterrors,thedistribution ofthemeshfocusesonlyonthoselocationsthatneedmoremes hintervalsor mid-degreeapproximatingpolynomials. hp –meshrenementalgorithm2demonstrates signicantrobustnessirregardlessofaccuracyrequireme ntorproblemtypeonthe examplesinthisChapter.Asastartingpoint, hp –meshrenementalgorithm2utilizes low-degreepolynomialsthatensureasparseNLP.Renement ofthemeshoccursin regionsoflargeerror,andmid-degreepolynomialsareused onintervalswithsmooth solutions.Therefore,thecomputationalefciencyofan h –methodiscombinedwiththe high-accuracyofa p –method. 137

PAGE 138

CHAPTER7 AEROASSISTEDORBITALINCLINATIONCHANGEMANEUVER InthisChapter,anaeroassistedorbitalinclinationchang emaneuverforasmall maneuverablespacecraftinlow-Earthorbit(LEO)isstudied. Itisknownthatusing atmosphericforce(namely,liftanddrag)maysignicantly reducetheamountoffuel requiredtochangetheinclinationofthevehicleascompare dtoanall-propulsive maneuver[ 76 ].ThesurveypapersofRef.[ 77 ]and[ 78 ]summarizeearlyworkon aeroassistedorbitaltransfers.Ref.[ 79 ]studiedmultiple-passaeroassistedgliding maneuverswithalargeinclinationchangeofalargespacecr aftbetweengeostationary orbit(GEO)andLEO.AkeyresultofRef.[ 79 ]isthattheaeroassistedmaneuverwas muchmorefuelefcientthananall-propulsivemaneuver.Fu rthermore,inRef.[ 79 ], itwasdemonstratedforhighlyconstrainedheatingrateson thevehiclethatthefuel requiredtoaccomplishthemaneuverdecreasedasthenumber ofatmosphericpasses increased. Whilepreviousresearchonaeroassistedorbitalinclinatio nchangemaneuvershas beenforlargespacecraft,inthisresearch,weconsideraLEO toLEOaeroassisted orbitaltransferofasmallmaneuverablespacecraft.Theae roassistedmaneuveris comparedagainstall-propulsivemaneuvers.Themultiplepassaeroassistedglide maneuversareanalyzedasfunctionsofthenumberofatmosph ericpasses,therequired inclinationchange,andthemaximumallowableheatingrate onthevehicle.Theresults ofthisChapteraresimilartotheresultsfoundinRef.[ 80 ]. 138

PAGE 139

7.1EquationsofMotion Theequationsofmotionoverasphericalnon-rotatingEartha regiveninspherical coordinatesas[ 19 79 81 ] r = v sin r = v cos r cos r cos = v cos r sin r v = D m r 2 sin r r = 1 v L cos m r 2 v 2 r cos r = 1 v L sin m cos r v 2 r cos r cos tan (7–1) where r isthegeocentricradius, isthelongitude, isthelatitude, v isthespeed, r is theightpathangle, istheheadingangle, D isthedragforce, L istheliftforce, is theEarthgravitationalparameter,and isthebankangle.Next,itisnotedthatdragand liftaremodeled,respectively,as[ 19 ] D = 1 2 v 2 AC D (7–2) L = 1 2 v 2 AC L (7–3) where istheatmosphericdensity, A isthevehiclereferencearea, C D isthecoefcient ofdrag,and C L isthecoefcientoflift.Theatmosphericdensityismodele dasan exponentialoftheform = 0 exp( h ), (7–4) where 0 istheatmosphericdensityatsealevel, istheinverseofdensityscaleheight, and h = r R e isthealtitudewhere R e istheradiusoftheEarth.Theaerodynamic 139

PAGE 140

Table7-1.PhysicalandAerodynamicConstants. Quantity Value 0 1.225 kg/m 3 1 = 7200 m A 0.32 m C D 0 0.032 K 1.4 C L 0.5699 m 0 818 kg I sp 310 s ^ Q 19987 W/cm 2 h atm 110 km R e 6378145 m 3.986 10 14 m 3 /s 2 g 0 9.80665 m/s 2 coefcientsaremodeledusingastandarddragpolarofthefo rm[ 81 ] C L = C L C D = C D 0 + KC 2 L (7–5) where K isthedragpolarconstant, C D 0 isthezero-liftdragcoefcient,and C L isthelift slope.Inthisresearchwechooseagenerichighlift-to-dra gratiovehiclewhosemodel istakenfromRef.[ 82 ].Theaerodynamicdataforthevehicle,alongwiththephysi cal constants,aregiveninTable 7-1 7.1.1ConstraintsOntheMotionoftheVehicle Thefollowingconstraintsareenforcedonthemotionofthev ehicleduringight. Theinitialandnalconditionscorrespondtoacircularlow -Earthorbitofaltitude h 0 Theinitialinclinationiszerowhilethenalinclinationi s i f .Duringatmosphericight, constraintsareenforcedattheentranceandexitoftheatmo sphereontheightpath anglesuchthatuponatmosphericentrythevehiclemustdesc endintotheatmosphere whileuponatmosphericexitthevehiclemustascend.Furthe rmore,theangleofattack, stagnationpointheatingrate,andaltitudeareallconstra inedduringatmosphericight. 140

PAGE 141

7.1.1.1Initialandnalconditions Theinitialconditionsforthisstudycorrespondtoacircul arequatorialorbitof altitude h 0 .Theinitialorbitalconditionsforthisproblemaregiveni ntermsoforbital elements[ 83 ]as a ( t 0 )= R e + h 0 e ( t 0 )=0 i ( t 0 )=0,n( t 0 )=0 ( t 0 )=0, ( t 0 )=0, (7–6) where a isthesemi-majoraxis, i istheinclination, istheargumentofperigee, e isthe eccentricity, n isthelongitudeofascendingnode,and isthetrueanomaly.Thevalues for n( t 0 ), ( t 0 ), and ( t 0 ) arearbitrarilysettozero.Furthermore,thenalconditio nsfor thisproblemareforacircularorbitofaltitude h f = h 0 withaninclination i f andaregiven as a ( t f )= R e + h f e ( t f )=0, i ( t f )= i f (7–7) Becausethenalorbitiscircularandthelocationofthespac ecraftintheterminalorbit isnotconstrained, ( t f ) isundenedand ( t f ) and n( t f ) arefree. 7.1.1.2Interiorpointconstraints Alongwiththeinitialandnalconditionsonthevehicle,the followinginterior-point constraintsareenforcedateveryatmosphericentrance: r ( t atm 0 ) = h atm + R e r ( t atm 0 ) 0, (7–8) where t atm 0 isthebeginningofanyatmosphericightsegmentand h atm istheedge ofthesensibleatmosphere.Theseconstraintsensurethata ttheentranceofany atmosphericightsegmentthevehicleisdescendingintoth eatmosphere.Furthermore, attheexitofanyatmosphericightsegment,thefollowingc onstraintsareimposedsuch 141

PAGE 142

thatthevehicleisascendinguponatmosphericexit: r ( t atm f ) = h atm + R e r ( t atm f ) 0, (7–9) where t (atm) f isthetimeattheendofanyatmosphericightsegment. 7.1.1.3Vehiclepathconstraints Thefollowinginequalitypathconstraintsareenforceddur ingatmosphericight.The angleofattackisconstrainedas 0 max (7–10) Next,thestagnationpointheatingrateisconstrained.Usi ngtheChapmanequation[ 84 ], thestagnationpointheatingrate, Q ,isgivenas Q = ^ Q ( = 0 ) 0.5 ( v = v c ) 3.15 (7–11) where v c = p = R e .Therefore,duringatmosphericight,thefollowingconst rainton heatingrateisimposed: 0 Q Q max (7–12) Finally,duringatmosphericight,thealtitudeisconstra inedas: 0 h h atm (7–13) or,alternatively, R e r h atm + R e (7–14) 7.1.2TrajectoryEventSequence Thetrajectoryeventsequenceforamaneuverwith n atmosphericpassesisnow described.FromtheinitialconditionsofEq.( 7–6 ),adeorbitvelocityimpulse, V 1 ,is appliedtosendthevehicleintotheatmosphere.Thevehicle thenglidesinandoutof theatmosphere n times.Uponthenalatmosphericexit,aboostvelocityimpu lse, V 2 142

PAGE 143

isappliedtogivethevehicleenoughenergytoreachthenal orbitcondition.Upon reachingthenalorbitalaltitude,analrecircularizati onvelocityimpulse, V 3 isapplied torecircularizetheorbit.Thetrajectoryeventsequencei snowdescribedindetail(see Figs. 7-1A – 7-1C ): (i) Event: V 1 isappliedtotheinitialorbitalconditionsinEq.( 7–6 ); (ii) Phase: anexo-atmosphericglidethatterminatesattheedgeofthes ensible atmosphere( h atm ); (iii) Phase: anatmosphericglidethatterminatesattheedgeofthesensi bleatmosphere ( h atm ); Thefollowingtwophasesarethenrepeatedfortheremaining n 1 atmosphericpasses: (iv) Phase: anexo-atmosphericglidethatterminatesattheedgeofthes ensible atmosphere( h atm ); (v) Phase: anatmosphericglidethatterminatesattheedgeofthesensi bleatmosphere ( h atm ); Uponthenalatmosphericexit,thefollowingeventsandpha seoccur: (vii) Event: V 2 isappliedatthemaximumaltitudeofthesensibleatmospher e; (viii) Phase: anexo-atmosphericglidethatterminatesatthenalorbita laltitude; (ix) Event: V 3 isappliedatthenalorbitalaltitudetore-circularizeth eorbit. 7.1.3PerformanceIndex Thecostfunctionalforthemultiple-passaeroglidemaneuv eristominimizethefuel consumedduringthemaneuver.Thiscorrespondstominimizi ngthenegativeofthenal mass,i.e., J = m ( t f ). (7–15) Assumingimpulsivethrust, V k V k ,themassaftereachimpulseisthencomputed fromtheGoddardrocketequationas m + = m exp( V g 0 I sp ), (7–16) 143

PAGE 144

ADeorbitandFirstAtmosphericFlightSegment. BIntermediateExo-AtmosphericandAtmosphericFlightSegments. CFinalExo-AtmosphericFlightSegment. Figure7-1.SchematicofTrajectoryDesignforAeroassistedO rbitalTransferProblem. 144

PAGE 145

where m and m + arethevaluesofthemassimmediatelybeforeandafterappli cation oftheimpulse.Therefore,thenalmassafterthreevelocit yimpulsesisgivenas m ( t f )= m ( t 0 ) exp( ( V 1 + V 2 + V 3 )) g 0 I sp (7–17) 7.2Re-FormulationofProblem Themannerinwhichthemultiple-passaeroassistedorbital transferissolved differsslightlyfromtheformulationgiveninSection 7.1 .Insteadofusing and asthe controls,thequantities u 1 and u 2 areusedascontrols,where u 1 and u 2 aredenedas u 1 = C L sin u 2 = C L cos (7–18) Intermsof u 1 and u 2 C L and aregivenas C L = q u 2 1 + u 2 2 (7–19) = tan 1 ( u 1 = u 2 ). (7–20) Itisnotedthat u 1 and u 2 areusedinsteadofthemoreobviouschoicesofcoefcient ofliftandbankanglebecause u 1 and u 2 avoidthewrappingproblemthatarises whenusinganglesascontrols[ 19 ].Intermsof u 1 and u 2 ,thedifferentialequations ofEq.( 7–1 )aregivenas r = v sin r = v cos r cos r cos = v cos r sin r v = v 2 AC D 2 m r 2 sin r r = 1 v v 2 A 2 m u 2 r 2 v 2 r cos r = 1 v v 2 A 2 m cos r u 1 v 2 r cos r cos tan (7–21) 145

PAGE 146

Next,becausethecoefcientofliftislinearin ,theconstrainton inEq.( 7–10 )is equivalentlywrittenas 0 C L C L (7–22) where C L denotesthemaximumallowablevalueof C L .Therefore,thecontrols u 1 and u 2 areconstrainedas C L ( u 1 u 2 ) C L (7–23) 0 u 2 1 + u 2 2 C 2 L (7–24) Inthisstudythemaximumallowableangleofattackis max 40 degwhichcorresponds to C L =0.4 Next,theheatingrateconstraintofEq.( 7–11 )isrewrittentobebetterbehaved computationally.Specically,insteadofconstraining Q ,weconstrainthenatural logarithmof Q log Q ,as 1 log Q log Q max (7–25) Usingtheaboveformulations,theoptimalcontrolproblemc orrespondingtothe aeroassistedorbitaltransferproblemwassolvedincanoni calunitswherelengthis measuredinEarthradii,timeismeasuredinunitsof p R 3 e = ,speedismeasuredinunits of p = R e ,andmassismeasuredinunitsof m 0 ,theinitialmassofthevehicle.This systemofunitsisusedbecauseitiswell-scaled. 7.3OptimalControlProblem Theoptimalcontrolproblemcorrespondingtotheminimum-f uelaeroassistedorbital transferproblemisnowstated.Minimizethecostfunctiona lofEq.( 7–15 )subjecttothe dynamicconstraintsofEqs.( 7–21 ),thepathconstraintsofEqs.( 7–14 ),( 7–23 ),( 7–24 ), and( 7–25 ),andtheinteriorpointconstraintsofEqs.( 7–8 )and( 7–9 ).Theoptimal controlproblemissolvedusingthe hp –meshrenementalgorithmofRef.[ 67 ]whichis implementedinthesoftwareGPOPS[ 26 67 ]andusestheNLPsolverSNOPT[ 11 68 ]. 146

PAGE 147

AllNLPderivativeswerecomputedusingtheMatlabautomatic differentiator INTLAB [ 85 ]. 7.4All-PropulsiveInclinationChangeManeuvers Thefollowingtwoall-propulsivetransfersarestudiedino rdertocomparedwiththe aeroassistedorbitalinclinationchangemaneuver:(1)asi ngleimpulsedirectchangein inclinationviarotationoftheinitialvelocityand(2)atw oimpulsebi-parabolictransfer. Figs. 7-2A and 7-2B showthesingleimpulseandbi-parabolictransfers,respec tively.It isseenthatthesingleimpulsetransferrotatestheorbital velocityvectorbyanangle i Next,therstimpulseofthebi-parabolictransfersendsth eorbittoanorbitofinnite radiuswherethechangeininclinationisfree.Uponreturnt otheinitialorbit,asecond impulseisappliedtorecircularizetheorbit.Forinitiala ndterminalcircularorbitsofthe samesizebutwithdifferentinclinations,thedirectimpul seandtotalbi-parabolicimpulse requiredtochangeinclinationbyanamount i aregiven,respectively,as[ 86 ] V =2 v 0 sin i 2 V =2( p 2 1) v 0 (7–26) where v 0 istheinitialorbitalspeed.Byinspection,itisseenthatth ebi-parabolic transferismoreefcientfor i > 49 degrees.Therefore,whencomparingtheoptimal aeroassistedorbitaltransferfuelconsumptionwithall-p ropulsivemaneuvers,the aeroassistedorbitaltransferfuelconsumptioniscompare dwiththesingleimpulse maneuverwhen i 49 degandiscomparedwiththebi-parabolictransferwhen i > 49 deg. 7.5Results Themultiple-passaeroassistedorbitaltransferoptimalc ontrolproblemsummarized inSection 7.3 wassolvedfor n =(1,2,3,4) atmosphericpasses.Furthermore, foreachvalueof n ,achangeofinclinationof i f =(20,40,60,80) degsubjectto Q max =(400,800,1200, 1 ) W/cm 2 isstudied. 147

PAGE 148

ASingleImpulseTransfer. BBi-parabolicTransfer. Figure7-2.SchematicofAll-PropulsiveSingleImpulseandBi-Pa rabolicTransfers. 148

PAGE 149

7.5.1UnconstrainedHeatingRate Thecostofthemultiple-passaeroassistedorbitaltransfe risrstshownwithno constraintonheatingrate.Fig. 7-3A showstheminimumimpulse, V min ,required tochangetheinclinationoftheorbitalvehicleforone,two ,threeandfourpass aeroassistedmaneuversandtheall-propulsivemaneuvers. First,itisinterestingto notethatbecausethebi-parabolicmaneuverhasaconstantc ost,thecostforthe 60degand80degall-propulsiveinclinationchangesareequ ivalent.Next,itisseen forverysmallinclinationchange(i.e.,1deg),thereisvir tuallynodifferenceincost betweenaeroassistedmaneuversandall-propulsivemaneuv ers.For i f =(20,40,60,80) deg,itisseenthattheaeroassistedmaneuverscandramatic allyreducecostswhen comparedwiththeall-propulsivemaneuvers.Thecostforar equiredinclinationchange isinsensitivetothenumberofatmosphericpasses.For n =(1,2,3,4) ,thecostsare nearlyidenticalforanygivenvalue i f .ThisisseeninFig. 7-3B .Inthisgure,thenal massfractionaftercompletionoftheaeroassistedmaneuve risshown.Itisseenthat fortherangeofinclinationchanges,alargefractionofthe nalmassisretainedforthe vehicleafterthemaneuveriscompleted.Fora20degreeincl inationchange,morethan 60%ofthenalmassisretained.Furthermore,foran80degre einclinationchange maneuver,approximately30%ofthenalmassisretainedaft ercompletionofthe maneuver. Next,Figs. 7-4A – 7-4D showtheimpulsesfor V 1 V 2 and V 3 alongwiththetotal lossofvelocityfromdrag, V D .First,itisnotedthatthedeorbitandrecircularization impulses, V 1 and V 3 ,areverysmall.Regardlessoftheamountofinclinationcha nge thatisdesired,theseimpulsesaremuchsmallerthantheboo stimpulse, V 2 .When comparingFigs. 7-4B and 7-4D ,itisseenthatthemajorityoftheboostimpulseisused torecovervelocitylostduetoatmosphericdrag.Therefore ,asexpected,thecostofthe aeroassistedmaneuverislargelyduetorecoveringenergyl ostduetodrag. 149

PAGE 150

Total V (m/s)InclinationChange(deg) One-Pass Two-Pass Three-Pass Four-Pass All-Propulsive 0 0 2000 4000 40 60 80 20 6000 8000 max AMinimumimpulse, V min vs.nalinclination, i f 1.2 0 1 1 3 2 4 0.8 0.6 0.4 0.2 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPassesm ( t f ) = m 0 1.4 BFinalmassfraction, m ( t f ) = m 0 ,vs.numberofatmosphericpasses, n Figure7-3. V min ,vs. i f ,and m ( t f ) = m 0 vs. n 150

PAGE 151

1 3 28 32 34 36 38 30 24 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses V 1 (m/s) A V 1 ,vs. n 1 3 1000 2000 3000 24 4000 5000 6000 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses V 2 (m/s) B V 2 ,vs. n 0 1 3 30 50 24 10 40 60 20 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses V 3 (m/s) C V 3 ,vs. n 1 3 1000 2000 3000 24 4000 5000 6000 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses V D (m/s) D V D ,vs. n Figure7-4. V 1 V 2 V 3 ,and V D Impulsesvs. n Figs. 7-5A showsthetotalinclinationchangeachievedbyatmospheric forcesasa functionofnumberofpassesandrequiredinclinationchang e.Itisseenthatforeach requiredinclinationchange,themajorityoftheinclinati onchangeisperformedbythe atmosphere.Itisinterestingtonotefor i f =80 degthatthefour-passmaneuverhas lessinclinationchangefromtheatmosphereascomparedwit htheone,twoorthree passmaneuvers.Nevertheless,thecostofthe i f =80 degmaneuverisnearlythesame asthecostoftheone,twoorthreepassmaneuvers.Next,Fig. 7-5B showsthetotal transfertimefor n =(1,2,3,4) for i f =(20,40,60,80) .Asexpected,thetotaltransfer 151

PAGE 152

0 1 3 24 40 60 80 100 120 20 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses i aero (deg) A i aero vs. n 0 1 3 3.5 1.52.5 2 4 200 400 100 300 i ( t f )=20deg i ( t f )=40deg i ( t f )=60deg i ( t f )=80deg NumberofAtmosphericPasses TrajectoryDuration(minutes) B T vs. n Figure7-5.Changein i byAtmosphere, i aero ,andtransfertime, T ,vs. n for Q = 1 timeincreasesasthenumberofatmosphericpassesincrease s.Finally,itisnotedfor aone-passmaneuver,theaeroassistedtransferrequiresap proximately100minutes regardlessofdesiredinclinationchange. 152

PAGE 153

7.5.2ConstrainedHeatingRate Theeffectofconstrainingthestagnationpointheatingrat eonthesmallspacecraft isnowexamined.Figs. 7-6A – 7-6D showsthetotalrequiredimpulse, V min asa functionof i f and Q max for i f =(20,40,60,80) Q max =( 1 ,1200,800,400) W/cm 2 and n =(1,2,3,4) .Itisseenthat,foranyparticularvalueof i f V min increasesas Q max decreases.Furthermore,itisnotedagainthatthetotalcos tisinsensitivetothe numberofatmosphericpasses.Next,forevenahighlyconstr ainedheatingrate,the aeroassistedmaneuverismorefuelefcientthananall-pro pulsivemaneuver.Theonly caseforwhichtheaeroassistedmaneuverhasagreatercostt hantheall-propulsive maneuverisfor i f =80 degreesand Q =400 W/cm 2 .Figs. 7-7A – 7-7D showthenal massfractionaftercompletionofthemaneuversfor i f =(20,40,60,80) deg.Again,it isseenthatas Q max decreases,morefuelisrequiredtoachieveadesiredinclin ation change. Finally,Figs. 7-8A – 7-8D showthemaximumvalueofthesealevelgravity-normalized specicforce, S max ,experiencedbythevehicleduringatmosphericight,wher e S is denedas S = v 2 A 2 mg 0 q C 2 D + C 2 L (7–27) Foranunconstrainedheatingrate, S max rangesfromapproximatelythreetotwenty-ve. Astheheatingrateconstraintistightened,themagnitudeof S max decreasessignicantly forallvaluesofinclinationchange.For Q =400 W/cm 2 S max 1 foranydesired inclinationchange.7.5.3OptimalTrajectoryStructure Thetrajectoryoftheaeroassistedorbitalinclinationcha ngemaneuverisnow examinedindetail.A40degreeinclinationchangemaneuver isusedtoexaminethe structureoftheoptimaltrajectory.Figs. 7-9A – 7-9C showthealtitudeversusvelocity, heatingrateversustime,andaltitudeversusightpathang lefor Q max =( 1 ,800,400) W/cm 2 for i f =40 degforaonepassmaneuver.Itisseenforanunconstrainedhe ating 153

PAGE 154

All-Propulsive 1 3 1000 2000 3000 24 4000 5000 NumberofAtmosphericPasses V min (m/s) Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) A V min vs. n for i f =20 deg. All-Propulsive 1 3 2000 24 4000 6000 8000 NumberofAtmosphericPasses V min (m/s) 10000 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) B V min vs. n for i f =40 deg. All-Propulsive 1 3 12000 24 4000 6000 8000 NumberofAtmosphericPasses V min (m/s) 10000 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) C V min vs. n for i f =60 deg. All-Propulsive 1 3 12000 14000 24 4000 6000 8000 NumberofAtmosphericPasses V min (m/s) 10000 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) D V min vs. n for i f =80 deg. Figure7-6. V min ,vs. n forVariousMaxHeatingRatesandFinalInclinations. rate,thevehicledivesdeeplyintotheatmosphereonceandt henascendsbackoutof theatmosphere.Forthecasesof Q max =(800,400) W/cm 2 ,itisseenthatthevehicle descendsuntilreachingthemaximumvalueoftheheatingrat econstraint.Thevehicle thenglidesatthemaximumvalueoftheheatingrateconstrai ntforashortduration beforere-ascendingtoneartheedgeoftheatmosphere.Thev ehiclethendescends againtoreachtheheatingrateconstraintbeforereascendi ngoutoftheatmosphere. Therefore,foraone-passmaneuverwithanactiveconstrain tontheheatingrate,the trajectoryessentiallybecomesatwo-passmaneuver.Figs. 7-10A and 7-10B showthe 154

PAGE 155

All-Propulsive 0 1 1 3 0.5 1.5 24 NumberofAtmosphericPassesm ( t f ) = m 0 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) A m ( t f ) = m 0 vs.nfor i f =20 deg. 1.2 All-Propulsive 0 1 1 3 24 0.8 0.6 0.4 0.2 NumberofAtmosphericPassesm ( t f ) = m 0 1.4 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) B m ( t f ) = m 0 vs.nfor i f =40 deg. All-Propulsive 0 1 1 3 24 0.8 0.6 0.4 0.2 NumberofAtmosphericPassesm ( t f ) = m 0 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) C m ( t f ) = m 0 vs.nfor i f =60 deg. All-Propulsive 0 1 1 3 24 0.8 0.6 0.4 0.2 NumberofAtmosphericPassesm ( t f ) = m 0 Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) D m ( t f ) = m 0 vs.nfor i f =80 deg. Figure7-7. m ( t f ) = m 0 ,vs. n forVarious Q max ,and i f =(20,40,60,80) deg. angleofattackandbankangleforthevehiclefortheone-pas smaneuver.Itisseenthat foreachvalueofheatingrateconstraint,thevehicleenter sandexitstheatmosphere withanearmaximumvalueofangleofattackofapproximately 40deg.Furthermore,it isseenthatthevehicleentersatanearly90degbankangle,r otates180degrees,and exitsatanearly-90degbankangle. Figs. 7-11A – 7-11D showthealtitudeversusvelocityforone,two,threeandfou r passmaneuversfor i f =40 degreesand Q max =400 W/cm 2 .Itisinterestingtonotethat theprolesarenearlyidenticalregardlessofnumberofpas ses.Fortheone-pass 155

PAGE 156

0 1 3 8 2 2 4 4 6 NumberofAtmosphericPasses Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 )S max (Dimensionless) A S max vs. n for i f =20 deg. 0 1 3 5 15 24 10 NumberofAtmosphericPasses Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 )S max (Dimensionless) B S max vs. n for i f =40 deg. 0 1 3 5 25 35 30 15 24 10 20 NumberofAtmosphericPasses Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 )S max (Dimensionless) C S max vs. n for i f =60 deg. 0 1 3 30 50 24 10 40 20 NumberofAtmosphericPasses Q max = 1 (W/cm 2 ) Q max =1200 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 )S max (Dimensionless) D S max vs. n for i f =80 deg. Figure7-8. S max ,vs. n ,forVarious Q max and i f =(20,40,60,80) deg. maneuver,thevehicledescendsandascendstonearlytheedg eofthesensible atmosphere.Atthispoint,thevehicledescendsagainbefore exitingtheatmosphere. Forthetwo-passmaneuver,thevehiclesimplysplitstheone -passmaneuverbetween eachpass.Forthethreeandfour-passmaneuvers,thesemane uverscollapseintoa two-passmaneuver.Thesecondatmosphericpassforthethre e-passmaneuverand thesecondandthirdatmosphericpassesforthefourpassman euvernevergodeeply intotheatmosphere.Regardlessofnumberofpasses,thesam egeneraltwo-pass proleisexhibited.InexaminingFigs. 7-12A – 7-12D ,thesimilarityofthemaneuvers 156

PAGE 157

12 Altitude(km) Velocity(km/s) 5 7 11 8 9 6 10 40 60 80 100 120 20 Q max = 1 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) A h ,vs. v Time(sec)_ Q (W/cm 2 ) 0 0 1000 1000 1500 1500 2000 2000 2500 2500 3000 3000 3500 500 500 Q max = 1 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) B Q ,vs. t for i f =40 -deg. Altitude(km) r (deg) 0 5 5 10 40 60 80 100 120 20 Q max = 1 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) C h vs. r Figure7-9. h vs. v Q vs. t ,and h vs. r for i f =40 -deg. regardlessofnumberofpassesisseenagaininexaminingthe heatingrateversustime intheatmosphere.Regardlessofnumberofpasses,thesameh eatingrateproleis seen.Thersttimethemaximumvalueoftheheatingrateisme t,thevehiclerides alongthemaximumvalueoftheheatingrateandthenreascend s.Thesecondtimethe vehicleachievesthemaximumvalueoftheheatingrateconst raint,itonlytouchesit beforereascending.Again,itisseenthattheinteriorpasse sforthethreeandfour-pass maneuverscollapseintopassesthatdonotdescendfarintot heatmosphere. 157

PAGE 158

Time(sec) AngleofAttack(Deg) 0 0 30 50 1000 1500 2000 2500 3000 3500 10 500 40 20 Q max = 1 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) A vs. t Time(sec) BankAngle(Deg) 0 0 100 200 1000 1500 2000 2500 3000 3500 200 500 100 Q max = 1 (W/cm 2 ) Q max =800 (W/cm 2 ) Q max =400 (W/cm 2 ) B vs. t Figure7-10. and vs. t for n =1 and i f =40 deg. 158

PAGE 159

Altitude(km) 7 7.5 8 8.5 9 9.5 40 60 80 100 120 Velocity(km/s) A h vs. v for n =1 and i f =40 deg. Altitude(km) 7 7.5 8 8.5 9 9.5 40 60 80 100 120 Velocity(km/s) FirstPass SecondPass B h vs. v for n =2 and i f =40 deg. Altitude(km) 7 7.5 8 8.5 9 9.5 40 60 80 100 120 Velocity(km/s) FirstPass SecondPass ThirdPass C h vs. v for n =3 and i f =40 deg. Altitude(km) 7 7.5 8 8.5 9 9.5 40 60 80 100 120 Velocity(km/s) FirstPass SecondPass ThirdPass FourthPass D h vs. v for n =2 and i f =40 deg. Figure7-11. h vs. v for n =(1,2,3,4) and i f =40 degand Q max =400 (W/cm 2 ). 7.6Discussion Theresultsofthisstudydemonstrateseveralinterestingf eaturesofusing aeroassistedglidemaneuverstochangetheinclinationofa smallspacecraftinLEO.For anunconstrainedheatingrate,itisdemonstratedthatthet otalfuelrequiredtochange theinclinationofthesmallspacecraftismuchlessthanfor anall-propulsivemaneuver. Astheheatingrateistightened,thecostoftheaeroassisted maneuversincrease,but, arestillmoreefcientthantheall-propulsivemaneuverse xceptforthesinglecaseof i f =80 degreesand Q max =400 W/cm 2 .Foranunconstrainedheatingrate,regardless 159

PAGE 160

_ Q (W/cm 2 ) 0 0 1000 1500 2000 2500 3000 200 400 500 500 100 300 TimeSinceAtmosphericEntry(sec) A Q vs. t for n =1 and i f =40 deg. Q (W/cm 2 ) 0 0 1000 1500 2000 2500 3000 200 400 500 500 100 300 FirstPass SecondPass TimeSinceAtmosphericEntry(sec) B Q vs. t for n =2 and i f =40 deg. Q (W/cm 2 ) 0 0 1000 1500 2000 2500 3000 200 400 500 500 100 300 FirstPass SecondPass ThirdPass TimeSinceAtmosphericEntry(sec) C Q vs. t for n =3 and i f =40 deg. Q (W/cm 2 ) 0 0 1000 2000 3000 200 400 4000 5000 500 100 300 FirstPass SecondPass ThirdPass FourthPass TimeSinceAtmosphericEntry(sec) D Q vs. t for n =4 and i f =40 deg. Figure7-12. Q ,vs. t for n =(1,2,3,4) i f =40 degand Q max =400 (W/cm 2 ). ofnumberofatmosphericpasses,thevehicledivesdeeplyin totheatmosphereonly onceandaffectsalloftheinclinationchangewiththissing leatmosphericdive.For theconstrainedheatingratecases,itisdemonstratedthat regardlessofnumberof atmosphericpasses,thevehicleessentiallyconductsatwo -passmaneuver.Inthe rstpass,thevehicledivesintotheatmosphereandridesth eheatingrateconstraint foraniteduration.Inthesecondpass,thevehiclemerelyt ouchestheheatingrate constraintbeforereascending. 160

PAGE 161

7.7Conclusions Astudyhasbeengivenonchangingtheinclinationofasmallh ighlift-to-dragratio spacecraftinlow-Earthorbitusingaeroassistedmaneuvers .Threevelocityimpulses (ade-orbit,boost,andre-circularization)areusedalong withliftanddragforcesinthe atmospheretochangetheinclinationofthevehicle.Itisfo undthattheoptimalsolution isinsensitivetothenumberofatmosphericpasses.Foranun constrainedheatingrate, thesolutionreducestoaone-passmaneuverregardlessofnu mberofpasses.For highlyconstrainedheatingrates,thesolutionreducestoa two-passmaneuver.Finally, itwasfoundthatforallbutthehighestinclinationchanges withthetightestheatingrate constraintsthatanaeroassistedmaneuverwasmorefuelef cientthananall-propulsive maneuver. 161

PAGE 162

CHAPTER8 CONCLUSIONS 8.1DissertationSummary Numericalmethodsforsolvingoptimalcontrolproblemsfal lintotwocategories: indirectanddirectmethods.Anindirectmethodsolvestheco ntinuous-timerst-order necessaryconditionsfromthecalculusofvariations.While indirectmethodsingeneral providehigh-accuracyapproximations,itisoftenhardto ndasolutiontothese necessaryconditions.Directcollocationmethods,ontheo therhand,areamore robustnumericalapproximationtechnique.Directcolloca tionmethodsdiscretizethe originalcontinuous-timeoptimalcontrolproblemsuchtha ttheinnite-dimensional continuous-timeproblemistranscribedintoanite-dimen sionaldiscreteproblem,or nonlinearprogrammingproblem(NLP).Directmethodshavemo stcommonlybeen appliedas h –methods.An h –methodusesmanylow-degreeapproximatingsubintervals andguaranteesthattheNLPissparse.Incontrastto h –methods, p –methodsuseafew xednumberofintervalswithavariabledegreepolynomiali neachinterval.A p –method resultsinadenseNLP,but,ifthesolutiontoaproblemissim pleandsmooth,high accuracyapproximationstothecontinuous-timeoptimalco ntrolproblemareobtained withsmallNLPs. Inthisdissertation,an hp –directcollocationmethodhasbeenpresented.The hp –methodallowsforavariablenumberofapproximatinginter valswithavariable degreeapproximationwithineachinterval.Theresultisam ethodthatisrobustwith respecttothesolutionstructureofaproblem(thatis,some optimalcontrolproblemsare betterapproximatedbyan h –biasedmethod,whileotheroptimalcontrolproblemsare betterapproximatedbya p –biasedmethod).Thediscreteapproximationineachinterv al isdenedbythesetofLegendre-GaussorLegendre-Gauss-Ra daupoints.These pointsareusedbecausetheyaredenedbyhighlyaccuratequ adraturerulesandwhen usedasthesupportpointsinaLagrangepolynomial,donotex hibitRungephenomenon 162

PAGE 163

forhigh-degreeapproximations.Next,inanalyzingtheKKTco nditionsoftheNLP,it isfoundthatthe hp –methodofthisresearchresultsinadiscreterepresentati onofthe continuous-timerst-orderoptimalityconditionsifthec ostateiscontinuous. Two hp –meshrenementalgorithmshavebeenpresented.Ineachalg orithm,the accuracyoftheapproximationineachmeshintervalisasses sedbyevaluatingthe differential-algebraicconstraintsbetweencollocation points.Ifanapproximatinginterval isinaccurate,theintervalismodiedbyeitherincreasing thedegreeofapproximation withintheinterval,orbysubdividingtheinterval.The hp –meshrenementalgorithms iterativelyupdatetheapproximatinggriduntilanapproxi mationisaccurateinevery approximatinginterval.The hp –meshrenementalgorithmsarecapableofutilizing thebestapproximationstrategyfordifferentregionsofas olutionbymodifyingthe mesheitherbyincreasingthedegreeoftheapproximatingpo lynomialorincreasing thenumberofmeshintervals.Thetwo hp –meshrenementalgorithmspresented inthisresearcharedistinctlydifferentfromoneanotherb ecausetherstalgorithm isa p –biasedalgorithmwhilethesecondalgorithmisan h –biasedalgorithm.Both algorithmsarecomparedwithoneanotheralongwithcompari sonsto p and h –methods overseveralexamples.Itisdemonstratedthatthe hp –methodsarerobusttoproblem typeandaccuratelysolvetheproblemspresentedinacomput ationallyefcientmanner. Furthermore, hp –meshrenementalgorithm2isthemostsuccessfulstrategy when comparedwith p h ,and hp –meshrenementalgorithm1overtherangeofproblems. Utilizingthecomputationalsparsityofan h –methodasastartingpoint, hp –mesh renementalgorithm2resultsinasparseNLP.Furthermore, byusingmid-degree collocationonlyonsmoothintervals,highaccuracyapprox imationsareobtainedina computationallyefcientmanner. hp –meshrenementalgorithm2providesaccurate approximationstocontinuous-timeoptimalcontrolproble mswithsimultaneouslysmall andsparseNLPs. 163

PAGE 164

Finally,anoptimalaeroassistedorbitalinclinationchan gemaneuverisstudied. Inthisstudy,asmallmaneuverablespacecraftinlow-Eartho rbit(LEO)usesthe aerodynamicforcesofliftanddragtochangetheorbitalinc lination.Themaneuver isstudiedforamultiplenumberofatmosphericpassesandva riableheatingrate constraints.Theaeroassistedmaneuversisdemonstratedt obemorefuelefcient thanall-propulsiveinclinationchangemaneuvers.Furthe rmore,itisfoundthatthe aeroassistedmaneuversareinsensitivetonumberofatmosp hericpasses. 8.2FutureWork 8.2.1ContinuedWorkonMeshRenement Two hp –meshrenementalgorithmshavebeenpresentedinthisrese arch.Itis believedthatautomaticmeshgenerationisarichsubjectan dmanymoreadvances maybemade.Typically,meshrenementstrategiesarestudi edfor h –methods.In an h –method,thereareonlytwodegreesoffreedominreningame sh,namely,the numberofmeshintervalsandtheplacementofthemeshpoints hp –meshrenement algorithmsaresignicantlymorecomplex,however,becaus etherearethreedegrees offreedomforreningamesh,namely,thenumberofmeshinte rvals,theplacementof themeshpoints, and thedegreeoftheapproximatingpolynomialwithineachinte rval. Becauseofthisextradegreeoffreedomingeneratinganappro ximatingmesh,the choiceofhowtoreneameshismorecomplicated,butalsomor einterestingand fruitful.Presentedinthisresearcharetwostaticrenemen talgorithmsthatonly increasethesizeoftheNLPoneveryiteration.Twointerest ingfutureconceptsto explorearenowstated.Therstconceptisarenementproce ssto decrease thesize oftheNLPifanyintervalsaremoreaccuratethanthedesired accuracylevel,further increasingthecomputationalefciencyoftheNLP.Theseco ndconceptistoadda dynamicmeshrenementcomponentthatcanbeusedtoexactly locatetheswitchin activitybetweenactiveandinactivepathconstraints.One ofthestrengthsofadirect collocationmethodisthatthelocationofactiveandinacti veconstraintarcsdoesnot 164

PAGE 165

needtobeknowninordertoapproximatethesolutionofthepr oblem.Ifthelocationsof thesearcscanbedetermined,however,theapproximationsc anbemadesignicantly moreaccurate.8.2.2DiscontinuousCostate Ithasbeendemonstratedinthisresearchthathigh-accurac yapproximationsare achievedbytheproposed hp –methodifthecostateiscontinuous.Ifmeshpointsare atthelocationofdiscontinuityinthecostate,thetransfo rmedadjointsystemisan inexactdiscreterepresentationofthecontinuous-timer st-ordernecessaryconditions. Foradynamicrenementalgorithmthatwillexactlylocatet heswitchinactivityof inequalitypathconstraints,itislikelythatthecostatem aybediscontinuousatmesh points.Itisnecessarytodeterminehowtomakethetransfor medadjointsystema discreterepresentationofthecontinuous-timerst-orde rnecessaryconditionsfora discontinuouscostatesolutionifmeshpointsareattheloc ationofdiscontinuityinthe costate. 165

PAGE 166

REFERENCES [1]Kirk,D.E., OptimalControlTheory:AnIntroduction ,DoverPublications,Mineola, NewYork,2004. [2]Bryson,A.E.andHo,Y.-C., AppliedOptimalControl ,HemispherePublishing,New York,1975. [3]Pontryagin,L.S., MathematicalTheoryofOptimalProcesses ,JohnWiley&Sons, NewYork,1962. [4]Betts,J.T.,“SurveyofNumericalMethodsforTrajectoryO ptimization,” Journalof Guidnance,Control,andDynamics ,Vol.21,No.2,March–April1998,pp.193–207. [5]vonStryk,O.andBulirsch,R.,“DirectandIndirectMethod sforTrajectory Optimization,” AnnalsofOperationsResearch ,Vol.37,1992,pp.357–373. [6]Oberle,H.J.andGrimm,W.,“BNDSCO:AProgramfortheNumeri calSolution ofOptimalControlProblems,”Tech.rep.,InstituteofFligh tSystemsDynamics, GermanAerospaceResearchEstablishmentDLR,IB515-89/22,O berpfaffenhofen, Germany,1990. [7]Press,W.H.,Flannery,B.P.,Teukolsky,S.,andVetterli ng,W.T., Numerical Recipes:TheArtofScienticComputing ,CambridgeUniversityPress,1990. [8]Keller,H.B., NumericalSolutionofTwoPointBoundaryValueProblems ,SIAM, 1976. [9]Dickmanns,E.D.andWell,K.H.,“ApproximateSolutionofOpt imalControl ProblemsUsingThirdOrderHermitePolynomialFunctions,” LectureNotesin ComputationalScience,Springer ,Vol.27,1975. [10]Russell,R.D.andShampine,L.F.,“ACollocationMethod forBoundaryValue Problems,” NumericalMathematics ,Vol.19,No.1,1972,pp.1–28. [11]Gill,P.E.,Murray,W.,andSaunders,M.A.,“SNOPT:AnSQPAlgori thmfor Large-ScaleConstrainedOptimization,” SIAMReview ,Vol.47,No.1,January 2002,pp.99–131. [12]Biegler,L.T.andZavala,V.M.,“Large-ScaleNonlinearPr ogrammingUsingIPOPT: AnIntegratingFrameworkforEnterprise-WideOptimization,” Computersand ChemicalEngineering ,Vol.33,No.3,March2008,pp.575–582. [13]Byrd,R.H.,Nocedal,J.,andWaltz,R.A.,“KNITRO:AnIntegr atedPackagefor NonlinearOptimization,” LargeScaleNonlinearOptimization ,SpringerVerlag, 2006,pp.35–59. 166

PAGE 167

[14]Kraft,D., OnConvertingOptimalControlProblemsintoNonlinearProgra mming Codes ,Vol.F15of NATOASISeries, in ComputationalMathematicalProgramming ,Springer-Verlag,Berlin,1985,pp.261–280. [15]Bock,H.G.andPlitt,K.J.,“AMultipleShootingAlgorithmfo rDirectSolutionof OptimalControlProblems,” IFAC9thWorldCongress ,Budapest,Hungary,1984. [16]Schwartz,A., TheoryandImplementationofNumericalMethodsBasedonRung eKuttaIntegrationforSolvingOptimalControlProblems ,Ph.D.thesis,Departmentof ElectricalEngineering,UniversityofCalifornia,Berkeley, 1996. [17]Brauer,G.L.,Cornick,D.E.,andStevenson,R.,“Capabili tiesandApplicationsof theProgramtoOptimizeandSimulateTrajectories,”Tech.Rep .NASA-CR-2770, NationalAeronauticsandSpaceAdministration,1977. [18]Hargraves,C.R.andParis,S.W.,“DirectTrajectoryOp timizationUsingNonlinear ProgrammingTechniques,” JournalofGuidance,Control,andDynamics ,Vol.10, No.4,July–August1987,pp.338–342. [19]Betts,J., PracticalMethodsforOptimalControlandEstimationUsingNo nlinear Programming ,SIAMPress,Philadelphia,2nded.,2009. [20]Enright,P.J.andConway,B.A.,“DiscreteApproximations toOptimalTrajectories UsingDirectTranscriptionandNonlinearProgramming,” JournalofGuidance, Control,andDynamics ,Vol.19,No.4,July–August1996,pp.994–1002. [21]Hager,W.W.,“Runge-KuttaMethodsinOptimalControla ndtheTransformed AdjointSystem,” NumerischeMathematik ,Vol.87,2000,pp.247–282. [22]Benson,D.A., AGaussPseudospectralTranscriptionforOptimalControl ,Ph.D. thesis,DepartmentofAeronauticsandAstronautics,Massach usettsInstituteof Technology,Cambridge,Massachusetts,2004. [23]Huntington,G.T.,Benson,D.A.,andRao,A.V.,“OptimalCo ngurationof TetrahedralSpacecraftFormations,” TheJournaloftheAstronauticalSciences Vol.55,No.2,April-June2007,pp.141–169. [24]Betts,J.T.andHuffman,W.P.,“SparseOptimalControlSof tware–SOCS,”Tech. Rep.MEA-LR-085,BoeingInformationandSupportServices,Seattl e,Washington, July1997. [25]vonStryk,O., User'sGuideforDIRCOLVersion2.1:ADirectCollocationMe thod fortheNumericalSolutionofOptimalControlProblems ,TechnischeUniversitat Darmstadt,Darmstadt,Germany,1999. 167

PAGE 168

[26]Rao,A.V.,Benson,D.A.,Darby,C.L.,Francolin,C.,Patte rson,M.A.,Sanders, I.,andHuntington,G.T.,“Algorithm902:GPOPS,AMatlabSoftwa reforSolving Multiple-PhaseOptimalControlProblemsUsingtheGaussPseud ospectral Method,” ACMTransactionsonMathematicalSoftware ,Vol.37,No.2,April–June 2010,pp.22:1–22:39. [27]Dontchev,A.L.,Hager,W.W.,andMalanowski,K.,“ErrorBou ndsfortheEuler ApproximationandControlConstrainedOptimalControlProbl em,” Numerical FunctionalAnalysisandApplications ,Vol.21,2000,pp.653–682. [28]Dontchev,A.L.,Hager,W.W.,andVeliov,V.M.,“Second-O rderRunge-Kutta ApproximationsInConstrainedOptimalControl,” SIAMJournalonNumerical Analysis ,Vol.38,2000,pp.202–226. [29]Dontchev,A.L.andHager,W.W.,“TheEulerApproximationi nStateConstrained OptimalControl,” MathematicsofComputation ,Vol.70,2001,pp.173–203. [30]Zhao,Y.andTsiotras,P.,“ADensity-FunctionBasedMes hRenementAlgorithm forSolvingOptimalControlProblems,”InfotechandAerospace Conference,AIAA Paper2009-2019,Seattle,Washington,April2009. [31]Jain,D.andTsiotras,P.,“TrajectoryOptimizationUs ingMultiresolutionTechniques,” JournalofGuidance,Control,andDynamics ,Vol.31,No.5,September-October 2008,pp.1424–1436. [32]Betts,J.T.andErb,S.O.,“OptimalLowThrustTrajectori estotheMoon,” SIAM J.Appl.Dyn.Syst ,Vol.2,2003,pp.144–170. [33]Elnagar,G.,Kazemi,M.,andRazzaghi,M.,“ThePseudospec tralLegendreMethod forDiscretizingOptimalControlProblems,” IEEETransactionsonAutomatic Control ,Vol.40,No.10,1995,pp.1793–1796. [34]Elnagar,G.andRazzaghi,M.,“ACollocation-TypeMetho dforLinearQuadratic OptimalControlProblems,” OptimalControlApplicationsandMethods ,Vol.18, No.3,1998,pp.227–235. [35]Benson,D.A.,Huntington,G.T.,Thorvaldsen,T.P.,andR ao,A.V.,“Direct TrajectoryOptimizationandCostateEstimationviaanOrtho gonalCollocation Method,” JournalofGuidance,Control,andDynamics ,Vol.29,No.6, November-December2006,pp.1435–1440. [36]Huntington,G.T., AdvancementandAnalysisofaGaussPseudospectralTranscriptionforOptimalControl ,Ph.D.thesis,MassachusettsInstituteofTechnology, Cambridge,Massachusetts,2007. [37]Rao,A.V.,Benson,D.A.,Darby,C.L.,Patterson,M.A.,Sande rs, I.,andHuntington,G.T., User'sManualfor GPOPS Version2.2 http://gpops.sourceforge.net/,June2009. 168

PAGE 169

[38]Kameswaran,S.andBiegler,L.T.,“ConvergenceRatesfor DirectTranscription ofOptimalControlProblemsUsingCollocationatRadauPoint s,” Computational OptimizationandApplications ,Vol.41,No.1,2008,pp.81–126. [39]Garg,D.,Patterson,M.A.,Darby,C.L.,Francolin,C.,H untington,G.T., Hager,W.W.,andRao,A.V.,“DirectTrajectoryOptimization andCostate EstimationofFinite-HorizonandInnite-HorizonOptimalC ontrolProblemsvia aRadauPseudospectralMethod,” ComputationalOptimizationandApplications ,PublishedOnline:October6,2009.DOI10.1007/s10589-009 -9291-0. http://www.springerlink.com/content/n851q6n343p9k60 k/. [40]Garg,D.,Patterson,M.A.,Hager,W.W.,Rao,A.V.,Benson, D.A.,andHuntington, G.T.,“AUniedFrameworkfortheNumericalSolutionofOptim alControlProblems UsingPseudospectralMethods,” Automatica ,PublishedOnline,August2010. DOI:10.1016/j.automatica.2010.06.048. [41]Fahroo,F.andRoss,I.,“ASpectralPatchingMethodforD irectTrajectory Optimization,” JournalofAstronauticalSciences ,Vol.48,No.2-3,Apr.-Sept. 2000,pp.269–286. [42]Fahroo,F.andRoss,I.,“CostateEstimationbyaLegendr ePseudospectral Method,” JournalofGuidance,Control,andDynamics ,Vol.24,No.2,2001, pp.270–277. [43]Fahroo,F.andRoss,I.,“DirectTrajectoryOptimizati onbyaChebyshev PseudospectralMethod,” JournalofGuidance,Control,andDynamics ,Vol.25, No.1,2002,pp.160–166. [44]Fahroo,F.andRoss,I.,“PseudospectralMethodsforIn nite-HorizonNonlinear OptimalControlProblems,” JournalofGuidance,Control,andDynamics ,Vol.31, No.4,2008,pp.927–936. [45]Gong,Q.,Kang,W.,andRoss,I.M.,“APseudospectralMeth odfortheOptimal ControlofFeedbackLinearizableSystems,” IEEETransactionsonAutomatic Control ,Vol.51,No.7,July2006,pp.1115–1129. [46]Gong,Q.,Ross,I.M.,Kang,W.,andFahroo,F.,“Connecti onsBetweenthe CovectorMappingTheoremandConvergenceofPseudospectral Methods,” ComputationalOptimizationandApplications ,Vol.41,No.3,December2008, pp.307–335. [47]Gong,Q.,Fahroo,F.,andRoss,I.M.,“SpectralAlgorithm forPseudospectral MethodsinOptimalControl,” JournalofGuidance,Control,andDynamics ,Vol.31, No.3,May-June2008. [48]Ross,I.andFahroo,F.,“PseudospectralKnottingMethod sforSolvingOptimal ControlProblems,” JournalofGuidance,Control,andDynamics ,Vol.27,No.3, 2004,pp.397–405. 169

PAGE 170

[49]Ross,I.M.andFahroo,F.,“LegendrePseudospectralAppr oximationsof OptimalControlProblems,” LectureNotesinControlandInformationSciences Springer-Verlag,2003,pp.327–342. [50]Ross,I.andFahroo,F.,“AdvancesinPseusospectralMeth odsforOptimalControl,” AIAAGuidance,Navigation,andControlConference ,AIAAPaper2008-7309, Honolulu,Hawaii,August2008. [51]Ross,I.andFahroo,F.,“ConvergenceoftheCostatesDo esNotImply ConvergenceoftheControls,” JournalofGuidance,Control,andDynamics Vol.31,No.4,2008,pp.1492–1497. [52]Ross,I.M.andFahroo,F., User'sManualforDIDO2001 :AMATLABApplication forSolvingOptimalControlProblems ,DepartmentofAeronauticsandAstronautics, NavalPostgraduateSchool,TechnicalReportAAS-01-03,Monter ey,California, 2001. [53]Canuto,C.,Hussaini,M.Y.,Quarteroni,A.,andZang,T. A., SpectralMethodsin FluidDynamics ,Spinger-Verlag,Heidelberg,Germany,1988. [54]Fornberg,B., APracticalGuidetoPseudospectralMethods ,CambridgeUniversity Press,1998. [55]Trefethen,L.N., SpectralMethodsUsingMATLAB ,SIAMPress,Philadelphia, 2000. [56]Babuska,I.andSuri,M.,“The p and hp VersionoftheFiniteElementMethod,an Overview,” ComputerMethodsinAppliedMechanicsandEngineering ,Vol.80, 1990,pp.5–26. [57]Babuska,I.andSuri,M.,“The p and hp VersionoftheFiniteElementMethod, BasicPrinciplesandProperties,” SIAMReview ,Vol.36,1994,pp.578–632. [58]Gui,W.andBabuska,I.,“The h p ,and hp VersionsoftheFiniteElementMethodin 1Dimension.PartI.TheErrorAnalysisofthe p Version,” NumerischeMathematik Vol.49,1986,pp.577–612. [59]Gui,W.andBabuska,I.,“The h p ,and hp VersionsoftheFiniteElementMethodin 1Dimension.PartII.TheErrorAnalysisofthe h and h p Versions,” Numerische Mathematik ,Vol.49,1986,pp.613–657. [60]Gui,W.andBabuska,I.,“The h p ,and hp VersionsoftheFiniteElementMethodin 1Dimension.PartIII.TheAdaptive h p Version,” NumerischeMathematik ,Vol.49, 1986,pp.659–683. [61]Betts,J.T.andHuffman,W.P.,“MeshRenementinDirect TranscriptionMethods forOptimalControl,” OptimalControlApplicationsandMethods ,Vol.19,1998, pp.1–21. 170

PAGE 171

[62]Atkinson,K.E., AnIntroductiontoNumericalAnalysis ,JohnWileyandSons,New York,NY,1989. [63]Davis,P.J.andRabinowitz,P., MethodsofNumericalIntegration ,AcademicPress, NewYork,1984. [64]Soueres,P.andBoissonnat,J.D., RobotMotionPlanningandControl ,Springer, EnglewoodCliffs,NJ,1998. [65]Meditch,J.,“OntheProblemofOptimalThrustProgrammin gforaSoftLunar Landing,” IEEETransactiononAutomaticControl ,Vol.9,No.4,1964,pp.477–484. [66]Darby,C.L.,Hager,W.W.,andRao,A.V.,“Anhp-AdaptivePse udospectralMethod forSolvingOptimalControlProblems,” OptimalControlApplicationsandMethods August2010. [67]Darby,C.L.,Hager,W.W.,andRao,A.V.,“DirectTraject oryOptimizationbyA VariableLow-Orderhp-AdaptivePseudospectralMethod,” JournalofSpacecraft andRockets ,AcceptedforPublicationNovember2010(inpress). [68]Gill,P.E.,Murray,W.,andSaunders,M.A., User'sGuideforSNOPTVersion7: SoftwareforLargeScaleNonlinearProgramming ,February2006. [69]Betts,J.T., PracticalMethodsforOptimalControlandEstimationUsingNo nlinear Programming ,SIAMPress,Philadelphia,2010. [70]Rao,A.V.andMease,K.D.,“DichotomicBasisApproachtosol ving Hyper-SensitiveOptimalControlProblems,” Automatica ,Vol.35,No.4,April 1999,pp.633–642. [71]Rao,A.V.andMease,K.D.,“EigenvectorApproximateDichot omicBasisMethod forSolvingHyper-SensitiveoptimalControlProblems,” OptimalControlApplications andMethods ,Vol.21,No.1,January–February2000,pp.1–19. [72]Rao,A.V.,“ApplicationofaDichotomicBasisMethodtoPer formanceOptimization ofSupersonicAircraft,” JournalofGuidance,Control,andDynamics ,Vol.23,No.3, May–June2000,pp.570–573. [73]Rao,A.V.,“RiccatiDichotomicBasisMethodforsolvingH yper-Sensitiveoptimal ControlProblems,” JournalofGuidance,Control,andDynamics ,Vol.26,No.1, January–February2003,pp.185–189. [74]Bryson,A.E.,Desai,M.N.,andHoffman,W.C.,“Energy-State Approximationin PerformanceOptimizationofSupersonicAircraft,” AIAAJournalofAircraft ,Vol.6, No.6,1969,pp.481–488. [75]Seywald,H.andKumar,R.R.,“FiniteDifferenceSchemefo rAutomaticCostate Calculation,” JournalofGuidance,Control,andDynamics ,Vol.19,No.4,1996, pp.231–239. 171

PAGE 172

[76]London,H.S.,“ChangeofSatelliteOrbitPlanebyAerodyna micManeuvering,” IAS 29thAnnualMeeting ,InternationalAstronauticalSociety,NewYork,January196 1. [77]Walberg,G.D.,“ASurveyofAeroassistedOrbitTransfer, ” JournalofSpacecraft andRockets ,Vol.22,No.1,1985,pp.3–18. [78]Mease,K.D.,“AerossistedOrbitalTransfer:CurrentStat us,” TheJournalofthe AstronauticalSciences ,Vol.36,No.1/2,1988,pp.7–33. [79]Rao,A.V.,Tang,S.,andHallman,W.P.,“NumericalOptim izationStudyof Multiple-PassAeroassistedOrbitalTransfer,” OptimalControlApplicationsand Methods ,Vol.23,No.4,July–August2002,pp.215–238. [80]Darby,C.L.andRao,A.V.,“Minimum-FuelLow-EarthOrbit AeroassistedOrbital TransferofSmallSpacecraft,” JournalofSpacecraftandRockets ,Acceptedfor PublicationFebruary2011(inpress). [81]Vinh,N.-X.,Busemann,A.,andCulp,R., HypersonicandPlanetaryEntryFlight Mechanics ,ElsevierScience,NewYork,1981. [82]Clarke,K., PerformanceOptimizationofaCommonAeroVehicleUsingaLeg endre PseudospectralMethod ,Ph.D.thesis,DepartmentofAeronauticsandAstronautics, MassachusettsInstituteofTechnology,Cambridge,Massac husetts,2003. [83]Bate,R.R.,Mueller,D.D.,andWhite,J.E., FundamentalsofAstrodynamics ,Dover Publications,1971. [84]Detra,R.W.,Kemp,N.H.,andRiddell,F.R.,“Addendumto HeatTransferof SatelliteVehiclesRe-enteringtheAtmosphere,” JetPropulsion ,Vol.27,1957, pp.1256–1257. [85]Rump,S.M.,“INTLAB–INTervalLABoratory,” DevelopmentsinReliableComputing ,editedbyT.Csendes,KluwerAcademicPublishers,Dordrecht, 1999,pp. 77–104. [86]Tewari,A., AtmosphericandSpaceFlightDynamics ,Birkh ¨ auser,Boston,2007. 172

PAGE 173

BIOGRAPHICALSKETCH ChristopherDarbywasbornin1983inGainesville,Florida. Hereceivedhis BachelorofSciencedegreeinmechanicalengineeringinDecem ber2006,hisMasterof SciencedegreeinmechanicalengineeringinMay2010,andhis DoctorofPhilosophy inmechanicalengineeringinMay2011.Eachdegreewasearned attheUniversity ofFlorida.Hisresearchinterestsincludenumericalappro ximationstodifferential equations,optimalcontroltheory,andnumericalapproxim ationstothesolutionof optimalcontrolproblems. 173