Michael A. Patterson, William W. Hager, and Anil V. Rao

24 Pages·2014·1.43 MB·English

OPTIMALCONTROLAPPLICATIONSANDMETHODS
Optim.ControlAppl.Meth.(2014)
PublishedonlineinWileyOnlineLibrary(wileyonlinelibrary.com).DOI:10.1002/oca.2114
A ph mesh refinement method for optimal control
MichaelA.Patterson1,WilliamW.Hager2 andAnilV.Rao1,*,†
1DepartmentofMechanicalandAerospaceEngineering,UniversityofFlorida,Gainesville,FL32611,USA
2DepartmentofMathematics,UniversityofFlorida,Gainesville,FL32611,USA
SUMMARY
A mesh refinement method is described for solving a continuous-time optimal control problem using
collocationatLegendre–Gauss–Radaupoints.Themethodallowsforchangesinboththenumberofmesh
intervals and the degree of the approximating polynomial within a mesh interval. First, a relative error
estimate is derived based on the difference between the Lagrange polynomial approximation of the state
andaLegendre–Gauss–Radauquadratureintegrationofthedynamicswithinameshinterval.Thederived
relativeerrorestimateisthenusedtodecideifthedegreeoftheapproximatingpolynomialwithinamesh
shouldbeincreasedorifthemeshintervalshouldbedividedintosubintervals.Thedegreeoftheapprox-
imatingpolynomialwithinameshintervalisincreasedifthepolynomialdegreeestimatedbythemethod
remainsbelowamaximumallowabledegree.Otherwise,themeshintervalisdividedintosubintervals.The
process of refining the mesh is repeated until a specified relative error tolerance is met. Three examples
highlightvariousfeaturesofthemethodandshowthattheapproachismorecomputationallyefficientand
producessignificantlysmaller meshsizesforagivenaccuracytolerancewhencomparedwithfixed-order
methods.Copyright©2014JohnWiley&Sons,Ltd.
Received16September2013; Revised29January2014; Accepted30January2014
KEYWORDS: optimalcontrol;collocation;Gaussianquadrature;variable-order;meshrefinement
1. INTRODUCTION
Over the past two decades, direct collocation methods have become popular in the numerical
solutionofnonlinearoptimalcontrolproblems.Inadirectcollocationmethod,thestateandcontrol
arediscretizedatasetofappropriatelychosenpointsinthetimeintervalofinterest.Thecontinuous-
time optimal control problem is then transcribed to a finite-dimensional nonlinear programming
problem (NLP), and the NLP is solved using a well-known software [1, 2]. Originally, direct col-
location methods were developed as h methods (e.g., Euler or Runge–Kutta methods) where the
time interval is divided into a mesh and the state is approximated using the same fixed-degree
polynomialineachmeshinterval.Convergence inan hmethodisthenachievedbyincreasingthe
number and placement of the mesh points [3–5]. More recently, a great deal of research has been
carried out in the class of direct Gaussian quadrature orthogonal collocation methods [6–20]. In
a Gaussian quadrature collocation method, the state is typically approximated using a Lagrange
polynomialwherethesupportpointsoftheLagrangepolynomialarechosentobepointsassociated
withaGaussianquadrature.Originally,Gaussianquadraturecollocationmethodswereimplemented
asp methodsusingasingleinterval.Convergenceofthep methodwasthenachievedbyincreas-
ing the degree of the polynomial approximation. For problems whose solutions are smooth and
well-behaved, a Gaussian quadrature collocation method has a simple structure and converges at
anexponentialrate[21–23].Themostwell-developedGaussianquadraturemethodsarethosethat
*Correspondence to: Anil V. Rao, Department of Mechanical and Aerospace Engineering, University of Florida,
Gainesville,FL32611,USA.
†E-mail:anilvrao@ufl.edu
Copyright©2014JohnWiley&Sons,Ltd.
M.A.PATTERSON,W.W.HAGERANDA.V.RAO
employeitherLegendre–Gausspoints[9,13]orLegendre–Gauss–Radau(LGR)points[14,15,17]
orLegendre–Gauss–Lobattopoints[6].
While h methods have a long history and p methods have shown promise in certain types of
problems, both the h and p approaches have limitations. Specifically, achieving a desired accu-
racy tolerance may require an extremely fine mesh (in the case of an h method) or may require
the use of an unreasonably large-degree polynomial approximation (in the case of a p method).
Inordertoreducesignificantlythesizeofthefinite-dimensionalapproximation,andthusimprove
computational efficiency of solving the NLP, hp collocation methods have been developed. In an
hp method, both the number of mesh intervals and the degree of the approximating polynomial
within each mesh interval is allowed to vary. Originally, hp methods were developed as finite
elementmethodsforsolvingpartialdifferentialequations[24–28].Inthepastfewyears,theprob-
lem of developing hp methods for solving optimal control problems has been of interest [29, 30].
Thisrecentresearchhasshownthatconvergenceusinghp methodscanbeachievedwithasignifi-
cantlysmallerfinite-dimensionalapproximationthanwouldberequiredwhenusingeitheranhora
pmethod.
Motivated by the desire to improve computational efficiency when solving the NLP while pro-
vidinghighaccuracyinthediscreteapproximationoftheoptimalcontrolproblem,inthispaper,we
describe a new ph Gaussian quadrature collocation method for solving continuous-time nonlinear
optimalcontrolproblems.Here,wehavedeliberatelychangedtheorderfromhptophbecauseour
schemetriestoachievetheprescribederrortolerancebyincreasingthedegreeofthepolynomialsin
theapproximatingspace,andifthisfails,byincreasingthenumberofmeshintervals.Themethodis
dividedintothreeparts.First,anapproachisdevelopedforestimatingtherelativeerrorinthestate
within each mesh interval. This relative error estimate is obtained by comparing the original state
variable to a higher-order approximation of the state. This relative error estimate is used to deter-
mineifthedegreeofthepolynomialapproximationshouldbeincreasedorifmeshintervalsshould
be added. The polynomial degree is increased if it is estimated that the ensuing mesh requires a
polynomial degree that is less than a maximum allowable degree. Otherwise, the mesh is refined.
Thisprocessisrepeatedonaseriesofmeshesuntilaspecifiedaccuracytoleranceismet.Thedeci-
sion to increase the polynomial degree or refine the mesh is based on the ratio of the maximum
relativeerrorandtheaccuracytoleranceandisconsistentwiththeknownexponentialconvergence
ofaGaussianquadraturemethodforaproblemwhosesolutionissmooth.
Variousmeshrefinementmethodsemployingdirectcollocationmethodshavebeendescribedin
recent years [5, 29–31]. Reference [31] describes a method that employs a differentiation matrix
to attempt to identify switches, kinks, corners, and other discontinuities in the solution, and uses
Gaussianquadraturerulestogenerateameshthatisdenseneartheendpointsofthetimeinterval
of interest. Reference [5] employs a density function and attempts to generate a fixed-order mesh
on which to solve the problem. References [32] and [33] (and the references therein) describe a
dual weighted residual method for mesh refinement and goal-oriented model reduction. The dual
weighted residual method uses estimates of a dual multiplier together with local estimates of the
residualstoadaptivelyrefineameshandcontroltheerrorinproblemsgovernedbypartialdifferen-
tial equations. References [29] and [30] describe hp adaptive methods where the error estimate is
basedonthedifferencebetweenanapproximationofthetimederivativeofthestateandtheright-
handsideofthedynamicsmidwaybetweenthecollocationpoints.Itisnotedthattheapproachof
References [29] and [30] creates a great deal of noise in the error estimate, thereby making these
approachescomputationallyintractablewhenahigh-accuracysolutionisdesired.Furthermore,the
errorestimateofReferences[29]and[30]doesnottakeadvantageoftheexponentialconvergence
rate of a Gaussian quadrature collocation method. Finally, in Reference [3], an error estimate is
developedbyintegratingthedifferencebetweenaninterpolationofthetimederivativeofthestate
and the right-hand side of the dynamics. The error estimate developed in Reference [3] is predi-
catedontheuseofafixed-ordermethod(e.g.,trapezoid,Hermite–Simpson,andRunge–Kutta)and
computesalow-orderapproximationoftheintegraloftheaforementioneddifference.Ontheother
hand,inthispaper,anestimateoftheerroronagivenmeshisobtainedbyvaryingthedegreeofthe
polynomialintheGaussianquadratureapproximation.
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
AphMESHREFINEMENTMETHODFOROPTIMALCONTROL
The method of this paper is fundamentally different from any of the previously developed
methods due to the fact that it takes advantage of the exponential convergence properties of a
Gaussianquadrature.Specifically,inthediscretizationusedinthispaper,thestateisapproximated
usingapiecewisepolynomial,andthedynamicsarecollocatedattheLGRquadraturepointsineach
interval. This leads to a finite-dimensional mathematical programming problem that is solved to
obtainestimatesforthecontrolandthestateatthecollocationpoints.Arelativeerrorestimateisthen
derived that uses the difference between an interpolated value of the state and an LGR quadrature
approximation to the integral of the state dynamics. This relative error estimate remains computa-
tionallytractablewhenahigh-accuracysolutionisdesiredandreducessignificantlythenumberof
collocationpointsrequiredtomeetaspecifiedaccuracytolerancewhencomparedwiththemethods
ofReference[29]or[30].Furthermore,differentfromtheerrorestimatedevelopedinReference[3],
theerrorestimateinthisresearchutilizestheLGRquadraturepointsinboththeinterpolationofthe
stateandtheintegrationofthestatealongtheinterpolatedsolution.Consequently,theapproachof
thisresearchenablestheuseofaGaussianquadraturemethodtoestimatetheerrors,therebyachiev-
ingconvergenceusingfewercollocationpointsthanmaybenecessaryusinganhmethodwiththe
error estimate of Reference [3]. It is noted, however, that in the case of an h LGR method (i.e., a
methodwhoseorderisthesameinallmeshintervals),theerrorestimateinthispaperissimilarto
theestimatederivedinReference[3]withtheexceptionthattheapproachofthispaperstillemploys
ahigher-orderintegrationmethodthantheapproachofReference[3].
While the mesh refinement method presented in this paper is based on only the state error and
seemstoignoretheerrorinthecostateaswellasignorethecontrolminimumprinciple,inearlier
work(suchasReferences[14–16]),acloseconnectionhasbeenestablishedbetweenthefirst-order
optimality conditions for the mathematical program and the continuous first-order optimality con-
ditions (i.e., the system dynamics, the costate equation, and the minimum principle). Specifically,
itisobservedthatatasolutionofthemathematicalprogrammingproblem,theminimumprinciple
holdsatthecollocationpoints.Thus,bysolvingthemathematicalprogram,theminimumprinciple
issatisfiedexactly.AshasbeenshowninReferences[14–16],thediscretizationofthestateequation
leads to an induced discretization for the costate equation that is expressed in terms of the mul-
tipliers of the mathematical programming problem. Hence, when the mathematical programming
problem is solved, in essence a two point boundary-value problem is solved in which the control
has been eliminated through the minimum principle. In this paper, we develop a mesh refinement
methodthatisbasedentirelyontheestimationoftheerrorinthestate.Whiletheerrorintheinduced
costateapproximationcouldbemonitoredusingthesameapproachasisusedtomonitorthestate
error,onefindsthattheerrorsconnectedwiththesystemdynamicsandtheinducedcostateequation
aretightlycoupledinthattheinducedcostateerrorislargeduringintervalswherethestateerroris
large.Asaresult,thebenefitofmonitoringtheerrorinthestateandtheinducedcostateismarginal
whencomparedwithmonitoringtheerrorinonlythestate.Thus,itismoreefficientandsimplerto
focusentirelyontheerrorinthestateequation.Theeffectivenessofourmeshrefinementstrategy
isstudiedonthreeexamplesthathavedifferentfeaturesintheoptimalsolutionand,thus,exercise
differentbenefitsofthenewphapproach.Itisfoundthatthemeshrefinementmethoddevelopedin
thispaperisaneffectiveyetsimplewaytogeneratemeshesandtoreducecomputationtimeswhen
comparedwithfixed-orderhmethods.
Thesignificanceofthisresearchisthreefold.First,wedevelopasystematicwaytoestimatethe
errorintheRadaudiscretizationofanoptimalcontrolproblemusingthefactthatRadaucollocation
isaGaussianquadratureintegrationmethod.Second,basedonthederivedestimateoftheerror,we
provideasimpleyeteffectivephmeshrefinementmethodthatallowsboththedegreeoftheapprox-
imating polynomial and the number of mesh intervals to vary. Third, we show on three nontrivial
examplesthatthephmeshrefinementmethoddevelopedinthispaperismorecomputationallyeffi-
cientandproducessmallermeshesforagivenaccuracytolerancewhencomparedwithtraditional
fixed-orderhmethods.
Thispaperisorganizedasfollows.InSection2,weprovideamotivationforournewphmethod.
InSection3,westatethecontinuous-timeBolzaoptimalcontrolproblem.InSection4,westatethe
integralformofthephLGRintegrationmethod[14–16]thatisusedasthebasisforthephmesh
refinement method developed in this paper. In Section 5, we develop the error estimate and our
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
M.A.PATTERSON,W.W.HAGERANDA.V.RAO
newph-adaptivemeshrefinementmethod.InSection6,weapplythemethodofSection5tothree
examplesthathighlightdifferentfeaturesofthemethod.InSection7,wedescribethekeyfeatures
of our approach and compare our method to recently developed hp-adaptive methods. Finally, in
Section8,weprovideconclusionsonourwork.
2. MOTIVATIONFORNEWph-ADAPTIVECOLLOCATIONMETHOD
In order to motivate the development of our new ph-adaptive Gaussian quadrature collocation
method,considerthefollowingtwofirst-orderdifferentialequationsontheinterval(cid:2) 2Œ(cid:2)1;C1(cid:3):
dy
1
Df .(cid:2)/D(cid:4)cos.(cid:4)(cid:2)/; y .(cid:2)1/Dy ; (1)
1 1 10
d(cid:2)
80; (cid:2)16(cid:2) <(cid:2)1=2
dy <
2 Df .(cid:2)/D (cid:4)cos.(cid:4)(cid:2)/; (cid:2)1=26(cid:2) 6C1=2; ; y .(cid:2)1/Dy : (2)
2 2 20
d(cid:2) :0; C1=2<(cid:2) 6C1
Thesolutionstothedifferentialequations(1)and(2)aregiven,respectively,as
y .(cid:2)/Dy Csin.(cid:4)(cid:2)/; (3)
1 10
8y ; (cid:2)16(cid:2) <(cid:2)1=2;
< 20
y .(cid:2)/D y C1Csin.(cid:4)(cid:2)/; (cid:2)1=26(cid:2) 6C1=2; (4)
2 20
:y C2; C1=2<(cid:2) 6C1:
20
Suppose now that it is desired to approximate the solutions to the differential equations (1) and
(2) using the following three different methods that employ the LGR [34] collocation method as
describedinvariousformsinReferences[14–17,19]:(i)apmethodwherethestateisapproximated
using an Nth degree polynomial on Œ(cid:2)1;C1(cid:3) and N is allowed to vary; (ii) an h method using
k k
K equallyspacedmeshintervalswhereK isallowedtovaryandafixedfourth-degreepolynomial
is employed within each mesh interval; and (iii) a ph method where both the number of mesh
intervals, K, and the degree of the approximating polynomial, N , within each mesh interval are
k
allowed to vary. Using any of the aforementioned approximations (p, h, or ph), within any mesh
intervalŒT ;T (cid:3),thefunctionsy .(cid:2)/; .i D1;2/areapproximatedusingthefollowingLagrange
k(cid:2)1 k i
polynomialapproximations:
y.k/.(cid:2)/(cid:3)Y.k/.(cid:2)/DNXC1Y.k/`.k/.(cid:2)/; `.k/.(cid:2)/DNYC1 (cid:2) (cid:2)(cid:2)l.k/ ; (5)
i i ij j j (cid:2).k/(cid:2)(cid:2).k/
jD1 lD1 j l
l¤j
where the support points for `.k/.(cid:2)/; j D 1;:::;N C 1, are the N LGR points [34]
j k k
.(cid:2).k/;:::;(cid:2).k//onŒT ;T (cid:3)suchthat(cid:2).k/ D T and(cid:2).k/ D T isanoncollocatedpointthat
1 N k(cid:2)1 k 1 k(cid:2)1 NC1 k
definestheendofmeshintervalk.WithinanyparticularmeshintervalŒT ;T (cid:3) (cid:4) Œ(cid:2)1;C1(cid:3),the
k(cid:2)1 k
approximationsofy.k/.(cid:2)/; i D1;2,aregivenatthesupportpoints(cid:2).k/ ; j D1;:::;N,as
i jC1
N
y.k/(cid:2)(cid:2).k/ (cid:3)(cid:3)Y.k/(cid:2)(cid:2).k/ (cid:3)DY.k/CXI.k/f (cid:2)(cid:2).k/(cid:3); .i D1;2/; .j D1;:::;N/; (6)
i jC1 i jC1 i1 jl i l
lD1
where Y.k/ is the approximation to y .(cid:2).k// at the start of the mesh interval and I.k/ .j;l D
i1 i 1 jl
1;:::;N / is the N (cid:5)N LGR integration matrix (see Reference [14] for details) defined on the
k k k
meshintervalŒT ;T (cid:3).
k(cid:2)1 k
Suppose now that we define the maximum absolute error in the solution of the differential
equationas
ˇ (cid:2) (cid:3)ˇ
E D max ˇY.k/(cid:2)y (cid:2).k/ ˇ; .i D1;2/:
i ˇ ij i j ˇ
j2Œ1;:::;Nk(cid:2)
k2Œ1;:::;K(cid:2)
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
AphMESHREFINEMENTMETHODFOROPTIMALCONTROL
Figure1(a)and(b)showthebase-10logarithmofE asafunctionofN forthep methodandas
1
a function of K for the h method. First, it is seen that because y .(cid:2)/ is a smooth function, the p
1
methodconvergesexponentiallyasafunctionofN whilethehmethodconvergessignificantlymore
slowly as a function of K. Figure 1(c) and (d) show E . Unlike y , the function y is continuous
2 1 2
but not smooth. As a result, the h method converges faster than the p method because no single
polynomial(regardlessofdegree)onŒ(cid:2)1;C1(cid:3)isabletoapproximatethesolutiontoequation(2)as
accurately as a piecewise polynomial. However, while the h method converges more quickly than
does the p method when approximating the solution of equation (2), it is seen that the h method
doesnotconvergeasquicklyasthepmethoddoeswhenapproximatingthesolutiontoequation(1).
In fact, when approximating the solution of equation (1), it is seen that the h method achieves an
(a) (b)
(c) (d)
(e)
Figure1. Base-10logarithmofabsoluteerrorsinsolutionsofequations(1)and(2)atLagrangepolynomial
supportpointsusingp,h,andphmethods.
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
M.A.PATTERSON,W.W.HAGERANDA.V.RAO
errorof(cid:3) 10(cid:2)7 forK D 24,whereasthep methodconvergesexponentiallyandachievesanerror
of(cid:3)10(cid:2)15forN D20.Asaresult,anhmethoddoesnotprovidethefastestpossibleconvergence
ratewhenapproximatingthesolutiontoadifferentialequationwhosesolutionissmooth.
Given the aforementioned p and h analysis, suppose now that the solution to equation (2) is
approximatedusingtheaforementionedphRadaumethod(i.e.,boththenumberofmeshintervals
and the degree of the approximating polynomial within each mesh interval are allowed to vary).
Assume further that the ph method is constructed such that the time interval Œ(cid:2)1;C1(cid:3) is divided
into three mesh intervals Œ(cid:2)1;(cid:2)1=2(cid:3), Œ(cid:2)1=2;C1=2(cid:3), and ŒC1=2;C1(cid:3), and Lagrange polynomial
approximationsoftheformofequation(5)ofdegreeN ,N ,andN ,respectively,areusedineach
1 2 3
meshinterval.Furthermore,supposeN ,N ,andN areallowedtovary.Becausethesolutiony .(cid:2)/
1 2 3 2
isaconstantinthefirstandthirdmeshintervals,itispossibletosetN DN D2andvaryonlyN .
1 3 2
Figure1(e)showstheerroriny .(cid:2)/,E Dmaxjy (cid:2)Y jusingtheaforementionedthree-interval
2 2ph 2 2
phapproach.Similartotheresultsobtainedusingthep methodwhenapproximatingthesolution
of equation (1), in this case, the error in the solution of equation (2) converges exponentially as a
functionofN .Thus,whileanhmethodmayoutperformapmethodonaproblemwhosesolution
2
is not smooth, it is possible to improve the convergence rate by using a ph-adaptive method. The
foregoing analysis provides a motivation for the development of the ph method described in the
remainderofthispaper.
3. BOLZAOPTIMALCONTROLPROBLEM
Without loss of generality, consider the following general optimal control problem in Bolza form.
Determinethestate,y.t/ 2 Rny,thecontrolu.t/ 2 Rnu,theinitialtime,t ,andtheterminaltime,
0
t ,onthetimeintervalt 2Œt ;t (cid:3)thatminimizethecostfunctional
f 0 f
Z tf
J D(cid:5).y.t /;t ;y.t /;t /C g.y.t/;u.t/;t/dt (7)
0 0 f f
t0
subjecttothedynamicconstraints
dy
Da.y.t/;u.t/;t/; (8)
dt
theinequalitypathconstraints
c 6c.y.t/;u.t/;t/6c ; (9)
min max
andtheboundaryconditions
b 6b.y.t /;t ;y.t /;t /6b : (10)
min 0 0 f f max
Thefunctions(cid:5),g,a,c,andbaredefinedbythefollowingmappings:
(cid:5) W Rny (cid:5)R(cid:5)Rny (cid:5)R ! R;
g W Rny (cid:5)Rnu (cid:5)R ! R;
a W Rny (cid:5)Rnu (cid:5)R ! Rny;
c W Rny (cid:5)Rnu (cid:5)R ! Rnc;
b W Rny (cid:5)R(cid:5)Rny (cid:5)R ! Rnb;
where all vector functions of time are treated as row vectors. In this presentation, it will be useful
to modify the Bolza problem given in equations (7)–(10) as follows. Let (cid:2) 2 Œ(cid:2)1;C1(cid:3) be a new
independentvariablesuchthat
t (cid:2)t t Ct
f 0 f 0
t D (cid:2) C : (11)
2 2
TheBolzaoptimalcontrolproblemofequations(7)–(10)isthendefinedintermsofthevariable(cid:2)
as follows. Determine the state, y.(cid:2)/ 2 Rny, the control, u.(cid:2)/ 2 Rnu, the initial time, t , and the
0
terminaltime,t ,onthetimeinterval(cid:2) 2Œ(cid:2)1;C1(cid:3)thatminimizethecostfunctional
f
t (cid:2)t Z C1
J D(cid:5).y.(cid:2)1/;t ;y.C1/;t /C f 0 g.y.(cid:2)/;u.(cid:2)/;(cid:2)It ;t /d(cid:2) (12)
0 f 0 f
2
(cid:2)1
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
AphMESHREFINEMENTMETHODFOROPTIMALCONTROL
subjecttothedynamicconstraints
dy t (cid:2)t
f 0
D a.y.(cid:2)/;u.(cid:2)/;(cid:2)It ;t /; (13)
0 f
d(cid:2) 2
theinequalitypathconstraints
c 6c.y.(cid:2)/;u.(cid:2)/;(cid:2)It ;t /6c ; (14)
min 0 f max
andtheboundaryconditions
b 6b(cid:4)y.(cid:2)1/;t ;y.C1/;t (cid:5)6b : (15)
min 0 f max
Suppose now that the time interval (cid:2) 2 Œ(cid:2)1;C1(cid:3) is divided into a mesh consisting of K mesh
intervals S D ŒT ;T (cid:3); k D 1;:::;K, where .T ;:::;T / are the mesh points. The mesh
k k(cid:2)1 k 0 K
K K
[ \
intervalsS .k D1;:::;K/havethepropertiesthat S DŒ(cid:2)1;C1(cid:3)and S D;,whilethe
k k k
kD1 kD1
meshpointshavethepropertythat(cid:2)1DT <T <T <(cid:6)(cid:6)(cid:6)<T DC1.Lety.k/.(cid:2)/andu.k/.(cid:2)/
0 1 2 K
bethestateandcontrolinS .TheBolzaoptimalcontrolproblemofequations(12)–(15)canthen
k
berewrittenasfollows.Minimizethecostfunctional
J D(cid:5)(cid:2)y.1/.(cid:2)1/;t ;y.K/.C1/;t (cid:3)C tf (cid:2)t0 XK Z Tk g(cid:2)y.k/.(cid:2)/;u.k/.(cid:2)/;(cid:2)It ;t (cid:3) d(cid:2);
0 f 0 f
2
kD1 Tk(cid:2)1
.k D1;:::;K/;
(16)
subjecttothedynamicconstraints
dy.k/.(cid:2)/ t (cid:2)t (cid:2) (cid:3)
D f 0a y.k/.(cid:2)/;u.k/.(cid:2)/;(cid:2)It ;t ; .k D1;:::;K/; (17)
0 f
d(cid:2) 2
thepathconstraints
(cid:2) (cid:3)
c 6c y.k/.(cid:2)/;u.k/.(cid:2)/;(cid:2)It ;t 6c ; .k D1;:::;K/; (18)
min 0 f max
andtheboundaryconditions
(cid:2) (cid:3)
b 6b y.1/.(cid:2)1/;t ;y.K/.C1/;t 6b : (19)
min 0 f max
Because the state must be continuous at each interior mesh point, it is required that the condition
y.T(cid:2)/Dy.TC/; .k D1;:::;K(cid:2)1/besatisfiedattheinteriormeshpoints.T ;:::;T /.
k k 1 K(cid:2)1
4. LEGENDRE–GAUSS–RADAUCOLLOCATIONMETHOD
Thephformofthecontinuous-timeBolzaoptimalcontrolprobleminSection3isdiscretizedusing
collocationatLGRpoints[14–17,19].IntheLGRcollocationmethod,thestateofthecontinuous-
timeBolzaoptimalcontrolproblemisapproximatedinS ; k 2Œ1;:::;K(cid:3),as
k
y.k/.(cid:2)/(cid:3)Y.k/.(cid:2)/DNXkC1Y.k/`.k/.(cid:2)/; `.k/.(cid:2)/DNYkC1 (cid:2) (cid:2)(cid:2)l.k/ ; (20)
j j j (cid:2).k/(cid:2)(cid:2).k/
jD1 lD1 j l
l¤j
where (cid:2) 2 Œ(cid:2)1;C1(cid:3); `.k/.(cid:2)/; j D 1;:::;N C 1, is a basis of Lagrange polynomials,
j k
(cid:2) (cid:3)
(cid:2).k/;:::;(cid:2).k/ are the LGR [34] collocation points in S D ŒT ;T /, and (cid:2).k/ D T is a
1 Nk k k(cid:2)1 k NkC1 k
noncollocatedpoint.DifferentiatingY.k/.(cid:2)/inequation(20)withrespectto(cid:2),weobtain
dY.k/.(cid:2)/ DNXkC1Y.k/d`.jk/.(cid:2)/: (21)
d(cid:2) j d(cid:2)
jD1
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
M.A.PATTERSON,W.W.HAGERANDA.V.RAO
The cost functional of equation (16) is then approximated using a multiple-interval LGR
quadratureas
J (cid:3)(cid:5)(cid:2)Y.1/;t ;Y.K/ ;t (cid:3)CXK XNk tf (cid:2)t0w.k/g(cid:2)Y.k/;U.k/;(cid:2).k/It ;t (cid:3); (22)
1 0 NKC1 f 2 j j j j 0 f
kD1jD1
where w.k/ .j D 1;:::;N / are the LGR quadrature weights [34] in S D ŒT ;T (cid:3), k 2
j k k k(cid:2)1 k
Œ1;:::;K(cid:3), U.k/; i D 1;:::;N , are the approximations of the control at the N LGR points in
i k k
meshintervalk 2Œ1;:::;K(cid:3); Y.1/istheapproximationofy.T /,andY.K/ istheapproximation
1 0 NKC1
ofy.T /(wherewerecallthatT D(cid:2)1andT DC1).Collocatingthedynamicsofequation(17)
K 0 K
attheN LGRpointsusingequation(21),wehave
k
NXkC1D.k/Y.k/(cid:2) tf (cid:2)t0a(cid:2)Y.k/;U.k/;(cid:2).k/It ;t (cid:3)D0; .i D1;:::;N /;
ij j 2 i i i 0 f k
jD1
where
d`.k/.(cid:2).k//
D.k/ D j i ; .i D1;:::;N ; j D1;:::;N C1/;
ij d(cid:2) k k
aretheelementsoftheN (cid:5).N C1/LGRdifferentiationmatrix[14]D.k/associatedwithS ; k 2
k k k
Œ1;:::;K(cid:3). While the dynamics can be collocated in differential form, in this paper, we choose
to collocate the dynamics using the equivalent implicit integral form (see References [14–16] for
details).TheimplicitintegralformoftheLGRcollocationmethodisgivenas
Y.k/ (cid:2)Y.k/(cid:2) tf (cid:2)t0 XNk I.k/a(cid:2)Y.k/;U.k/;(cid:2).k/It ;t (cid:3)D0; .i D1;:::;N /; (23)
iC1 1 2 ij i i i 0 f k
jD1
whereI.k/; .i D1;:::;N ;j D1;:::;N ;k D1;:::;K/istheN (cid:5)N LGRintegrationmatrix
ij k k k k
inmeshintervalk 2Œ1;:::;K(cid:3);itisobtainedbyinvertingasubmatrixofthedifferentiationmatrix
formedbycolumns2throughN C1:
k
h i(cid:2)1
I.k/ D D.k/(cid:6)(cid:6)(cid:6)D.k/ :
2 NkC1
ItisnotedforcompletenessthatI.k/D.k/ D(cid:2)1(seeReferences[14–16]),where1isacolumnvector
1
oflengthN ofallones[14–16].Next,thepathconstraintsofequation(18)inS ; k 2Œ1;:::;K(cid:3),
k k
areenforcedattheN LGRpointsas
k
(cid:2) (cid:3)
c 6c Y.k/;U.k/;(cid:2).k/It ;t 6c ; .i D1;:::;N /: (24)
min i i i 0 f max k
Theboundaryconditionsofequation(19)areapproximatedas
(cid:2) (cid:3)
b 6b Y.1/;t ;Y.K/ ;t 6b : (25)
min 1 0 NKC1 f max
Itisnotedthatcontinuityinthestateattheinteriormeshpointsk 2 Œ1;:::;K (cid:2)1(cid:3)isenforcedvia
thecondition
Y.k/ DY.kC1/; .k D1;:::;K(cid:2)1/; (26)
NkC1 1
wherethesamevariableisusedforbothY.k/ andY.kC1/.Hence,theconstraintofequation(26)
NkC1 1
iseliminatedfromtheproblembecauseitistakenintoaccountexplicitly.TheNLPthatarisesfrom
the LGR collocation method is then to minimize the cost function of equation (22) subject to the
algebraicconstraintsofequations(23)–(25).
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
AphMESHREFINEMENTMETHODFOROPTIMALCONTROL
5. ph-ADAPTIVEMESHREFINEMENTMETHOD
We now develop a ph-adaptive mesh refinement method using the LGR collocation method
describedinSection4.Wecallourmethodaphmethodbecausewefirsttrytoadjustthepolyno-
mialdegreetoachieveconvergence,andifthisfails,weadjustthemeshspacing.Theph-adaptive
meshrefinementmethoddevelopedinthispaperisdividedintotwoparts.InSection5.1,themethod
forestimatingtheerrorinthecurrentsolutionisderived,andinSection5.4,thep-then-hstrategy
isdevelopedforrefiningthemesh.
5.1. Errorestimateineachmeshinterval
In this section, an estimate of the relative error in the solution within a mesh interval is derived.
BecausethestateistheonlyquantityintheLGRcollocationmethodforwhichauniquelydefined
function approximation is available, we develop an error estimate for the state. The error estimate
isobtainedbycomparingtwoapproximationstothestate,onewithhigheraccuracy.Thekeyidea
is that for a problem whose solution is smooth, an increase in the number of LGR points should
yieldastatethatmoreaccuratelysatisfiesthedynamics.Hence,thedifferencebetweenthesolution
associatedwiththeoriginalsetofLGRpointsandtheapproximationassociatedwiththeincreased
numberofLGRpointsshouldyieldanestimatefortheerrorinthestate.
Assume that the NLP of equations (22)–(25) corresponding to the discretized control problem
has been solved on a mesh S D ŒT ;T (cid:3); k D 1;:::;K, with N LGR points in mesh inter-
k k(cid:2)1 k k
val S . Suppose that we want to estimate the error in the state at a set of M D N C 1 LGR
k(cid:2) (cid:3) k k
points (cid:2)O.k/;:::;(cid:2)O.k/ , where (cid:2)O.k/ D (cid:2).k/ D T , and that (cid:2)O.k/ D T . Suppose further
1 Mk 1 1 k(cid:2)1 MkC1 k
(cid:2) (cid:3)
thatthevalues ofthestateapproximationgiven inequation(20)atthepoints (cid:2)O.k/;:::;(cid:2)O.k/ are
1 Mk
(cid:2) (cid:3)
denoted Y.(cid:2)O.k//;:::;Y.(cid:2)O.k// . Next, let the control be approximated in S using the Lagrange
1 Mk k
interpolatingpolynomial
U.k/.(cid:2)/D XNk U.k/`O.k/.(cid:2)/; `O.k/.(cid:2)/D YNk (cid:2) (cid:2)(cid:2)l.k/ ; (27)
j j j (cid:2).k/(cid:2)(cid:2).k/
jD1 lD1 j l
l¤j
andletthecontrolapproximationat(cid:2)O.k/ bedenotedU.(cid:2)O.k//,1 6 i 6 M .Weusethevalueofthe
i i k
right-handsideofthedynamicsat.Y.(cid:2)O.k//;U.(cid:2)O.k//;(cid:2)O.k//toconstructanimprovedapproximation
i i i
ofthestate.LetYO.k/ beapolynomialofdegreeatmostM thatisdefinedontheintervalS .Ifthe
k k
derivativeofyO.k/ matchesthedynamicsateachoftheRadauquadraturepoints(cid:2)O.k/; 1 6 i 6 M ,
i k
thenwehave
YO.k/(cid:2)(cid:2)O.k/(cid:3)DY.k/.(cid:2) /Ctf (cid:2)t0 XMk IO.k/a(cid:2)Y.k/(cid:2)(cid:2)O.k/(cid:3);U.k/(cid:2)(cid:2)O.k/(cid:3);(cid:2)O.k/(cid:3); jD2;:::;M C1;
j k(cid:2)1 2 jl l l l k
lD1
(28)
whereIO.k/; j;l D 1;:::;M ,istheM (cid:5)M LGRintegrationmatrixcorrespondingtotheLGR
jl k k k
(cid:2) (cid:3)
pointsdefinedby (cid:2)O.k/;:::;(cid:2)O.k/ .UsingthevaluesY.(cid:2)O.k//andYO.(cid:2)O.k//,l D 1;:::;M C1,the
1 Mk l l k
absoluteandrelativeerrorsintheithcomponentofthestateat.(cid:2)O.k/;:::;(cid:2)O.k/ /arethendefined,
1 MkC1
respectively,as
(cid:2) (cid:3) ˇ (cid:2) (cid:3) (cid:2) (cid:3)ˇ
E.k/ (cid:2)O.k/ DˇYO.k/ (cid:2)O.k/ (cid:2)Y.k/ (cid:2)O.k/ ˇ;
i l ˇ i l i l ˇ
(cid:2) (cid:3) (cid:6) (cid:7)
ei.k/(cid:2)(cid:2)Ol.k/(cid:3) D 1C mEai.xk/ (cid:2)Ol.ˇˇkY/.k/(cid:2)(cid:2)O.k/(cid:3)ˇˇ; l Di D1;:1:;::;:M:;knyC;1; : (29)
ˇ i j ˇ
j2Œ1;:::;MkC1(cid:2)
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca
M.A.PATTERSON,W.W.HAGERANDA.V.RAO
ThemaximumrelativeerrorinS isthendefinedas
k
e.k/ D max e.k/.(cid:2)O.k//: (30)
max i l
i2Œ1;:::;ny(cid:2)
l2Œ1;:::;MkC1(cid:2)
5.2. Rationaleforerrorestimate
TheerrorestimatederivedinSection5.1issimilartotheerrorestimateobtainedusingthemodified
EulerRunge–KuttaschemetonumericallysolveadifferentialequationyP.t/ D f.y.t//.Thefirst-
orderEulermethodisgivenas
y Dy Chf.y /; (31)
jC1 j j
where h is the step size and y is the approximation to y.t/ at t D t D jh. In the second-order
j j
modified Euler Runge–Kutta method, the first stage generates the following approximation yN to
y.t /:
jC1=2
yN Dy C 1hf.y /:
j 2 j
The second stage then uses the dynamics evaluated at yN to obtain an improved estimate yO of
jC1
y.t /:
jC1
yO Dy Chf.yN/: (32)
jC1 j
The original Euler scheme starts at y and generates y . The first-stage variable yN is the
j jC1
interpolant of the line (first-degree polynomial) connecting .t ;y / and .t ;y / evaluated at
j j jC1 jC1
thenewpointt .Thesecondstagegiveninequation(32)usesthedynamicsattheinterpolant
jC1=2
yN toobtainanimprovedapproximationtoy.t /.BecauseyO isasecond-orderapproximation
jC1 jC1
toy.t /andy isafirst-orderapproximationtoy ,theabsolutedifferencejyO (cid:2)y j
jC1 jC1 jC1 jC1 jC1
isanestimateoftheerroriny inamannersimilartotheabsoluteerrorestimateE.k/.(cid:2)O.k//.l D
jC1 i l
1;:::;M C1/derivedinequation(29).
k
The effectiveness of the derived error estimate derived in Section 5.1 can be seen by revisiting
themotivatingexamplesofSection2.Figure2(a)and(b)showthepandherrorestimates,respec-
tively, E and E , in the solution to equation (1); (c) and (d) show the p and h error estimates,
1p 1h
respectively, E and E , in the solution to equation (2); and (e) shows the ph error estimates,
2p 2h
E , in the solution to equation (2). It is seen that the error estimates are nearly identical to the
2ph
actualerror.Therelativeerrorestimategiveninequation(30)isusedinthenextsectionasthebasis
formodifyinganexistingmesh.
5.3. Estimationofrequiredpolynomialdegreewithinameshinterval
SupposeagaintheLGRcollocationNLPofequations(22)–(25)hasbeensolvedonameshS ; k D
k
1;:::;K.Supposefurtherthatitisdesiredtomeetarelativeerroraccuracytolerance(cid:6)ineachmesh
intervalS ; k D 1;:::;K.Ifthetolerance(cid:6) isnotmetinatleastonemeshinterval,thenthenext
k
stepistorefinethecurrentmesh,eitherbydividingthemeshintervalorincreasingthedegreeofthe
approximatingpolynomialwithinthemeshinterval.
ConsiderameshintervalS ; q 2Œ1;:::;K(cid:3),whereN LGRpointswereusedtosolvetheNLP
q q
ofequations(22)–(25),andagain,let(cid:6)bethedesiredrelativeerroraccuracytolerance.Supposefur-
therthattheestimatedmaximumrelativeerror,e.q/,hasbeencomputedasdescribedinSection5.1
max
andthate.q/ >(cid:6)(i.e.,theaccuracytolerance(cid:6)isnotsatisfiedintheexistingmeshinterval).Finally,
max
letN andN beuser-specifiedminimumandmaximumboundsonthenumberofLGRpoints
min max
withinanymeshinterval.Accordingtotheconvergencetheorysummarizedin[35,36],theerrorin
aglobalcollocationschemebehaveslikeO.N2:5(cid:2)k/,whereN isthenumberofcollocationpoints
withinameshintervalandk isthenumberofcontinuousderivativesinthesolution[21,22].Ifthe
solutionissmooth,thenwecouldtakek DN.Hence,ifN wasreplacedbyN CP,thentheerror
bounddecreasesbyatleastthefactorN(cid:2)P.
Copyright©2014JohnWiley&Sons,Ltd. Optim.ControlAppl.Meth.(2014)
DOI:10.1002/oca