Solving optimal control problems with MATLAB 使用 MATLAB 解決最優(yōu)控制問題

PAGE

10

SolvingoptimalcontrolproblemswithMATLAB—Indirectmethods

XuezhongWang?

1 Introduction

Thetheoryofoptimalcontrolhasbeenwelldevelopedforoverfortyyears.Withtheadvancesofcomputertechnique,optimalcontrolisnowwidelyusedinmulti-disciplinaryapplicationssuchasbiologicalsystems,communi-cationnetworksandsocio-economicsystemsetc.Asaresult,moreandmorepeoplewillbenefitgreatlybylearningtosolvetheoptimalcontrolprob-lemsnumerically.Realizingsuchgrowingneeds,booksonoptimalcontrolputmoreweightonnumericalmethods.Inretrospect,[1]wasthefirstandthe“classic”bookforstudyingthetheoryaswellasmanyinterestingcases(time-optimal,fuel-optimalandlinearquadraticregulator(LQR)problems).Necessaryconditionsforvarioussystemswerederivedandexplicitsolutionsweregivenwhenpossible.Later,[2]provedtobeaconciseyetexcellentbookwithmoreengineeringexamples.Oneofthedistinguishfeaturesofthisbookisthatitintroducedseveraliterativealgorithmsforsolvingproblemsnumer-ically.Morerecently,[3]usesMATLABtosolveproblemswhichiseasierandmoreprecise.However,thenumericalmethodscoveredinthesebooksareinsufficientforthewiderangeofproblemsemergingfromvariousfields.Especially,forthoseproblemswithfreefinaltimeandnonlineardynamics.

ThistutorialshowscommonroutinesinMATLABtosolvebothfixedandfreefinaltimeproblems.Specifically,thefollowingsubjectsarediscussedwithexamples:

HowtouseSymbolicMathToolboxtoderivenecessaryconditionsandsolveforexplicitsolutions?

Howtosolveafixed-final-timeoptimalcontrolproblemwithsteepestdescentmethod?

?ISE.Dept.,NCSU,Raleigh,NC27695(xwang10@)

HowtoreformulatetheoriginalproblemasaBoundaryValueProblem(BVP)andsolveitwithbvp4c?

Howtoget‘goodenough’solutionswithbvp4candunderwhatcondi-tionswillbvp4cfailtofindasolution?

Itshouldbenotedthatalltheroutines(exceptthesteepestdescentmethod)discussedinthistutorialbelongtothe“indirectmethods”category.Thismeansthatconstraintsonthecontrolsandstatesarenotconsidered.1Inotherwords,thecontrolcanbesolvedintermsofstatesandcostatesandtheproblemisequivalentlytoaBVP.Thereferenceofalltheexamplesusedinthistutorialarestatedsuchthattheresultscanbecomparedandverified.

2 Optimalcontrolproblemswithfixed-final-time

Inmostbooks[1][2],itisfree-final-timeproblemthatbeingtackledfirsttoderivethenecessaryconditionsforoptimalcontrol.Fixed-final-timeprob-lemsweretreatedasanequivalentvariationwithonemorestatefortime.However,fornumericalmethods,fixed-final-timeproblemsarethegeneralformandwesolvethefree-final-timeproblembyconvertingitintoafixed-final-timeproblem.Thereasonissimple:whendealingwithoptimalcontrolproblems,itisinevitabletodonumericalintegration(eitherbyindirectordirectmethods).Therefore,atimeintervalmustbespecifiedforthesemeth-ods.

Thefirstproblemisfrom[2],Example5.1-1frompage198to202.TheroutinedeployedinthisexampleshowshowtoderivenecessaryconditionsandsolveforthesolutionswithMATLABSymbolicMathToolbox.

Example1Thesystem

x˙1(t)=x2(t) (1)

x˙2(t)=?x2(t)+u(t) (2)

/

Considertheperformancemeasure:

tf1

J(u)=

u2(t)dt

0 2

1Thelastexampleimplementsasimpleconstraintonthecontroltoshowlimitationsofbvp4csolver.

withboundaryconditions:

［ ］

x(0)=0 x(2)= 52T

Considertheperformancemeasure:

J(u)=

1

(x1(2)?5)2+

1

(x2(2)?2)2+

/21

u2(t)dt

2 2 02

withboundaryconditions:

x(0)=0 x(2)isfree

Considertheperformancemeasureasin(a)withfollowingboundarycon-ditions:

x(0)=0 x1(2)+5x2(2)=15

ThefirststepistoformtheHamiltonianandapplynecessaryconditionsforoptimality.WithMATLAB,thiscanbedoneasfollows:

%Stateequationssymsx1x2p1p2u;Dx1=x2;

Dx2=-x2+u;

%Costfunctioninsidetheintegralsymsg;

g=0.5*u^2;

%Hamiltoniansymsp1p2H;

H=g+p1*Dx1+p2*Dx2;

%CostateequationsDp1=-diff(H,x1);

Dp2=-diff(H,x2);

%solveforcontroludu=diff(H,u);

sol_u=solve(du,’u’);

TheMATLABcommandsweusedherearediffandsolve.diffdifferen-tiatesasymbolicexpressionandsolvegivessymbolicsolutiontoalgebraic

equations.FormoredetailsaboutSymbolicMathToolbox,pleasereferto[5].Applyingthenecessaryconditionsforoptimalitywegettwoequations:2

i

p˙?

?H

?H

i

=??x?

(3)

?u?=0 (4)

Thefirstequationgivescostateequations.Fromthesecondequation,wesolveforcontroluintermsofstatesandcostates.Thesecondstepistosubstituteufrom(4)backtothestateandcostateequationstogetasetof2nfirst-orderordinarydifferentialequations(ODE’s).Asolution(with2narbitrarycoefficients)canbeobtainedbyusingthedsolvecommandwithoutanyboundaryconditions.Thesymbolicsolutionlooksdifferentfromtheonein[2].Bysimplifyingtheexpressionandrearrangethearbitrarycoefficients,itisnotdifficulttoseethatthetwoarethesame.

%SubstituteutostateequationsDx2=subs(Dx2,u,sol_u);

%convertsymbolicobjectstostringsforusing’dsolve’eq1=strcat(’Dx1=’,char(Dx1));

eq2=strcat(’Dx2=’,char(Dx2));eq3=strcat(’Dp1=’,char(Dp1));eq4=strcat(’Dp2=’,char(Dp2));

sol_h=dsolve(eq1,eq2,eq3,eq4);

Asstatedin[2],thedifferencesofthethreecasesinthisproblemaremerelytheboundaryconditions.For(a),thearbitrarycoefficientscanbedeterminedbysupplyingthe2nboundaryconditionstodsovle:

% casea:(a)x1(0)=x2(0)=0;x1(2)=5;x2(2)=2;conA1=’x1(0)=0’;

conA2=’x2(0)=0’;

conA3=’x1(2)=5’;

conA4=’x2(2)=2’;

sol_a=dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);

AgainthesolutionsgivenbyMATLABand[2]lookdifferentfromeachother.YetFigure1showsthatthetwoareinfactequivalent.Forallthefiguresinthisproblem,*representstatetrajectoryfrom[2]whilesymbolicsolutionfromMATLABisplottedwithacontinuousline.

2Weuse*toindicatetheoptimalstatetrajectoryorcontrol.

Figure1:TrajectoriesforExample1(a)

Forcase(b),wesubstitutet0=0,tf=2inthethegeneralsolutiontogetasetof4algebraicequationswith4arbitrarycoefficients.Thenthefourboundaryconditionsaresuppliedtodeterminethesecoefficients.sol_bisastructureconsistsoffourcoefficientsreturnedbysolve.WegetthefinalsymbolicsolutiontotheODE’sbysubstitutingthesecoefficientsintothegeneralsolution.Figure2showstheresultsofbothsolutionsfromMATLABand[2].

Itshouldbenotedthatwecannotgetthesolutionbydirectlysupplyingtheboundaryconditionswithdsolveaswedidincase(a).Becausethe

lattertwoboundaryconditions:p?(2)=x?(t)?5,p?(2)=x?(t)?2replaces

1 1 2 2

costatesp1,p2withstatesx1,x2intheODE’swhichresultedinmoreODE’s

thanvariables.

% caseb:(a)x1(0)=x2(0)=0;p1(2)=x1(2)-5;p2(2)=x2(2)-2;

eq1b=char(subs(sol_h.x1,’t’,0));

eq2b=char(subs(sol_h.x2,’t’,0));

eq3b=strcat(char(subs(sol_h.p1,’t’,2)),...

’=’,char(subs(sol_h.x1,’t’,2)),’-5’);

eq4b=strcat(char(subs(sol_h.p2,’t’,2)),...

’=’,char(subs(sol_h.x2,’t’,2)),’-2’);sol_b=solve(eq1b,eq2b,eq3b,eq4b);

Figure2:TrajectoriesforExample1(b)

%SubstitutethecoefficientsC1=double(sol_b.C1);

C2=double(sol_b.C2);C3=double(sol_b.C3);C4=double(sol_b.C4);

sol_b2=struct(’x1’,{subs(sol_h.x1)},’x2’,{subs(sol_h.x2)},...

’p1’,{subs(sol_h.p1)},’p2’,{subs(sol_h.p2)});

Case(c)isalmostthesameascase(b)withslightlydifferentbound-aryconditions.From[2],theboundaryconditionsare:x?(2)+5x?(2)=

1 2

15,p?(2)=5p?(2).Figure3showstheresultsofbothsolutionsfromMAT-

2 1

LABand[2].

%casec:x1(0)=x2(0)=0;x1(2)+5*x2(2)=15;p2(2)=5*p1(2);

eq1c=char(subs(sol_h.x1,’t’,0));

eq2c=char(subs(sol_h.x2,’t’,0));

eq3c=strcat(char(subs(sol_h.p2,’t’,2)),...

’-(’,char(subs(sol_h.p1,’t’,2)),’)*5’);eq4c=strcat(char(subs(sol_h.x1,’t’,2)),...

’+(’,char(subs(sol_h.x2,’t’,2)),’)*5-15’);sol_c=solve(eq1c,eq2c,eq3c,eq4c);

%SubstitutethecoefficientsC1=double(sol_c.C1);

Figure3:TrajectoriesforExample1(c)

C2=double(sol_c.C2);C3=double(sol_c.C3);C4=double(sol_c.C4);

sol_c2=struct(’x1’,{subs(sol_h.x1)},’x2’,{subs(sol_h.x2)},...

’p1’,{subs(sol_h.p1)},’p2’,{subs(sol_h.p2)});

Ifthestateequationsarecomplicated,itisusuallyimpossibletoderiveexplicitsolutions.Wedependonnumericalmethodstofindthesolutions.Inthenextexample,wewillillustratetwonumericalroutines:steepestdescentmethodandconverttoaBVP.Theproblemcanbefoundin[2]page338to339,Example6.2-2.

( )1 1 2 1

Example2Thestateequationsforacontinuousstirred-tankchemicalre-actoraregivenas:

x1(t)+2

x˙(t)=?2[x(t)+0.25]+[x(t)+0.5]exp 25x1(t) ?[x(t)+0.25]u(t)

x1(t)+2

x˙2(t)=0.5?x2(t)?[x2(t)+0.5]exp(25x1(t))

(5)

Theflowofacoolantthroughacoilinsertedinthereactoristocontrolthefirst-order,irreversibleexothermicreactiontakingplaceinthereactor.x1(t)=T(t)isthedeviationfromthesteady-statetemperatureandx2(t)=

C(t)isthedeviationfromthesteady-stateconcentration.u(t)isthenormal-izedcontrolvariable,representingtheeffectofcoolantflowonthechemicalreaction.

Theperformancemeasuretobeminimizedis:

J=

[x2(t)+x2(t)+Ru2(t)]dt,

/0.78

1 2

0

withtheboundaryconditions

f

x(0)=［0.050］T,x(t)isfree

First,wewillimplementthesteepestdescentmethodbasedontheschemeoutlinedin[2].Thealgorithmconsistsof4steps:

Subdividetheinterval[t0,tf]intoNequalsubintervalsandassumeapiecewise-constantcontrolu(0)(t)=u(0)(tk),t∈[tk,tk+1]k=0,1,···,N?1

Applyingtheassumedcontrolu(i)tointegratethestateequationsfromt0totfwithinitialconditionsx(t0)=x0andstorethestatetrajectoryx(i).

Applyingu(i)andx(i)tointegratecostateequationsbackward,i.e.,from[tf,t0].The“initialvalue”p(i)(tf)canbeobtainedby:

p(i)(t

)=?h(x(i)(t

)).

f

?x

f

Evaluate?H(i)(t)/?u,t∈[t0,tf]andstorethisvector.

1?H(i)1

If

1

1?u

≤γ (6)

(i)2

1

1

1?H 1

?u

/tfr1?H(i)11Tr1?H(i)11

t0

?u

?u

≡

1

1

1

1

dt (7)

thenstoptheiterativeprocedure.Hereγisapreselectedsmallpositiveconstantusedasatolerance.

If(6)isnotsatisfied,adjustthepiecewise-constantcontrolfunctionby:

u(i+1)

(tk)=u

(i)

(tk)?τ

?H(i)

?u(tk), k=0,1,···,N?1

Replaceu(i)byu(i+1)andreturntostep2.Here,τisthestepsize.

ThemainloopinMATLABisasfollows:

fori=1:max_iteration

%1)startwithassumedcontroluandmoveforward[Tx,X]=ode45(@(t,x)stateEq(t,x,u,Tu),[t0tf],...

initx,options);

%2)Movebackwardtogetthetrajectoryofcostatesx1=X(:,1);x2=X(:,2);

[Tp,P]=ode45(@(t,p)costateEq(t,p,u,Tu,x1,x2,Tx),...[tft0],initp,options);

p1=P(:,1);

%Important:costateisstoredinreverseorder.Thedimensionof

%costatesmayalsodifferentfromdimensionofstates

%Useinterploatetomakesurexandpisalignedalongthetimeaxisp1=interp1(Tp,p1,Tx);

%CalculatedeltaHwithx1(t),x2(t),p1(t),p2(t)dH=pH(x1,p1,Tx,u,Tu);

H_Norm=dH’*dH;

%Calculatethecostfunction

J(i,1)=tf*(((x1’)*x1+(x2’)*x2)/length(Tx)+...

0.1*(u*u’)/length(Tu));

end

%ifdH/du<epslon,exitifH_Norm<eps

%DisplayfinalcostJ(i,1)

break;else

%adjustcontrolfornextiterationu_old=u;

u=AdjControl(dH,Tx,u_old,Tu,step);end;

Becausethestepsizeofode45isnotpredetermined,interpolationisusedtomakesurex(t),p(t)andu(t)arealignedalongthetimeaxis.

%Stateequations

functiondx=stateEq(t,x,u,Tu)dx=zeros(2,1);

u=interp1(Tu,u,t);%Interploatethecontrolattimet

dx(1)=-2*(x(1)+0.25)+(x(2)+0.5)*exp(25*x(1)/(x(1)+2))...

-(x(1)+0.25).*u;

dx(2)=0.5-x(2)-(x(2)+0.5)*exp(25*x(1)/(x(1)+2));

%Costateequations

functiondp=costateEq(t,p,u,Tu,x1,x2,xt)dp=zeros(2,1);

x1=interp1(xt,x1,t); %Interploatethestatevarialbesx2=interp1(xt,x2,t);

u=interp1(Tu,u,t); %Interploatethecontrol

dp(1)=p(1).*(u+exp((25*x1)/(x1+2)).*((25*x1)/(x1+2)^2-...

25/(x1+2))*(x2+1/2)+2)-...

2*x1-p(2).*exp((25*x1)/(x1+2))*((25*x1)/(x1+2)^2-...

25/(x1+2))*(x2+1/2);

dp(2)=p(2).*(exp((25*x1)/(x1+2))+1)-...

p(1).*exp((25*x1)/(x1+2))-2*x2;

Instep4,wecalculate?H/?u|t=2Ru(t)?p1(t)[x1(t)+0.25]oneachtimesubinterval(functionpH)andcomparel?H/?ul2withthepreselectedγ.Ifl?H/?ul2>γ,adjustu(functionAdjControl).

%PartialderivativeofHwithrespecttoufunctiondH=pH(x1,p1,tx,u,Tu)

%interploatethecontrolu=interp1(Tu,u,tx);

R=0.1;

dH=2*R*u-p1.*(x1+0.25);

%Adjustthecontrol

functionu_new=AdjControl(pH,tx,u,tu,step)

%interploatedH/du

pH=interp1(tx,pH,tu);u_new=u-step*pH;

Thestepsizeτissetasaconstantinthisexamplebysomeposthocinvestigation.Abetterstrategyistoselectτwithalinesearchmethodwhichwillmaximizethereductionofperformancemeasurewithgiven?H(i)/?uineachiteration.Figure4showstheoptimalstatetrajectoryandcontroloverthetime.ThevalueofperformancemeasureasafunctionofiterationnumberisshowninFigure5.In[2],twomorenumericalalgorithmswereintroduced:variationofextremalsandquasilinearization.ThesetwomethodsbasicallyreformulateandsolvetheoriginalproblemasaBoundaryValueProblem(BVP).InMATLAB,aBVPistypicallysolvedwithbvp4c.[4]isanexcellentreferenceonusingbvp4c.Forfix-final-timeproblems,ucanalwaysbesolved

Figure4:Example2Steepestdescentmethod

Figure5:Performancemeasurereduction

withrespecttoxandpbyapplyingthenecessaryconditions(3)and(4).Andwewillhave2nODE’sand2nboundaryconditions.

Figure6:Solutionfrombfv4cforProblem2

%Initialguessforthesolution

solinit=bvpinit(linspace(0,0.78,50),...

[000.50.5]);

options=bvpset(’Stats’,’on’,’RelTol’,1e-1);globalR;

R=0.1;

sol=bvp4c(@BVP_ode,@BVP_bc,solinit,options);t=sol.x;

y=sol.y;

%Calculateu(t)fromx1,x2,p1,p2

ut=(y(3,:).*(y(1,:)+1/4))/(2*0.1);

n=length(t);

%Calculatethecost

J=0.78*(y(1,:)*y(1,:)’+y(2,:)*y(2,:)’+...

ut*ut’*0.1)/n;

TheODE’sandboundaryconditionsaretwomajorconsiderationsinusingbvp4c.AruleofthumbisthatthenumberofODE’smustequaltothenumberofboundaryconditionssuchthattheproblemissolvable.OncetheoptimalcontrolproblemisconvertedtoaBVP,itisverysimpletosolve.

%

%ODE’sforstatesandcostates

functiondydt=BVP_ode(t,y)globalR;

t1=y(1)+.25;

t2=y(2)+.5;

t3=exp(25*y(1)/(y(2)+2));t4=50/(y(1)+2)^2;

u=y(3)*t1/(2*R);

dydt=[-2*t1+t2*t3-t2*u

0.5-y(2)-t2*t3

-2*y(1)+2*y(3)-y(3)*t2*t4*t3+y(3)*u+y(4)*t2*t4*t3

-2*y(2)-y(3)*t3+y(4)*(1+t3)];

%

%Theboundaryconditions:

%x1(0)=0.05,x2(0)=0,tf=0.78,p1(tf)=0,p2(tf)=0;

functionres=BVP_bc(ya,yb)res=[ya(1)-0.05

ya(2)-0

yb(3)-0yb(4)-0];

Inthisexample,bvp4cworksperfectly.Itisfasterandgivesbetterre-sults,i.e.asmallerperformancemeasureJcomparingtothesteepestdescentmethod(seeFigure6).Inthefollowingsection,wewillsolelyusebvp4cwhennumericalsolutionsareneeded.

3 Optimalcontrolproblemswithfree-final-time

Nowwearepreparedtodealwithfree-final-timeproblems.WewillusebothSymbolicMathToolboxandbvp4cinthenextexample,whichcanbefindfrom[3]onpage77,Example2.14.

Example3Givenadoubleintegralsystemas:

x˙1(t)=x2(t) (8)

x˙2(t)=u(t) (9)

Theperformancemeasureis:

J=

1/tf

u2(t)dt

20

findtheoptimalcontrolgiventheboundaryconditionsas:

1f 2f

x(0)=［12］T,x(t)=3,x(t)isfree

TousetheSymbolicMathToolbox,theroutineisverysimilartoProblem

1.WefirstsupplytheODE’sandboundaryconditionsonstatesandcostatestodsolve.Theonlydifferenceisthatthefinaltimetfitselfisnowavariable.Asaresult,thesolutionisafunctionoftf.Next,weintroducefourmorevariables,namelyx1(tf),x2(tf),p1(tf),p2(tf)intothesolutionobtainedabove.Withoneadditionalboundaryconditionfrom

H(x?(t

),u?(t),p?(t),t)+?h(x?(t),t

)=0

f f f f ?t f f

2

Forthisproblem,h≡0andwehavep1(tf)x2(tf)?0.5p2(tf)=0.Nowwehave5algebraicequationswith5unknowns.Andsolvecomesinhandytosolvethisproblem.Figure7showstheresultsfromMATLABandtheanalyticalsolutionin[3].AlthoughSymbolicMathToolboxworksfineinthisexample,itshouldbepointedoutthatinmostproblems,itisimpossibletogetexplicitsolutions.Forexample,thereisnoexplicitsolutionsforExample1eventhoughthestateequationsaresimilartothoseofExample3.

sol=dsolve(’Dx1=x2,Dx2=-p2,Dp1=0,Dp2=-p1’,...’x1(0)=1,x2(0)=2,x1(tf)=3,p2(tf)=0’);

eq1=subs(sol.x1)-’x1tf’;eq2=subs(sol.x2)-’x2tf’;eq3=subs(sol.p1)-’p1tf’;eq4=subs(sol.p2)-’p2tf’;

eq5=sym(’p1tf*x2tf-0.5*p2tf^2’);

%%

sol_2=solve(eq1,eq2,eq3,eq4,eq5);tf=sol_2.tf;

x1tf=sol_2.x1tf;x2tf=sol_2.x2tf;

x1=subs(sol.x1);x2=subs(sol.x2);p1=subs(sol.p1);p2=subs(sol.p2);

Becauseofthelimitationsofsymbolicmethod,numericalmethodsaremoreusefulindealingwithmoregeneralproblems.However,whenwetrytouseanumericalmethodsuchasbvp4c,weimmediatelyencounteredwith

Figure7:Example3Symbolicmethod

aproblem:thetimeintervalisnotknown.Onecommontreatment[4][7]forsuchasituationistochangetheindependentvariablettoτ=t/T,the

augmentedstateandcostateequationswillthenbecomex?˙

=Tf(x,q,τ).3

Nowtheproblemisposedonfixedinterval[0,1].Thiscanbeimplementedinbvp4cbytreatingTasanauxiliaryvariable.Thefollowingcodesnippetshowsthedetails.

solinit=bvpinit(linspace(0,1),[2;3;1;1;2]);

sol=bvp4c(@ode,@bc,solinit);y=sol.y;

time=y(5)*sol.x;ut=-y(4,:);

%

%ODE’sofaugmentedstatesfunctiondydt=ode(t,y)

dydt=y(5)*[y(2);-y(4);0;-y(3);0];

%

%boundaryconditions:x1(0)=1;x2(0)=2,x1(tf)=3,p2(tf)=0;

3fdenotestheODE’sforstateandcostates.

% p1(tf)*x2(tf)-0.5*p2(2)^2functionres=bc(ya,yb)

res=[ya(1)-1;ya(2)-2;yb(1)-3;yb(4);yb(3)*yb(2)-0.5*yb(4)^2];

Alternatively,wecanaccomplishthisbytreatingTasaparameterforbvp4c[4]Thedifferencebetweenthetwoliesintheparameterlistofbvpinit,andfunctiondefinitionoftheODE’sandboundaryconditions.Figure8showstheresultfromnumericalmethodwhichisthesameastheanalyticalsolution.

solinit=bvpinit(linspace(0,1),[2;3;1;1],2);

sol=bvp4c(@ode,@bc,solinit);y=sol.y;

time=sol.parameters*sol.x;ut=-y(4,:);

%

%ODE’sofaugmentedstatesfunctiondydt=ode(t,y,T)

dydt=T*[y(2);-y(4);0;-y(3)];

%

%boundaryconditions:x1(0)=1;x2(0)=2,x1(tf)=3,p2(tf)=0;

% p1(tf)*x2(tf)-0.5*p2(2)^2functionres=bc(ya,yb,T)

res=[ya(1)-1;ya(2)-2;yb(1)-3;yb(4);yb(3)*yb(2)-0.5*yb(4)^2];

Nowitisthetimetotalkaboutthelimitationsofbvp4c.Thoughtitworkswellsofar,wemustbearinmindthatthequalityofthesolutionfrombvp4cisheavilydependentontheinitialguess.Abadinitialguessmayresultininaccuratesolutions,ornosolutions,orsolutionswhichmakenosense.Forexample,youmaytrytheinitialguessp1(0)=0andcomparetheresultswiththeanalyticalsolutiontoseethedifference.BylookingattheODE’sclosely,wefindthatp˙1(t)=0.Whensuppliedwithaninitialguessof0,bvp4cfailstofindthesolutionduetothestatesingularity.

Thelastproblemistofurtherillustratethecapability,aswellasthelimitation,ofbvp4cinsolvingoptimalcontrolproblems.Thisexamplecanbefindfrom[6]whichisatime-optimalproblemforthedoubleintegralsystemwithasimplecontrolconstraint.

Figure8:Example3Numericalmethod

x˙1(t)

=x2(t)

(10)

x˙2(t)

=u(t)

(11)

Minimizethefinaltime:

T=

dt

0

todrivethestatestotheorigin:

Example4Givenadoubleintegralsystemas:

/tf

f

x(0)=［22］T,x(t)=［00］T,|u|≤1

4

2

μx2+x2

?

3

ToformulatethisproblemasaBVP,controldomainsmoothingtechniquemustbeappliedfirst[6].TheresultedODE’sarelistedbelow:

1

5

x˙ =x(x

+ μx3

\ (12)

2

5

μx2+x2

x˙ =x(? x4

3

\ (13)

4

x˙3 =0 (14)

x˙4

=

?x5x3

(15)

x˙5

=

0

(16)

Theboundaryconditionsare:

x1(0)=2,x2(0)=2,x1(1)=0,x2(1)=0,x3(1)2+x4(1)2=1.

3

4

Auxiliaryvariablex5=T,andu=x4/?μx2+x2.Itiswellknownthat

theoptimalcontrolforthisproblemisthepiecewisecontinuous“bang-bang”control,whichmeansthatuwilltakeeitherthemaximumortheminimumvaluewithintheadmissiblecontrolset,i.e.+1or?1inthisproblem.Becauseitisalinearsecondordersystem,therewillbeonlyoneswitchpointanduisnotcontinuousattheswitchpoint.TheMATLABcodeisasfollows

globalmu;mu=0.5;

solinit=bvpinit(linspace(0,1,10),[2;2;0;-1],4);

sol=bvp4c(@ode,@bc,solinit);

...

%thesolutionforonevalueofmuisusedasguessforthenext.fori=2:3

ifi==2

mu=0.3;

else

mu=.1;end;

%Aftercreatingfunctionhandles,thenewvalueofmu

%willbeusedinnestedfunctions.sol=bvp4c(@ode,@bc,sol);

...

end

%

functiondydt=ode(t,y,T)globalmu;

term=sqrt(mu*y(3)^2+y(4)^2);dydt=T*[y(2)+mu*y(3)/termy(4)/term

0

-y(3)];

%

%boundaryconditions,with4statesand1parameters,5conditionsare

%needed:x1(0)=2,x2(0)=2;x1(1)=0;x2(1)=0;x3(1)^2+x4(1)^2=1;

functionres=bc(ya,yb,T)res=[ya(1)-2

ya(2)-2

Figure9:Example4Statetrajectory

yb(1)

yb(2)yb(3)^2+yb(4)^2-1];

Asμ→0,uappearstobemoreandmorestiffattheswitchpoint,approximatingthe“bang-bang”control.Threevaluesweretested:μ=0.5,μ=0.3,μ=0.1.Weusescontinuationtosolvethisproblem,bywhichmeansthesolutionforμ=0.5isusedasguessfortheBVPwithμ=0.3.ThetrajectoryandcontrolcorrespondingtodifferentμvaluesareshowninFigure9andFigure10.

Wecanclearlyseethatuisapproxi