Christopher L. Darby, William W. Hager, and Anil V. Rao

13 Pages·2011·0.52 MB·English

JOURNALOFSPACECRAFTANDROCKETS
Vol.48,No.3,May–June2011
Direct Trajectory Optimization Using a Variable Low-Order
Adaptive Pseudospectral Method
ChristopherL.Darby,∗WilliamW.Hager,†andAnilV.Rao‡
UniversityofFlorida,Gainesville,Florida32611
DOI:10.2514/1.52136
Avariable-orderadaptivepseudospectralmethodispresentedforsolvingoptimalcontrolproblems.Themethod
developedinthispaperadjustsboththemeshspacingandthedegreeofthepolynomialoneachmeshintervaluntila
specifiederrortoleranceissatisfied.Inregionsofrelativelyhighcurvature,convergenceisachievedbyrefiningthe
mesh, while in regions of relatively low curvature, convergence is achieved by increasing the degree of the
polynomial.Anefficientiterativemethodisthendescribedforaccuratelysolvingageneralnonlinearoptimalcontrol
problem.Usingfourexamples,theadaptivepseudospectralmethoddescribedinthispaperisshowntobemore
efficientthaneitheraglobalpseudospectralmethodorafixed-ordermethod.
Nomenclature x(cid:2)t(cid:4) = stateontimedomaint2(cid:5)t ;t (cid:7)
0 f
C = pathconstraintfunction x(cid:2)(cid:1)(cid:4) = stateontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7)
y(cid:2)t(cid:4) = componentofposition
D = N(cid:1)(cid:2)N(cid:3)1(cid:4)Radaupseudospectraldifferentiation
y(cid:2)t(cid:4) = generalvectorfunctionontimedomaint2(cid:5)t ;t (cid:7)
matrix 0 f
y(cid:2)(cid:1)(cid:4) = generalvectorfunctionontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7)
E = maximumabsolutesolutionerror
u = maximumallowablevalueofcontrol
F = magnitudeofdragforce,N max
d u = minimumallowablevalueofcontrol
F = magnitudeofgravityforce,N min
g U(cid:2)(cid:1)(cid:4) = controlapproximationontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7)
F = magnitudeofliftforce,N
l u(cid:2)t(cid:4) = controlontimedomaint2(cid:5)t ;t (cid:7)
F = magnitudeofthrustforce,N 0 f
t u(cid:2)t(cid:4) = controlontimedomaint2(cid:5)t ;t (cid:7)
g = integrandofcostfunctional 0 f
u(cid:2)(cid:1)(cid:4) = controlontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7)
h = altitude
I = numberofgridrefinementiterations v = speed
w = jthLegendre–Gauss–Radauquadratureweight
J = continuous-timecostfunctional j
(cid:2) = angleofattack,radanddeg
K = numberofmeshintervals
(cid:3) = flight-pathangle,radanddeg
L = magnitudeofliftforce,N
L(cid:2)(cid:1)(cid:4) = Lagrangepolynomialontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7) (cid:4) = accuracytolerance
(cid:5) = loadfactor
m = dimensionofcontinuous-timecontrol
(cid:6) = orientationangleorlongitude,rad
N = polynomialdegree
(cid:7) = trajectorycurvature
N = numberofapproximationpoints
a (cid:8) = Earthgravitationalparameter,m3=s2
N = numberofcollocationpoints
c (cid:9) = bankangle,rad
N = polynomialdegreeinmeshintervalk
k (cid:1) = timeontimeinterval(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7)
N = numberofnonzeroconstraintJacobianentries
z (cid:1) = Mayercost
n = dimensionofcontinuous-timestate
(cid:10) = latitude,radanddeg
r = geocentricradius,m
(cid:1) = boundaryconditionfunction
R = equatorialradiusofEarth,m
e = azimuthangle,radanddeg
s = dimensionofpathconstraintfunction
T = CPUtime,s
t = timeontimeintervalt2(cid:5)t ;t (cid:7)
t = terminaltime 0 f I. Introduction
f
t = initialtime OVER the past two decades, direct collocation methods have
0
t(cid:8) = controlswitchtime becomepopularinthenumericalsolutionofnonlinearoptimal
s
X(cid:2)(cid:1)(cid:4) = stateapproximationontimedomain(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7) control problems. In a direct collocation method, the state is
x(cid:2)t(cid:4) = componentofposition approximatedusingasetoftrial(basis)functionsandthedynamics
are collocated at specified set of points in the time interval. Most
commonly, direct collocation methods for optimal control are
PresentedinanabbreviatedformatasPaper2010-8272atthe2010AIAA/ employed using so-called h methods where a fixed low-degree
AASAstrodynamicsSpecialistConference,Toronto,Ontario,Canada,2–5
polynomial(e.g.,seconddegreeorthirddegree)stateapproximation
August 2010; received 23 August 2010; revision received 18 November
2010;acceptedforpublication19November2010.Copyright©2010byAnil isusedandtheproblemisdividedintomanyintervals.Convergence
V. Rao, William W. Hager, and Christopher L. Darby. Published by the of the numerical discretization is then achieved by increasing the
AmericanInstituteofAeronauticsandAstronautics,Inc.,withpermission. number of mesh intervals [1–3]. Grid refinement techniques are
Copiesofthispapermaybemadeforpersonalorinternaluse,oncondition commonlyusedtoobtainaspecifiedsolutionaccuracybyincreasing
thatthecopierpaythe$10.00per-copyfeetotheCopyrightClearanceCenter, thenumberofmeshintervalsinregionsofthetrajectorywherethe
Inc.,222RosewoodDrive,Danvers,MA01923;includethecode0022-4650/ errors are largest. Excellent examples of h methods for solving
11and$10.00incorrespondencewiththeCCC. optimalcontrolproblemsaregivenin[2–6].
∗Ph.D.Student.DepartmentofMechanicalandAerospaceEngineering;
In recent years, pseudospectral methods for solving optimal
cdarby@ufl.edu.
†Professor.DepartmentofMathematics;hager@ufl.edu. control problems have increased in popularity [7–19]. In a
‡Assistant Professor. Department of Mechanical and Aerospace pseudospectralmethod,thecollocationpointsarebasedonaccurate
Engineering; anilvrao@ufl.edu. Associate Fellow AIAA (Corresponding quadraturerulesandthebasisfunctionsaretypicallyChebyshevor
Author). Lagrangepolynomials.Incontrasttoanhmethod,apseudospectral
433
434 DARBY,HAGER,ANDRAO
methodistypicallyemployedasapmethodwhereasinglemesh startswithacoarsemeshanduseslow-degreepolynomialsineach
interval is used, and convergence is achieved by increasing the meshinterval.Thedegreeofthepolynomialwithinameshintervalis
degree,p,ofthepolynomial.Forproblemswherethesolutionsare thenincreasedonlyifthesolutioninthatmeshintervalissufficiently
infinitelysmoothandwellbehaved,apseudospectralmethodhasa smooth. Consequently, in the algorithm developed in this paper
simplestructureandconvergesatanexponentialrate[20–22].The more closely resembles an h method because the degree of the
most well-developed pseudospectral methods are the Gauss approximatingpolynomialinameshintervaliskeptrelativelysmall.
pseudospectral method [12,13], the Radau pseudospectral method Themethodof[33]maybepreferableoverthemethodofthispaper
(RPM)[17–19],andtheLobattopseudospectralmethod[7]. forproblemsforwhichthesolutionsareknownaprioritobesmooth
Whilepseudospectralmethodshavetypicallybeenappliedasp withtheexceptionofafewisolatedpoints,whilethemethodofthis
methods [7,8,12,13,18,19], relying on convergence using global paperispreferableforproblemsforwhichthesolutionstructureis
polynomialshasseverallimitations.First,evenforsmoothproblems, notknownapriori.
an accurate approximation may require an unreasonably large- Wealsocompareourmethodto[3],wherelow-degreeRunge–
degreeglobalpolynomial.Inaddition,theconvergencerateofap Kuttaintervalsareused.In[3],themeshisrefinedgloballybasedon
method on a problem for which the formulation or solution is theintegralofthecurvature.Inourhpscheme,ontheotherhand,the
nonsmooth may be extremely slow, resulting in a poor approxi- refinementtechniqueof[3]isappliedlocallytoeachmeshinterval
mation even when a high-degree polynomial is used. A second where refinement is deemed necessary. As a result, mesh interval
limitation of a p method is that the use of a high-degree global changesaremorelocalinthemethodofthispaperascomparedwith
polynomialresultsinanonlinearprogrammingproblem(NLP)for themethodof[3].Furthermore,thealgorithmofthispaperallowsfor
whichthedensitygrowsquicklyasafunctionofthenumberofNLP adifferentdegreepolynomialapproximationineachmeshinterval.
variables.ThisgrowthintheNLPdensityisduetothedenseblocks Finally,thegrowthrateofthemeshinouralgorithmisafunctionof
that arise from the global pseudospectral differentiation matrix. the error in the current solution and the desired user prescribed
Thus,apmethodbecomescomputationallyintractableforproblems accuracy.
whenanexcessivelylarge-degreepolynomialisrequired.Whileh Thispaperisorganizedasfollows.InSec.II,wedefinetheBolza
methodsarecomputationally moretractablethanpmethods,they optimal control problem of interest in this paper. In Sec. III, we
mayrequirealargenumberofmeshintervalsinordertoachievean formulate a multiple-interval discretization of the Bolza optimal
acceptableaccuracybecauseexponentialconvergenceislostusing control problem using the RPM. In Sec. IV, we provide two
anhmethod.Moreover,whentheNLPissufficientlylarge,itmaybe motivatingexamplesfortheuseofanhp-adaptivepseudospectral
difficulttocomputeasolution. method. In Sec. V, we present our variable-order hp-adaptive
To simultaneously improve accuracy and computational pseudospectralmethod.InSec.VI,weprovidefourexamplesofthe
efficiency using pseudospectral methods, we combine the best method using optimal control problems of varying complexity.
featuresofbothanhmethodandapmethodtoformaso-calledhp Finally,inSecs.VIIandVIII,weprovideadiscussionoftheresults
method.Asitsnameimplies,anhpmethodisonewherethenumber andconclusions.
of mesh intervals, the mesh interval widths, and the polynomial
degree in each mesh interval is determined algorithmically.
II. BolzaOptimalControlProblem
Previously,hpmethodshavebeendevelopedinthecontextoffinite
elementsinmechanicsandspectralmethodsinfluiddynamics.In Considerthefollowingfairlygeneraloptimalcontrolproblemin
particular,[23–27]describethemathematicalpropertiesofh,p,and Bolzaform.Minimizethecostfunctional
hp methods for finite elements. Galvao et al. [28] showed the Z
applicationofanhp-adaptiveleast-squaresspectralelementmethod J(cid:9)(cid:1)(cid:5)x(cid:2)t (cid:4);t ;x(cid:2)t (cid:4);t (cid:7)(cid:3) tfg(cid:5)x(cid:2)t(cid:4);u(cid:2)t(cid:4);t(cid:7)dt (1)
0 0 0 f
(LS-SEM) for solving hyperbolic partial differential equations. t0
Heinrichs [29] developed an adaptive spectral least-squares
collocation scheme for Burgers’s equation. Dorao and Jakobsen subjecttothedynamicconstraints
[30] showed how to apply of an hp-adaptive LS-SEM to the dx
populationbalanceequation,whileDoraoetal.[31]developedan (cid:9)f(cid:2)x(cid:2)t(cid:4);u(cid:2)t(cid:4);t(cid:4) (2)
dt
hp-adaptive spectral element solver for reactor modeling. An
overview of hp-adaptive spectral element methods for solving theinequalitypathconstraints
problemsincomputationalfluiddynamicscanbefoundin[32].
In this paper, we develop a new hp-adaptive pseudospectral C(cid:5)x(cid:2)t(cid:4);u(cid:2)t(cid:4);t(cid:7)(cid:10)0 (3)
method for solving optimal control problems. In the method
developed in this paper, the accuracy of the solution is improved
andtheboundaryconditions(i.e.,theeventconstraints)
either by increasing the degree of the polynomial within a mesh
interval or by refining the mesh. The decision to increase the (cid:10)(cid:5)x(cid:2)t (cid:4);t ;x(cid:2)t (cid:4);t (cid:7)(cid:9)0 (4)
polynomial degree or refine the mesh is based on the relative 0 0 0 f
curvatureineachmeshinterval.Iftheratiobetweenthemaximum where x(cid:2)t(cid:4)2Rn is the state, u(cid:2)t(cid:4)2Rm is the control,
curvature and the average curvature is sufficiently large on an C(cid:5)x(cid:2)t(cid:4);u(cid:2)t(cid:4);t(cid:7)2Rs represents the control and state constraints,
interval,thentheaccuracyisincreasedbyrefiningthemeshandby andtistime.Forconvenience,allquantitieswillbetreatedasrow
usinglow-degreepolynomialsinthenewlycreatedmeshintervals. vectors;thatis,ify(cid:2)t(cid:4)2Rq,theny(cid:2)t(cid:4)(cid:11)(cid:5)y (cid:2)t(cid:4); (cid:12)(cid:12)(cid:12); y (cid:2)t(cid:4)(cid:7).
1 q
Otherwise, the accuracy is improved by increasing the degree of Suppose now that the time interval t2(cid:5)t ;t (cid:7) is divided into a
approximating polynomial within a mesh interval. The method is 0 f
meshconsistingofKmeshintervals[t ,t ],k(cid:9)1;... ;K,where
demonstratedonfourexamplesofvaryingcomplexityandisfound k(cid:6)1 k
(t ;...;t )arethemeshpoints;themeshpointshavetheproperty
tobeaviablemethodforefficientlyandaccuratelysolvingcomplex 0 K
thatt <t <t <(cid:12)(cid:12)(cid:12)<t (cid:9)t .Ineachmeshintervalt2(cid:5)t ;t (cid:7),
optimalcontrolproblems. 0 1 2 K f k(cid:6)1 k
wemakethechangeofvariables
Itisnotedthat[33]alsodevelopsanhp-adaptivepseudospectral
methodforsolvingoptimalcontrolproblems.Themethodof[33] 2t(cid:6)(cid:2)t (cid:3)t (cid:4)
(cid:1)(cid:9) k k(cid:6)1 ; (cid:2)t <t (cid:4) (5)
startswithaglobalapproximation.Theaccuracyisthenimprovedby t (cid:6)t k(cid:6)1 k
k k(cid:6)1
increasingthedegreeofthepolynomialwithinameshintervaluntilit
is deemed necessary to refine the mesh. The mesh is refined only Using the change of variables given in Eq. (5), the interval t2
when exponential convergence is lost within an existing mesh (cid:5)t ;t (cid:7)istransformedto(cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7).Moreover,wehave
k(cid:6)1 k
interval.Asaresult,themethodof[33]morecloselyresemblesap
methodbecauseitcanresultinrelativelylarge-degreepolynomials d(cid:1) (cid:9) 2 ; (cid:2)k(cid:9)1;...;K(cid:4) (6)
withinameshinterval.Ontheotherhand,themethodofthispaper dt t (cid:6)t
k k(cid:6)1
DARBY,HAGER,ANDRAO 435
Next,letx(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)andu(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)bethestateandcontrolinthekthmesh NXk(cid:3)1 t (cid:6)t
interval as a function of (cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7). Using Eq. (6), the Bolza X(cid:2)k(cid:4)D(cid:2)k(cid:4)(cid:6) k k(cid:6)1f(cid:2)X(cid:2)k(cid:4);U(cid:2)k(cid:4);(cid:1)(cid:2)k(cid:4);t ;t (cid:4)(cid:9)0;
optimalcontrolproblemofEqs.(1–4)canbewrittenasfollows.First, j(cid:9)1 j ij 2 i i i k(cid:6)1 k
thecostfunctionalofEq.(1)canbewrittenas
(cid:2)i(cid:9)1;... ;N (cid:4) (14)
k
J(cid:9)(cid:1)(cid:5)x(cid:2)1(cid:4)(cid:2)(cid:6)1(cid:4);t0;x(cid:2)K(cid:4)(cid:2)(cid:3)1(cid:4);tf(cid:7) where
(cid:3)XK tk(cid:6)2tk(cid:6)1Z (cid:3)1g(cid:5)x(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);u(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);(cid:1);tk(cid:6)1;tk(cid:7)d(cid:1) (7) D(cid:2)ijk(cid:4)(cid:9)L_(cid:2)jk(cid:4)(cid:2)(cid:1)i(cid:2)k(cid:4)(cid:4); (cid:2)i(cid:9)1;...;Nk(cid:4) (cid:2)j(cid:9)1;... ;Nk(cid:3)1(cid:4)
k(cid:9)1 (cid:6)1 (15)
Next,thedynamicsofEq.(2)andthepathconstraintsofEq.(3)are is the N (cid:1)(cid:2)N (cid:3)1(cid:4) Radau pseudospectral differentiation matrix
k k
givenintermsof(cid:1)inmeshintervalk2(cid:5)1;...;K(cid:7),respectively,as [18]inthekthmeshinterval.Furthermore,thepathconstraintsof
Eq.(3)inmeshintervalk2(cid:5)1;... ;K(cid:7)areenforcedattheN LGR
dx(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4) (cid:9)tk(cid:6)tk(cid:6)1f(cid:5)x(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);u(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);(cid:1);t ;t (cid:7) (8) pointsas k
d(cid:1) 2 k(cid:6)1 k
C(cid:2)X(cid:2)k(cid:4);U(cid:2)k(cid:4);(cid:1)(cid:2)k(cid:4);t ;t (cid:4)(cid:10)0; (cid:2)i(cid:9)1;... ;N (cid:4) (16)
i i i k(cid:6)1 k k
0 (cid:13)C(cid:5)x(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);u(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4);(cid:1);tk(cid:6)1;tk(cid:7) (9) Finally,theeventconstraintsareapproximatedas
Furthermore,theboundaryconditionsofEq.(4)aregivenas (cid:1)(cid:2)X(cid:2)1(cid:4);t ;X(cid:2)K(cid:4) ;t (cid:4)(cid:9)0 (17)
1 0 NK(cid:3)1 K
(cid:1)(cid:5)x(cid:2)1(cid:4)(cid:2)(cid:6)1(cid:4);t ;x(cid:2)K(cid:4)(cid:2)(cid:3)1(cid:4);t (cid:7)(cid:9)0 (10)
0 f Continuityacrossthemeshpointsismaintainedbythecondition
Finally,itisrequiredthatthestatebecontinuousattheinterfaceof X(cid:2)k(cid:4) (cid:9)X(cid:2)k(cid:3)1(cid:4) (18)
meshintervals;thatis,x(cid:2)t(cid:6)k(cid:4)(cid:9)x(cid:2)t(cid:3)k(cid:4),k(cid:9)1;... ;K(cid:6)1. Nk(cid:3)1 1
TheNLPthatarisesfromtheRadaupseudospectralapproximationis
thentominimizethecostfunctionofEq.(12)subjecttothealgebraic
III. RadauPseudospectralMethod
constraintsofEqs.(14–18).Inourcomputerimplementationofthe
The hp method in this paper is developed by discretizing the NLP,weusethesamevariableforbothX(cid:2)k(cid:4) andX(cid:2)k(cid:3)1(cid:4).Hence,the
multiple-intervalformoftheBolzaoptimalcontrolproblem,givenin Nk(cid:3)1 1
constraint (18) can be eliminated since it is explicitly taken into
Sec.II,usingthepreviouslydevelopedRPM,asdescribedin[18]
account.
WhilewechoosetheRPMtodiscretizetheoptimalcontrolproblem,
ItisusefultoobserveafewkeypropertiesoftheRPM.First,when
with only slight modifications, the hp-adaptive method described
used as a single-interval p method (see Fig. 1a), the state is
this paper can be developed using other pseudospectral methods
approximated attheN collocation pointsplustheterminal point,
(e.g., the Gauss [12,13,15] or the Lobatto [7] pseudospectral k
method). An advantage of using the Radau scheme is that the
continuity conditions x(cid:2)t(cid:6)(cid:4)(cid:9)x(cid:2)t(cid:3)(cid:4) across mesh points are 10
k k
particularlyeasytoimplement.
9
In the RPM, the state of the continuous Bolza problem is
approximatedwithineachmeshintervalk2(cid:5)1;... ;K(cid:7)as 8
7
NXk(cid:3)1 NYk(cid:3)1(cid:1)(cid:6)(cid:1)(cid:2)k(cid:4)
x(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:14)X(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9) X(cid:2)k(cid:4)L(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4); L(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9) l
j j j (cid:1) (cid:6)(cid:1)(cid:2)k(cid:4) 6
j(cid:9)1 l(cid:9)1 j l
l≠j 5
(11)
4
where (cid:1)2(cid:5)(cid:6)1;(cid:3)1(cid:7), L(cid:2)jk(cid:4)(cid:2)(cid:1)(cid:4), j(cid:9)1;...;Nk(cid:3)1, is a basis of 3
Lagrange polynomials, ((cid:1)(cid:2)k(cid:4);...;(cid:1)(cid:2)k(cid:4)) are the Legendre–Gauss–
1 Nk 2
Radau (LGR) [18,34] collocation points in mesh interval k, and -1 -0.5 0 0.5 1
(cid:1)(cid:2)k(cid:4) (cid:9)(cid:3)1isanoncollocatedpoint.ThecostfunctionalofEq.(7)is
thNekn(cid:3)1approximatedusingamultiple-intervalLGRquadratureas a) Nc and Na for for a p–type radau collocation method
8
J(cid:14)(cid:1)(cid:2)X(cid:2)1(cid:4);t ;X(cid:2)K(cid:4) ;t (cid:4)
1 0 NK(cid:3)1 K
XK XNk (cid:1)t (cid:6)t (cid:2) 7
(cid:3) k k(cid:6)1 w(cid:2)k(cid:4)g(cid:2)X(cid:2)k(cid:4);U(cid:2)k(cid:4);(cid:1)(cid:2)k(cid:4);t ;t (cid:4) (12)
2 j j j j k(cid:6)1 k 6
k(cid:9)1 j(cid:9)1
whereU(cid:2)k(cid:4),i(cid:9)1;...;N ,aretheapproximationsofthecontrolat 5
i k
theN LGRpointsinthekthmeshinterval,X(cid:2)1(cid:4)istheapproximation
k 1 4
ofx(cid:2)t (cid:4),andX(cid:2)K(cid:4) istheapproximationofx(cid:2)t (cid:4).Differentiating
X(cid:2)k(cid:4)(cid:2)(cid:1)0(cid:4)inEq.(1N1K)(cid:3)w1ithrespectto(cid:1),weobtain f 3
dX(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4) (cid:11)X_(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9)NXk(cid:3)1X(cid:2)k(cid:4)L_(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4) (13) 2-1 -0.5 0 0.5 1
d(cid:1) j j
j(cid:9)1
b) N and N for an h–type radau collocation method
c a
CollocatingthedynamicsattheNkLGRpointsusingEq.(13),we Fig.1 NumberofcollocationpointsNcandnumberofapproximation
have pointsNaforap-andh-Radaucollocationmethod.
436 DARBY,HAGER,ANDRAO
whilethecontrolisapproximatedatonlytheN collocationpoints. andthecontrolinequalityconstraint
k
Foramulti-intervaldiscretizationwhereaRadauschemeisusedfor
eachintervalinthemesh,thecontrolisapproximatedateverypoint umin(cid:10)u(cid:10)umax (26)
wherethestateisapproximatedwiththesingleexceptionofthefinal
where u (cid:9)0, u (cid:9)3, g(cid:9)1:5, and t is free. The optimal
pointt(cid:9)t (seeFig.1b).Finally,wecommentonthefactthatthe min max f
f solutiontothisexampleis
RPMdescribedinthissectioncanbewrittenequivalentlyineither
differentialorimplicitintegralform[18,19].Wechoosetousethe
(cid:2)h(cid:8)(cid:2)t(cid:4);v(cid:8)(cid:2)t(cid:4);u(cid:8)(cid:2)t(cid:4)(cid:4)
differentialform(thatis,theformdescribedinthispaper)becauseit
hascomputationaladvantagesovertheintegralform.See[18,19]for ((cid:2)(cid:6)3t2(cid:3)v t(cid:3)h ;(cid:6)3t(cid:3)v ;0(cid:4); t<t(cid:8)
further details on the equivalence between the differential and (cid:9) h 4 0 0 2 0 i s
3t2(cid:3)(cid:2)(cid:6)3t(cid:8)(cid:3)v (cid:4)t(cid:3)3t(cid:8)2(cid:3)h ;3t(cid:3)(cid:2)(cid:6)3t(cid:8)(cid:3)v (cid:4);3 ; t>t(cid:8)
integralsformsoftheRPM. 4 s 0 2s 0 2 s 0 s
(27)
IV. MotivatingExamples wheret(cid:8)is
s
In this section, we consider examples that motivate the
t(cid:8) v
developmentofanhp-adaptivepseudospectralmethodforsolving t(cid:8)(cid:9) f(cid:3) 0 (28)
optimalcontrolproblems.Thefirstexamplehasasmoothsolution, s 2 3
whilethesecondexamplehasanonsmoothsolution.Thegoalofthe
with
motivatingexamplesistoprovideinsightintowhenanhmethodora
pmethodmaybemoreappropriate. q(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)(cid:3)
t(cid:8)(cid:9)2v (cid:3)4 1v2(cid:3)3h (29)
f 3 0 3 2 0 2 0
A. Example1:ProblemwithaSmoothSolution FortheboundaryconditionsgiveninEq.(25),wehavet(cid:8)(cid:9)1:4154
s
Considerthefollowingoptimalcontrolproblem[35].Minimize andt(cid:8)(cid:9)4:1641.AsseeninEq.(27),theoptimalcontrolu(cid:8)(cid:2)t(cid:4)for
f
thecostfunctional thisexampleisbang–bang.Asaresult,thecontrolisdiscontinuous
andthestatecomponentsh(cid:2)t(cid:4)andv(cid:2)t(cid:4)arenonsmooth[inthiscase,
J(cid:9)t (19)
f h(cid:2)t(cid:4)ispiecewisequadraticandv(cid:2)t(cid:4)ispiecewiselinear,asseenin
Eq.(27)].
subjecttothedynamicconstraints
TheconvergenceratesforthisproblemusingaRadaupmethod
x_(cid:9)cos(cid:2)(cid:6)(cid:4); y_(cid:9)sin(cid:2)(cid:6)(cid:4); (cid:6)_(cid:9)u (20) and an equally spaced Radau h-2 method are now analyzed.
ExaminingFig.2c,itisseenthatthepmethodandtheh-2method
converge erratically and much more slowly than was the case for
andtheboundaryconditions
example1.Althoughtheconvergenceofbothmethodsappearstobe
x(cid:2)0(cid:4)(cid:9)0; y(cid:2)0(cid:4)(cid:9)0; (cid:6)(cid:2)0(cid:4)(cid:9)(cid:6)(cid:11) x(cid:2)t (cid:4)(cid:9)0; erratic, it turns out that the asymptotic convergence rate for both
f
methods is log linear. In other words, the error is bounded by
y(cid:2)tf(cid:4)(cid:9)0; (cid:6)(cid:2)tf(cid:4)(cid:9)(cid:11) (21) the product of a constant with the inverse of the degree of the
approximating polynomial (in the case of the p method) or the
ThesolutiontotheoptimalcontrolproblemofEqs.(19–21)is productofaconstantandtheinverseofthenumberofmeshintervals
(inthecaseofthehmethod).Moreover,theconvergencerateforthe
(cid:5)x(cid:8)(cid:2)t(cid:4);y(cid:8)(cid:2)t(cid:4);(cid:6)(cid:8)(cid:2)t(cid:4);u(cid:8)(cid:2)t(cid:4)(cid:7)(cid:9)(cid:5)(cid:6)sin(cid:2)t(cid:4);(cid:6)1(cid:3)cos(cid:2)t(cid:4);t(cid:6)(cid:11);1(cid:7) (22)
hmethodinFig.2cisstronglyinfluencedbythedistancebetweenthe
where t(cid:8)(cid:9)2(cid:11). It is seen that the solution to the optimal control switchtimeofthecontinuous-timeoptimalsolutionandtheclosest
problemfgiveninEqs.(19–21)issmooth.Wewillapproximatethis meshpoint,andthisdistancedependsinacomplicatedwayonthe
numberofmeshintervals.
solution using a Radau pseudospectral p method, and we will
To see more clearly the underlying log linear asymptotic
comparetheerrortothatofauniformlyspacedsecond-degreeRadau
convergencerateofthepandh-2methods,considerthefollowing
h method. Figures 2a and 2b show the maximum state error as a
equivalent fixed-time formulation of the optimal control problem
function of the number of collocation points and the number of
giveninEqs.(23–26).Minimizethecostfunctional
nonzero constraint Jacobian entries, respectively, for both the p
pmreotbhloedmanisdstmheosoethco,nthde-dpegmreeethhomdectohnovde.rBgeescaautsaentheexspoolnuetinotniatlortahteis, J(cid:9)h2(cid:2)tf(cid:4)(cid:3)v2(cid:2)tf(cid:4)(cid:3)Z tfudt (30)
whilethesecond-degreehmethodconvergesatasignificantlyslower 0
rateascomparedwiththepmethod.Consequently,inthisexample, subjecttothedynamicconstraintsgiveninEq.(24),theboundary
convergenceisachievedmuchmorerapidlyusingapmethod.In conditions
particular,foranygivendesiredsolutionaccuracy,theoverallsizeof
theNLPresultingfromthepmethodismuchsmallerthanthesize (cid:5)h(cid:2)0(cid:4);v(cid:2)0(cid:4)(cid:7)(cid:9)(cid:2)h ;v (cid:4)(cid:9)(cid:2)10;(cid:6)2(cid:4) (31)
0 0
from the h method, thus making the NLP of the p method much
easiertosolve. thecontrolinequalityconstraintgiveninEq.(26),andtf(cid:9)t(cid:8)f[where
t(cid:8)isgivenbyEq.(29)].Themodifiedoptimalcontrolproblemwas
f
B. Example2:ProblemwithaNonsmoothSolution solvedfor2j (j(cid:9)1;... ;8)collocationpointswhere,inthiscase,
eachmeshhastwicethenumberofintervalsasthepreviousmesh.As
Consider the following optimal control problem of a soft lunar
aresult,thedistancebetweentheclosestmeshpointandtheoptimal
landing,asgivenin[36].Minimizethecostfunctional
switchpointt(cid:8)decreasestozeromonotonicallyasafunctionofj.In
s
Z tf contrast,whenallpossibleequallyspacedmeshesareconsidered(as
J(cid:9) udt (23) showninFig.2c),thedistancebetweentheclosestgridpointandthe
0 switchingpointdoesnotdecreasemonotonicallyasafunctionofthe
subjecttothedynamicconstraints numberofmeshintervals.Figure2dshowsthebase-10logarithmof
themaximumabsoluteerrorasafunctionofthebase-2logarithmof
h_(cid:9)v; v_(cid:9)(cid:6)g(cid:3)u (24) thenumberofcollocationpoints.AsseeninFig.2d,theasymptotic
convergence rate of both methods is log linear. Because the
theboundaryconditions asymptotic error in both the p and h-2 method are similar, their
efficiency depends on the relative time it takes to solve discrete
(cid:5)h(cid:2)0(cid:4);v(cid:2)0(cid:4);h(cid:2)t (cid:4);v(cid:2)t (cid:4)(cid:7)(cid:9)(cid:2)h ;v ;h ;v (cid:4)(cid:9)(cid:2)10;(cid:6)2;0;0(cid:4) (25)
f f 0 0 f f problems of similar dimensions. As seen in Figs. 2e and 2f, the
DARBY,HAGER,ANDRAO 437
2 2
0 0
-2 -2
-4 -4
-6 -6
-8 -8
0 20 40 60 80 100 0 500 1000 1500 2000 2500 3000
a) log E vs. N using a p–Method and a second– degree b) log E vs. N using a p–Method and a second– degree
10 c 10 z
h–Method motivating example 1 h–Method motivating example 1
0.5 0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1.5 -2
-2 -2.5
0 20 40 60 80 0 2 4 6 8
c) log E vs. N using a p–Method and a second– degree d) log E vs. log N using a p–Method and a second– degree
10 c 10 2 c
h–Method for the optimal control problem of eqs. (23)– h–Method for the alternate form of motivating example 2
(26)
0.5 0.4
0
0.3
-0.5 (s)
e
m
Ti 0.2
-1 U
P
C
0.1
-1.5
-2 0
0 5000 10000 15000 0 20 40 60 80
e) log E vs. N using a p–Method and a second– degree e) T vs. N using a p–Method and a second– degree
10 z c
h–Method for the motivating example 2 h–Method for motivating example 2
Fig.2 h-methodandp-methoderrorsformotivatingexamples1and2.
sparsityandthesolutiontimeoftheh-2methodgrowmuchmore attheappropriatelocations.Thesolutionofageneraloptimalcontrol
slowly than the sparsity and solution time of the p method. problem will require a combination of choosing the appropriate
Consequently, for this example, where the optimal solution has a polynomial degree and the appropriate locations of the mesh
bang–bangcontrol,theh-2methodismuchmorecomputationally intervals.Inaregionwherethesolutioniswellbehaved,itisexpected
efficient. thatconvergencewillbeachievedfastestbyincreasingthedegreeof
the approximating polynomial. In a region where the solution is
either not smooth or is not accurately approximated using a
C. DiscussionofMotivatingExamples reasonablylow-degreepolynomial,refiningthemeshatthelocations
Theresultsofthemotivatingexamplespointtothefactthat,in ofthenonsmoothnessandusingapolynomialofanappropriate(but
general,thecharacterofthesolutiontoanoptimalcontrolproblem relativelylow)degreeineachmeshintervalwouldappeartobethe
willchangeindifferentregionsofthesolution.Inexample1,itwas best choice. Because the behavior of the solution to an arbitrary
seenthatrapidconvergencewasachievedsimplybyincreasingthe optimalcontrolproblemisnotknownapriori,itisusefultoallowfor
degreeofthepolynomial.Ontheotherhand,inexample2,itwas boththepolynomialdegreetobedifferentineachmeshintervaland
seenthatwhenthedegreeoftheapproximatingpolynomialisfixed, forthenumberandlocationsofthemeshintervalstobedetermined.
therateofconvergencecanbeimprovedbyplacingthemeshpoints Combiningtheselasttwofeaturesleadstoaso-calledhpmethod.In
438 DARBY,HAGER,ANDRAO
this paper, we develop an hp method for solving optimal control areinaregionwherethejthpathconstraintisactive,b(cid:2)k(cid:4)mayhave
lj
problems using pseudospectral methods where higher-order positivevalues.Ifanyofthesepositivevaluesexceeds(cid:4),thepath
polynomialsareusedinregionswherethesolutionissmoothand constraintisnotsatisfiedtothedesiredaccuracytolerance.
largernumbersoflow-degreepolynomialsareusedinregionswhere Ourerrorestimationprocessfocusesontheerrorintheconstraints.
thesolutionisnotwellapproximatedbypolynomialsofareasonable Theaccuracyinthefirst-orderoptimalityconditionsinthecontrol
degree. problemisimplicitlytakenintoaccountwhenwespecifyanerror
toleranceforthefirst-orderoptimalityconditionsintheNLPsolver.
V. Variable-Orderhp-AdaptiveAlgorithm
B. IncreasingPolynomialDegreeorRefiningMesh
A. AssessmentofApproximationErrorinaMeshInterval
Iftheaccuracyinmeshintervalkneedstobeimproved,thefirst
Consider again a K-interval Bolza optimal control problem, as
definedinSec.II,onthetimeintervalt2(cid:5)t ;t (cid:7).Let[t ,t ]bethe stepistodetermineifthemeshshouldberefinedorifthedegreeof
0 f k(cid:6)1 k theapproximationshouldbeincreased.SupposethatX(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)isthe
kthintervalandsupposethatthesolutioniscollocatedatN pointson m
k componentofthestateapproximationinthekthmeshintervalthat
thisinterval.Thegoalofthehp-adaptivealgorithmistoimprovethe
accuracy of the solution in a computationally efficient manner by correspondstothemaximumvalueofeithera(cid:2)lik(cid:4).Thecurvatureofthe
determiningifaparticularmeshintervalinthecurrentmeshhasmeta mthcomponentofthestateinmeshintervalkisgivenby
specified accuracy tolerance. If a mesh interval has not met the
jX(cid:3)(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)j
accuracytolerance,thenthenumberanddistributionofcollocation (cid:7)(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9) m (36)
pointsneedstobemodifiedeitherbyincreasingthedegreeofthe j(cid:5)1(cid:3)X_(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)2(cid:7)3=2j
m
approximatingpolynomialinthemeshintervaland/orrefiningthe
mesh. Let (cid:7)(cid:2)k(cid:4) and (cid:7)(cid:2)(cid:2)k(cid:4) be the maximum and mean value of (cid:7)(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4),
max
Let (cid:4) be an accuracy tolerance for the discretized differential- respectively,where(cid:7)(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)iscomputedusingEq.(36).Furthermore,
algebraic constraints (that is, the discretized dynamic and path letr betheratioofthemaximumtothemeancurvature:
k
constraints). The kth mesh interval is considered to be within the
accuracy tolerance if the maximum violation in the differential- r (cid:9)(cid:7)(cid:2)mka(cid:4)x (37)
algebraicequationson[tk(cid:6)1,tk]isbelow(cid:4).Unlikethestate,which k (cid:7)(cid:2)(cid:2)k(cid:4)
hasapolynomialapproximation,asgiveninEq.(11),thecontrolhas
nosuchdefinedapproximation.Whileanyinterpolatingfunctionfor If rk<rmax, where rmax>0 is a user-defined parameter, the
whichthevaluematchestheapproximationforthecontrolattheN curvatureisconsidereduniforminthismeshintervalandalarger-
k
LGRcollocationpointsinmeshintervalkisacceptable,asfarasthe degreepolynomialisusedtoobtainabetterapproximationtomesh
sdoelgurteioenLoafgrthanegNeLpPoliysncoomnciaelrnisedu,siendthtoisappappreorx,itmheatfeoltlhoewcionngtNroklthin- itnhteerrevsatlokf.tIhferkm>esrhmianx,tetrhveanltahnedretheeximstessahliasrgreeficnuerdv.aturerelativeto
meshintervals1;...;K(cid:6)1:
1. DeterminationofNewPolynomialDegreeWithinaMeshInterval
U(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9)NXk(cid:3)1U(cid:2)k(cid:4)L^(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4); L^(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9)NYk(cid:3)1(cid:1)(cid:6)(cid:1)i(cid:2)k(cid:4) (32) SupposetheresultofSec.V.Bisthatthepolynomialdegreewithin
i i j (cid:1) (cid:6)(cid:1)(cid:2)k(cid:4) meshintervalkshouldbeincreased.LetMdenotetheinitialdegree
i(cid:9)1 ii(cid:9)≠1j j i ofthepolynomialinthisinterval.ThedegreeNkofthepolynomialin
meshintervalkisgivenbytheformula
Becausethefinaltimet isnotcollocated,thefollowing(N (cid:6)1)th-
f k
degreeLagrangepolynomialisusedtoapproximatethecontrolinthe N (cid:9)M(cid:3)ceil(cid:5)log (cid:2)e(cid:2)k(cid:4) =(cid:4)(cid:4)(cid:7)(cid:3)A (38)
finalmeshintervalK: k 10 max
whereA>0isanarbitraryintegerconstantdescribedinSec.V.D,
U(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9)XNk U(cid:2)k(cid:4)L~(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4); L~(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)(cid:9)YNk (cid:1)(cid:6)(cid:1)i(cid:2)k(cid:4) (33) e(cid:2)mka(cid:4)xisthemaximumofa(cid:2)lik(cid:4)andb(cid:2)ljk(cid:4)overl,i,andj,andceilisthe
i i j (cid:1) (cid:6)(cid:1)(cid:2)k(cid:4) operatorthatroundstothenexthighestinteger.
i(cid:9)1 i(cid:9)1 j i
i≠j It is seen from Eq. (38) that the increase in the degree of the
polynomialinthekthmeshintervalisbasedontheratiobetweenthe
ItisnotedinEq.(32)that,fork2(cid:5)1;K(cid:6)1(cid:7),thesupportpointsofthe maximum error and the specified error tolerance. For a smooth
LagrangebasisL^(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)aretheNkLGRpointsplusthenoncollocated function,apseudospectralmethodexhibitexponentialconvergence.
point (cid:1)(cid:9)(cid:3)1, whereas in Eq. (33), the support points of the Hence,thegrowthinthepolynomialdegreeisrelatedtothelogofthe
LagrangebasisL~(cid:2)k(cid:4)(cid:2)(cid:1)(cid:4)aretheN LGRpoints. error.Itishopedthataunitincreaseinthepolynomialdegreewill
k
To estimate the error in the differential-algebraic equations, we improvetheerrorbyoneorderofmagnitude.
insertthecurrentstateandcontrolpolynomialsintheequationsand
evaluate them at L points (cid:2)t(cid:2)(cid:2)1k(cid:4);... ;t(cid:2)(cid:2)Lk(cid:4)(cid:4)2(cid:5)tk(cid:6)1;tk(cid:7) in each mesh 2. DeterminationofNumberandPlacementofNewMeshPoints
intervalk2(cid:5)1;...;K(cid:7):
SupposetheresultofSec.V.Bisthatthekthmeshintervalshould
berefined.Thefollowingprocedureis thenusedtodetermine the
jx_(cid:2)k(cid:4)(cid:2)t(cid:2)(cid:2)k(cid:4)(cid:4)(cid:6)f(cid:2)k(cid:4)(cid:2)X(cid:2)k(cid:4);U(cid:2)k(cid:4);t(cid:2)(cid:2)k(cid:4)(cid:4)j(cid:9)a(cid:2)k(cid:4) (34)
i l i l l l li locations of the new mesh points. First, the new number of mesh
intervals,denotedn ,iscomputedas
k
C(cid:2)jk(cid:4)(cid:2)X(cid:2)lk(cid:4);U(cid:2)lk(cid:4);t(cid:2)(cid:2)lk(cid:4)(cid:4)(cid:9)b(cid:2)ljk(cid:4) (35) nk(cid:9)ceil(cid:5)Blog10(cid:2)e(cid:2)mka(cid:4)x=(cid:4)(cid:4)(cid:7) (39)
where1(cid:10)l(cid:10)L,1(cid:10)i(cid:10)n,and1(cid:10)j(cid:10)s.Ifeveryelementofa(cid:2)k(cid:4) whereBisanarbitraryintegerconstant,describedinSec.V.D.The
li
andb(cid:2)k(cid:4)islessthan(cid:4)onthecurrentmeshinterval,thenthespecified locationsofthenewmeshpointsaredeterminedusingtheintegralof
lj
errortolerance(cid:4)hasbeensatisfiedonthisinterval.Ifanyelementof acurvaturedensityfunctioninamannersimilartothatgivenin[3].
a(cid:2)k(cid:4) or b(cid:2)k(cid:4) is greater than (cid:4), then either the mesh interval will be Specifically,let(cid:12)(cid:2)(cid:1)(cid:4)bethedensityfunction[3]givenby
li lj
refined or the degree of the polynomial on this interval will be (cid:12)(cid:2)(cid:1)(cid:4)(cid:9)c(cid:7)(cid:2)(cid:1)(cid:4)1=3 (40)
increased.
ItisnotedthatEq.(35)measuresthequalitywithwhichthepath wherecisaconstantchosen,sothat
constraintofEq.(3)issatisfiedinmeshintervalkatapointl.Ifweare
Z
inaregionofthetrajectorywherethejthpathconstraintisinactive, (cid:3)1
(cid:12)(cid:2)(cid:13)(cid:4)d(cid:13)(cid:9)1
bl(cid:2)jk(cid:4)willbenegativeand,therefore,willalwaysbelessthan(cid:4).Ifwe (cid:6)1
DARBY,HAGER,ANDRAO 439
LetF(cid:2)(cid:1)(cid:4)bethecumulativedistributionfunctiongivenby (whereBcontrolsthegrowthinthenumberofmeshintervals).Itis
importanttounderstandthatatradeoffexistsbetweenusinglargeand
Z
(cid:1) smallvaluesofeitherAorB.IfeitherAorBissufficientlylarge,the
F(cid:2)(cid:1)(cid:4)(cid:9) (cid:12)(cid:2)(cid:13)(cid:4)d(cid:13) (41)
algorithm will use fewer iterations to converge to an acceptable
(cid:6)1
solution,butthesizeofthemeshorthenumberofcollocationpoints
Thenknewmeshpointsarechosensothat maygrowquicklybetweeniterations.IfeitherAorBissmall,the
meshwillgrowmuchmoreslowly,butthealgorithmmayrequire
i(cid:6)1
F(cid:2)(cid:1)(cid:4)(cid:9) ; 1(cid:10)i(cid:10)n (cid:3)1 manymoreiterationstoachieveanacceptablesolution.
i nk k Theparameterrmaxisusedtodeterminewhethertosubdividethe
currentintervalorincreasethedegreeofapproximatingpolynomial.
In[3],itisshownthatdistributingthemeshpointsinthiswayis,in Forr <1,thestrategytorefinethemeshwillalwaysoccurand,
max
some sense, optimal for the piecewise linear approximation of a simply, an h method will be used. If r (cid:9)1, then for each
curvewhentheerrorismeasuredintheL1norm.Ifnk(cid:9)1,thenno inaccuratemeshinterval,thefirstrefinemenmtaxattemptwillalwaysbe
subintervalsarecreated.Therefore,theminimumvaluefornkshould toincreasethedegreeofapproximatingpolynomialintheinterval.
beatleasttwo. Byusingafinitevalueofr greaterthanone,thealgorithmwill
max
eithersubdividethecurrentmeshintervalorincreasethedegreeof
C. hp-AdaptiveAlgorithm polynomialapproximation.
Theparameter(cid:4)specifiestheaccuracyofthedynamicconstraints
We now describe our hp-adaptive algorithm that uses the
andpathconstraintsoftheNLPapproximationbetweencollocation
components described in Secs. V.A, V.B, V.B.1, and V.B.2. The
points.Specifically,thesmallerthevalueof(cid:4),thegreatertherequired
algorithmstartsbyformingacoarsemesh,usingpolynomialsofa
fixeddegreeoneachinterval,andsolvingtheresultingNLP.Based accuracyinthedynamics.Inthisstudy,weassumethattheproblem
ontherulesgivennext,foreachintervalinthemesh,weeitherrefine hasbeenscaledsothat(cid:4)islessthanone.Ifthecontrolproblemis
themeshorincreasethedegreeofthepolynomials.Whenwerefine poorlyscaled,thenthecorrespondingNLPmaynotbesolvable.
themesh,weuseafixeddegreeMforthepolynomialsoneachnew
meshinterval.Whenweincreasethedegreeofthepolynomials,we VI. Examples
compute the new degree using the formula given in Eq. (38). We
Thehp-adaptivemethodofSec.Visnowdemonstratedonfour
continuetherefinementprocessuntilthespecifiederrortoleranceis examples.Allexamplesweresolvedusingamodifiedversionofthe
satisfied.Inmoredetail,therefinementprocessproceedsasfollows:
open-sourcesoftwareGPOPS[16]withtheNLPsolverSNOPT[37],
1)SolvetheNLPusingthecurrentmesh.
usinga2.5GHZCore2DuoMacbookProrunningMacOS-X10.5.8
Begin:Fork(cid:9)1;...;K,
withMATLABR2009b.Notethatalgorithmicdecisionsdependon
2)Ife(cid:2)mka(cid:4)x (cid:10)(cid:4),thencontinue(proceedtonextk). thechoiceoftheNLPsolver,andifadifferentsolverwasusedthat
3) If either rk(cid:13)rmax or Nk>M, then refine the kth mesh better exploits the sparsity of the discretized optimal control
intervalintonksubintervals,wherenkisgivenbyEq.(39).Set problem,thenalgorithmicdecisionscouldchange.AllCPUtimes
thedegreeofthepolynomialsoneachofthesubintervalstobe shownincludeonlythetimerequiredtosolvetheNLPonallgrids
Mandproceedtothenextk. anddonotincludethesetuptimerequiredtogeneratethegridforthe
4)Otherwise,setthedegreeofthepolynomialsonthekth subsequent mesh iteration. For all examples, the following values
subintervaltobeNk,whichisgiveninEq.(38). werechosenfortheparametersdescribedinSec.V:A(cid:9)1,B(cid:9)2,
End:Fork(cid:9)1;... ;K. L(cid:9)30,andr (cid:9)2.ToensurethattheNLPissolvedtothesame
max
5)ReturntoStep1. accuracyastheoptimalcontrolproblem,theSNOPTfeasibilityand
optimalitytolerancesaresettoavaluesmallerthan(cid:4).Next,noupper
D. Remarks limitwasplacedonthedegreeofthepolynomial,ascomputedby
Inthisalgorithm,wealternatebetweenameshrefinementstepin Eq.(38),butalimitofnk(cid:10)3Bwasplacedonthegrowthofthemesh,
as computed by Eq. (39). To obtain consistent results in all
which the degree of the polynomials are set to M on each of the
refined intervals and a step where we increase the degree of the computations,thesamestoppingcriterionisusedforallthemethods.
Intheparticularcaseofapmethod,ifasolutionisunacceptable,the
polynomial according to Eq. (38). The goal of the hp-adaptive
degreeoftheglobalpolynomialisincreasedarbitrarilyby10untila
algorithmofthispaperistouseasmanylow-degreemeshintervalsas
solutionisachieved.Furthermore,inordertoensureafinesampling
possibleandtoonlyuselarger-degreepolynomialsinregionswhere
betweenthecollocationpoints,avalueL(cid:9)1000ischosenforap
thesolutionissmooth.Evenforsubdomainsofthetrajectorywhere
method. If a solution is not attained upon reaching a 200th-deg
larger-degree polynomials may be appropriate, the degree of the
polynomialapproximation,thenitisdecidedthatthepmethodwas
polynomial approximation is still small when compared with a p
unabletoobtainasolution.Theparameterr controlsthedecision
method. The reason for controlling the degree of the polynomial max
toeitherrefinethemeshorincreasethedegreeofthepolynomial.For
approximation is threefold. First, as discussed, the number of
thehmethod,wesetr <1toforcethealgorithmtoperformonly
nonzeroentriesintheNLPconstraintJacobiangrowsasthesquareof max
meshrefinement(thatis,thenumberofmeshintervalscanchange,
thenumberofcollocationpointsineachinterval.Thus,ifthedegree
butthedegreeofthepolynomialinameshintervalcannotchange
of the approximating polynomial in a mesh interval becomes too
fromitsinitialvalue).Theterminologyh-xdenotesanhmethodof
large,theNLPbecomesincreasinglydenseandtheNLPsolverhasto
degreexineverymeshinterval,whilehp(cid:6)xdenotesanhpmethod
work much harder to compute a solution. Second, the difference
withaminimumpolynomialdegreeofxineachinterval.Finally,for
betweenthenumberofnonzeroconstraintJacobianentriesforanh
examples 1 and 2, the maximum absolute error shown is the
method and hp method is small, even when the total number of
maximum absolute difference between the NLP solution and the
collocationpointsislarge,providedthatthemaximumpolynomial
exact solution at the state approximation points (that is, the LGR
degreeinanymeshintervalissmall.Finally,duetotheexponential
pointsoneachmeshintervalplustheterminalpoint).
convergenceofapseudospectralmethodinameshintervalwherethe
solution is smooth, it should be possible to obtain an accurate
solutionusingarelativelysmallnumberofcollocationpointsinthat A. Example1
meshinterval.Foralloftheexamplesstudiedinthispaper,aseventh- ConsideragainthemotivatingexamplegivenbyEqs.(19–21).It
degreepolynomialwasthelargestevercomputedbythealgorithm wasseeninSec.IV.Athatthemostrapidconvergencewasachieved
describedinSec.V.C. usingapmethod.Supposenowthatweconsidertheperformanceof
Next,thealgorithmusesthearbitrarilychosenparametersA,as thep,h-2,andhp-2methodsonthisexample.Inthecaseofap
describedinEq.(38)(whereAcontrolsthegrowthofthenumberof method,theinitialmeshwastwocollocationpoints,whilefortheh
collocationpointsinameshinterval)andB,asdescribedinEq.(39) methodandthehpmethod,theinitialmeshwas10meshintervals
440 DARBY,HAGER,ANDRAO
Table1 Summaryofaccuracyandspeedforexample1,usingvarious solution,thelocationofmaximumcurvatureinthestatewilloccurat
collocationstrategiesandaccuracytolerances thesametimeasthediscontinuityinthecontrol.Itisseenthat,except
for the second grid iteration, the maximum curvature in the state
(cid:4) E Method T,s N K I N
c z convergestothetruelocationofthediscontinuityinthecontrolasthe
10(cid:6)3 2:40(cid:1)10(cid:6)7 p 0.015 12 1 2 650 mesh is refined. Table 2 shows the performance of the algorithm
5:96(cid:1)10(cid:6)3 hp-2 0.0091 20 10 1 482 usingp,h-2,andhp-2methods.First,itisinterestingtonotethatthe
5:96(cid:1)10(cid:6)3 h-2 0.0091 20 10 1 482
pmethodwasneverabletoproducesolutionsforthegivenvaluesof
10(cid:6)4 2:40(cid:1)10(cid:6)7 p 0.015 12 1 2 650
(cid:4).Whenthepmethodwasforcedtostopat200collocationpoints,it
4:70(cid:1)10(cid:6)6 hp-2 0.023 40 10 2 1,202
6:66(cid:1)10(cid:6)4 h-2 0.040 64 32 3 1,538 resultedinaCPUtimeof4.84sand82,002nonzeroJacobianentries.
10(cid:6)5 2:40(cid:1)10(cid:6)7 p 0.015 12 1 2 650 Inanycase,theCPUtimeandnumberofnonzeroNLPconstraint
1:21(cid:1)10(cid:6)7 hp-2 0.030 50 10 2 1,652 Jacobian entries using the p method are significantly greater than
1:39(cid:1)10(cid:6)5 h-2 0.089 160 80 3 3,842 thoseobtainedusingeithertheh-2methodorhp-2method.Using
10(cid:6)6 2:40(cid:1)10(cid:6)7 p 0.015 12 1 2 650 eitherthehmethodorthehp-2methodwith(cid:4)(cid:9)10(cid:6)1,theaccuracy
2:59(cid:1)10(cid:6)9 hp-2 0.030 60 10 2 2,162 tolerance was satisfied on the initial mesh. For (cid:4)(cid:9)(cid:2)10(cid:6)2;
5:22(cid:1)10(cid:6)7 h-2 0.47 480 240 4 11,522 10(cid:6)3;10(cid:6)4(cid:4),thecomputationalperformanceofthehp-2andtheh-2
10(cid:6)7 2:40(cid:1)10(cid:6)7 p 0.015 12 1 2 650
methodswassimilar.
4:87(cid:1)10(cid:6)11 hp-2 0.037 70 10 2 2,732
1:98(cid:1)10(cid:6)8 h-2 4.63 1,568 784 5 37,634
C. Example3
Consider the following optimal control problem, which is a
variationoftheminimumtimetoclimbofasupersonicaircraft[38–
withtwocollocationpointsperinterval.Table1providesasummary
40].Minimizethecostfunctional
oftheresultsobtainedforthisexample.Whenusingthepmethod,
thealgorithmincreasedthepolynomialdegreeonceinordertomeet
J(cid:9)t (42)
the accuracy tolerance of all values of (cid:4). In the case of the hp-2 f
method,itseenfor(cid:4)(cid:9)10(cid:6)3thatthealgorithmdidnotiterate(thatis,
subjecttothedynamicconstraints
the solution on the first grid met the accuracy tolerance).
Furthermore, for (cid:4)<10(cid:6)3, the mesh was never refined; only the F (cid:6)F F
h_(cid:9)vsin(cid:3); v_(cid:9) t d(cid:6)F sin(cid:3); (cid:3)_(cid:9) g(cid:2)(cid:5)(cid:6)cos(cid:3)(cid:4)
degree of the polynomial was increased in order to satisfy the m g v
accuracytolerance.Moreover,itisseenfor(cid:4)(cid:9)10(cid:6)3and10(cid:6)4that
(43)
computational performance of the h-2 and hp-2 methods was
similar.Finally,weobservefor(cid:4)(cid:9)10(cid:6)3and10(cid:6)4thatthenumberof theboundaryconditions
nonzeroentriesintheNLPconstraintJacobianusingeithertheh-2
andhp-2methodswasessentiallythesame. h(cid:2)0(cid:4)(cid:9)0 m; h(cid:2)t (cid:4)(cid:9)19995 m; v(cid:2)0(cid:4)(cid:9)129:31m=s
f
Supposenowthatweanalyzethecomputationalperformanceof
v(cid:2)t (cid:4)(cid:9)295:09m=s; (cid:3)(cid:2)0(cid:4)(cid:9)0 deg; (cid:3)(cid:2)t (cid:4)(cid:9)0 deg (44)
thedifferentmethodsfor(cid:4)<10(cid:6)3downto(cid:4)(cid:9)10(cid:6)7.First,itisseen f f
that for all values of (cid:4), from 10(cid:6)3 to 10(cid:6)7, that computational
efficiencyusingthepmethodandthehp-2methodareessentially andtheinequalitypathconstraints
the same, with the p method being slightly more efficient for the h(cid:13)0; (cid:3) (cid:10)(cid:3) (45)
higheraccuracytolerances.Next,itisseenthat,as(cid:4)isdecreasedfrom max
10(cid:6)5to10(cid:6)7,thehp-2methodismorecomputationallyefficientthan Toforcetheflight-pathanglepathconstraint(cid:3) (cid:10)(cid:3) tobeactiveon
max
theh-2method.Wenoteagainthat,inthecaseofthehp-2method, theoptimalsolution,wechoose(cid:3) (cid:9)45 deg.Additionaldetails
max
only the degree of the polynomial was increased in each mesh forthevehicleusedinthisproblemcanbefoundin[39,40].
interval (that is, mesh refinement was never performed using the This example was solved using the p, h-x, and hp-x methods,
hp-2method),andtheaccuracytolerancewassatisfiedonthesecond where the h-x and hp-x methods were initialized with a mesh
grid iteration. The fact that mesh refinement was never employed consistingof10uniformmeshintervalswithxcollocationpointsin
usingthehp-2methodisconsistentwiththefactthatthesolutionto each mesh interval and the p solutions were initialized using 20
thisproblemissmooth.Interestingly,themaximumvalueofrk on collocationpointsonthetimeinterval.Figures4a–4dshowtheflight-
anymeshintervalofthesolutionobtainedontheinitialgridwas1.32, pathangleoniterations(1,2,4,and8)forthehp-4methodusing
wellbelowthechosenvalueofrmax(cid:9)2.Bycontrast,usingtheh-2 (cid:4)(cid:9)10(cid:6)3. It is seen that the initial mesh does not accurately
method with (cid:4)<10(cid:6)3, the algorithm chose to refine the mesh approximatetheoptimalsolutionineitherthetakeoffregion(from
substantially. As a result, for higher accuracy tolerances, the h-2 t(cid:9)0tot(cid:14)10s)orintheregionnearpathconstraint(thatis,from
method was less computationally efficient when compared with t(cid:14)20 stot(cid:14)35 s).Asthemeshrefinementprogresses,however,it
hp-2method,whilethepmethodshowedminimalimprovementin isseenthatthesolutioninthesetworegionsimprovesdramatically,
computationalefficiencyoverthehp-2method. leadingtomuchhigheraccuracyonthefinalmeshwherethemajority
ofmeshpointsandcollocationpointsareplacedintheregionfrom
t(cid:9)0tot(cid:14)40 s.Finally,intheregiont2(cid:5)(cid:14)50;t (cid:7),itisobserved
f
B. Example2 that the mesh does not change from the second mesh iteration
ConsideragaintheexamplegiveninEqs.(23–26).Supposewe onwards.
solve this problem using an hp-2 method with an initial mesh Table3showsthecomputationtimesandnumberofnonzeroNLP
consistingof10uniform-widthmeshintervals.Figures3a–3dshow constraint Jacobian entries for this example using the p, h-2, and
theevolutionofthecontrolapproximationonvariousmesheswith hp-xmethods.Itisseenthatasolutionisobtainedusingapmethod
(cid:4)(cid:9)10(cid:6)3,wheretheinterpolationisperformedusingEqs.(32)and for(cid:4)(cid:9)(cid:2)10(cid:6)1;10(cid:6)2;10(cid:6)3(cid:4).Nosolutionwasobtainedfor(cid:4)(cid:9)10(cid:6)4,
(33).Onthefirstmesh,thediscontinuityinthecontrolisnotcaptured asthepmethodwasstoppedat200collocationpoints.Also,using
accurately. On the subsequent meshes, however, more collocation eithertheh-2methodorthehp-xmethods,itisseenthattheinitial
points are placed near the location of the discontinuity, thus mesh satisfies the accuracy tolerance for (cid:4)(cid:9)10(cid:6)1. For
improving the accuracy with which the switch in the control is (cid:4)(cid:9)(cid:2)10(cid:6)1;10(cid:6)2(cid:4),theperformanceoftheh-2,hp-x,andpmethods
captured.Ontheinitialgrid,r (cid:9)2:07fortheintervalcontainingthe aresimilar.For(cid:4)(cid:9)10(cid:6)3,itisseenthatthepmethodissignificantly
k
discontinuity,forcingsubdivisionofthisinterval.Figure3eshows lesscomputationallyefficient,whiletheh-2andhp-xmethodshave
the difference between the time at which the maximum curvature similar computational performance. No p-method solution was
occurs,denotedtmax,andtheoptimalswitchtimeofthecontrolt(cid:8)asa foundfor(cid:4)(cid:9)10(cid:6)4.TheresultingCPUtimeandnonzeroJacobian
(cid:7) s
functionoftheiterationnumber.Itisexpectedthat,foranaccurate entrieswere1100.4sand123,602,respectively.When(cid:4)(cid:9)10(cid:6)4,a
DARBY,HAGER,ANDRAO 441
3 3
2.5 Collocation 2.5 Collocation
Interpolation Interpolation
2 Optimal 2 Optimal
u 1.5 u 1.5
1 1
0.5 0.5
0 0
0 1 2 3 4 5 0 1 2 3 4 5
t t
a) Mesh 1 b) Mesh 2
3 3
2.5 Collocation 2.5 Collocation
Interpolation Interpolation
2 Optimal 2 Optimal
u 1.5 u 1.5
1 1
0.5 0.5
0 0
0 1 2 3 4 5 0 1 2 3 4 5
t t
c) Mesh 4 d) Mesh 8
25
ds) 20
n
o
c
millise 15
(
10
5
0
0 2 4 6 8
MeshIterationNumber
e) Difference between the time at which maximum
curvature Occurs, tmax, and the optimal switch time
K
in the control, t*, as a function of mesh iteration number
s
Fig.3 Controlforexample2onmeshesa)1,b)2,c)4,andd)8,usinganhp-2methodwithanaccuracytoleranceof(cid:2)(cid:1)10(cid:2)3,ande)t(cid:3)max(cid:2)ts(cid:3)vsmesh
iterationnumber.
significantimprovementisseenusinghigher-orderhp-xmethods,as
compared with the h-2 method. In comparing the h-2 and hp-2
methods, it is important to observe that allowing the flexibility to Table2 Summaryofaccuracyandspeedforexample2,usingvarious
increasethedegreeofthepolynomialfromitsminimumallowable collocationstrategiesandaccuracytolerance
value greatly reduces the total number of collocation points.
(cid:4) E Method T,s N K I N
Furthermore, the hp-4 method showed the best performance for c z
(cid:4)(cid:9)10(cid:6)4.Thus,whilethehp-4methodresultsinthedenserNLP,it 10(cid:6)1 —— p —— —— —— —— ——
alsoturnsoutthatthehp-4methodhasboththefewestnumberof 2:74(cid:1)10(cid:6)2 hp-2 0.097 20 10 1 282
collocation points and the fewest number of nonzero Jacobian 2:74(cid:1)10(cid:6)2 h-2 0.097 20 10 1 282
10(cid:6)2 —— p —— —— —— —— ——
entries.Thisexampleagaindemonstratestherelativeineffectiveness
4:96(cid:1)10(cid:6)3 hp-2 0.12 40 20 8 562
ofthepmethodovereithertheh-2methodoranhp-xmethodwhen
7:01(cid:1)10(cid:6)3 h-2 0.10 28 14 5 394
anactivepathconstraintispresent. 10(cid:6)3 —— p —— —— —— —— ——
3:20(cid:1)10(cid:6)3 hp-2 0.13 50 24 8 718
D. Example4 6:32(cid:1)10(cid:6)3 h-2 0.096 36 18 5 506
10(cid:6)4 —— p —— —— —— —— ——
Considerthefollowingoptimalcontrolproblem[1]ofmaximizing
6:97(cid:1)10(cid:6)4 hp-2 0.29 86 42 11 1222
the crossrange during the atmospheric entry of a reusable launch
7:20(cid:1)10(cid:6)4 h-2 0.20 72 36 8 1010
vehicle.Minimizethecostfunctional
442 DARBY,HAGER,ANDRAO
50 50
Collocation Collocation
40 40
Interpolation Interpolation
30 30
g) g)
de 20 de 20
( (
γ γ
10 10
0 0
-10 -10
0 50 100 150 200 0 50 100 150 200
t (s) t (s)
a) Mesh 1 b) Mesh 2
50 50
Collocation Collocation
40 Interpolation 40 Interpolation
30 30
g) g)
de 20 de 20
( (
γ γ
10 10
0 0
-10 -10
0 50 100 150 200 0 50 100 150 200
t (s) t (s)
c) Mesh 4 d) Mesh 8
Fig.4 Flight-pathangle(cid:4)(deg)vstforexample3onmeshiterations1,2,4,and8,usingthehp-4methodwithanaccuracytoleranceof(cid:2)(cid:1)10(cid:2)3.
J(cid:9)(cid:6)(cid:10)(cid:2)t (cid:4) (46) successfulforallvaluesof(cid:4).Furthermore,for(cid:4)(cid:9)(cid:2)10(cid:6)4;10(cid:6)5(cid:4),the
f
solutionsusinganyofthemethodswereessentiallythesame.Asthe
subjecttothedynamicconstraints accuracytoleranceistightened,however,theh-2methodrequiresa
very large number of collocation points to meet the accuracy
r_(cid:9)vsin(cid:3); (cid:6)_(cid:9)vcos(cid:3)sin ; (cid:10)_(cid:9)vcos(cid:3)cos tolerance. As a result, a significantly larger computation time is
rcos(cid:10) r requiredtosolvetheNLPusingtheh-2methodwhencomparedwith
F F cos(cid:9) (cid:1)F v(cid:2) usingthehp-xmethods.Inaddition,whilethepmethodresultedin
v_(cid:9)(cid:6) d(cid:6)F sin(cid:3); (cid:3)_(cid:9) l (cid:6) g(cid:6) cos(cid:3) solutionswiththefewestnumberofcollocationpoints,theNLPthat
m g mv v r
resultsfromtheuseofhigh-degreepolynomialsismuchmoredense
F sin(cid:9) vcos(cid:3)sin tan(cid:10) thantheNLPsarisingfromthehp-xandhmethods.Thisincreased
_(cid:9) l (cid:3) (47)
mvcos(cid:3) r NLP density in turn increases the computational time required to
solvetheproblem.Inparticular,itisseenfor(cid:4)(cid:9)10(cid:6)6thatthehp-x
andtheboundaryconditions
r(cid:2)0(cid:4)(cid:9)79248(cid:3)R m; r(cid:2)t (cid:4)(cid:9)24384(cid:3)R m Table3 Summaryofaccuracyandspeedforexample3,using
e f e variouscollocationstrategiesandaccuracytolerances
(cid:6)(cid:2)0(cid:4)(cid:9)0 deg; (cid:6)(cid:2)t (cid:4)(cid:9)free; (cid:10)(cid:2)0(cid:4)(cid:9)0 deg
f (cid:4) Method T,s N K I N
c z
(cid:10)(cid:2)t (cid:4)(cid:9)free; v(cid:2)0(cid:4)(cid:9)7803m=s; v(cid:2)t (cid:4)(cid:9)762m=s
f f 10(cid:6)1 p 0.20 30 1 2 3,424
(cid:3)(cid:2)0(cid:4)(cid:9)(cid:6)1 deg; (cid:3)(cid:2)t (cid:4)(cid:9)(cid:6)5 deg; (cid:2)0(cid:4)(cid:9)90 deg hp-4 0.16 40 10 1 1,202
f hp-3 0.083 30 10 1 812
(cid:2)t (cid:4)(cid:9)free (48) hp-2 0.064 20 10 1 482
f
h-2 0.064 20 10 1 482
Furtherdetailsofthisproblem,includingtheaerodynamicmodel, 10(cid:6)2 p 0.37 40 1 3 5,522
hp-4 0.32 48 12 3 1,442
canbefoundin[1].Itisnotedthat,unlikeeitheroftheprevioustwo
hp-3 0.22 35 11 3 977
exampleswhereeitherthecontrolwasdiscontinuousoraninequality
hp-2 0.51 38 17 10 962
path constraint was active, this problem has a much smoother
h-2 0.42 46 23 7 1,106
solution. 10(cid:6)3 p 21.45 110 1 10 38,282
This examplewas solved using the p, h-x, and hp-x methods, hp-4 1.59 95 23 8 2,903
wheretheh-xandhp-xmethodswereinitializedwithaninitialmesh hp-3 2.38 118 38 14 3,248
consistingof20uniformmeshintervalswithxcollocationpointsin hp-2 3.03 117 43 14 3,191
each mesh interval and the p solutions were initialized using 20 h-2 1.35 138 69 7 3,314
collocationpointsonthetimeinterval.Figures5a–5fshowthestate 10(cid:6)4 p —— —— —— —— ——
hp-4 7.33 192 46 10 5,906
obtainedusingthehpalgorithmofSec.V,wheretheinterpolationis
hp-3 8.09 234 72 10 6,608
performed using Eq. (11). Table 4 shows the computational
hp-2 29.14 245 95 20 6,617
performance using the p, h-x, and hp-x methods. Because the
h-2 18.88 334 167 13 8,018
optimal solution to this example is smooth, the p method was