Constrained 1d optimal control (Prob. 7.15 / Ex. 7.5)

• revised Oct. 26, 2021
	– improve text formatting (xdot -> (OverDot[x])) and clean  up other comments

Clear["Global`*"];

Optimal control problem has running cost    L=(1/2)((x^2)+(u^2))    integrated over t ∈(0,∞).  
 	Eq. of motion is         (OverDot[x]) = -x + u,     x(0) = x0.

Unconstrained case, leads to u(t) = -K* x(t), with K* = (Sqrt[2])-1.

j[τ_,xτ_]:=Module[{eq,init,x,λ,t,xs,λs,us,hs},
eq={
x'[t]==- x[t]-λ[t], 
λ'[t]==  λ[t]-x[t]
};
init={x[τ]==xτ,λ[∞]==0};
{xs,λs}=DSolveValue[{eq,init},{x,λ},t]//FullSimplify;
us[t_]:=-λs[t];
hs[t_]:=(1/2)(((xs[t]-λs[t]))^2)-((λs[t])^2)//FullSimplify;
{xs,λs,us,hs}]
{xu,λu,uu,hu}=j[τ,xτ];
Grid[{{
xu[t],Kstar=((uu[t])/(-xu[t]))//FullSimplify
}},Spacings->4]

Plot[{xu[t],λu[t],uu[t]}/.{τ->1,xτ->1},{t,0,10},PlotRange->All]

Constrained case

jc[τ_,x0_]:=Module[{eqx,initx,eqλ,initλ,x,λ,t,xs,λs,us,hs},
eqx=x'[t]==- x[t]-1;
initx=x[0]==x0;
xs=DSolveValue[{eqx,initx},x,t];

eqλ=λ'[t]==  λ[t]-xs[t];
initλ=λ[τ]==((Sqrt[2])-1)xs[τ];  (* wrong condition??  λ not continuous ?  *)
λs=DSolveValue[{eqλ,initλ},λ,t]//FullSimplify;

us[t_]:=-1;
hs[t_]:=(1/2)(((xs[t])^2)+1)-λs[t]xs[t]-λs[t]//FullSimplify;
{xs,λs,us,hs}]
{xc,λc,uc,hc}=jc[τ,x0];
{xc[t],λc[t],uc[t],hc[t]}  (* solution for λ(t) is probably wrong, since final condition may be wrong, but important part is ok *)

hc1=hc[t]/.x0->5;
Plot[hc1,{τ,0,2}]

FindMinimum[hc1,{τ,1}]

FindRoot[hc1,{τ,1}]

Compute overall cost function jtot

j1[τ_]:=Integrate[(1/2)((xc[t])^2),{t,0,τ}]+(τ/2)//FullSimplify;j1[τ]

(1/2)(1+(((-1+(Sqrt[2])))^2))//Simplify

xc[τ]

j2[τ_]:=Integrate[(2-(Sqrt[2]))((xu[t])^2)/.xτ->xc[τ],{t,τ,∞}]//Simplify;  j2[τ]

jtot[τ_]:=j1[τ]+j2[τ]//Simplify;jtot[τ]

jtot1=jtot[τ]/.x0->5//FullSimplify

Grid[{{Plot[jtot1,{τ,0.55,0.58}],Log[((1+5)/(2+(Sqrt[2])))]//N,FindMinimum[{jtot1,τ≥0.5},{τ,1}]}},Spacings->2]

Simplify/@{(∂_τ)jtot1, (∂_τ)jtot1/.τ->Log[((1+5)/(2+(Sqrt[2])))]}  (* check the minimization *)

So the minimum (inflection pt., actually) occurs just at the crossover point where u = –1, 
as one might expect.

Now do it numerically the “right” way, with a single equation, using u = -Clip[λ].  
Unfortunately, as DSolveValue does not seem to work with the piecewise ODE 
(although it claims it can...), I do it numerically.  
But the numerical solution confirms the above analysis.

j3[x0_]:=Module[{eq,init,x,λ,t,xs,λs,us,τ},
τ=10;
eq={
x'[t]==- x[t]-Clip[λ[t]],   (* value between ±1 *)
λ'[t]==  λ[t]-x[t]
};
init={x[0]==x0,λ[τ]==0};
{xs,λs}=NDSolveValue[{eq,init},{x,λ},{t,0,τ}]//FullSimplify;
us[t_]:=-Clip[λs[t]];
{xs,λs,us}]
{x,λ,u}=j3[5];

Plot[{x[t],λ[t],u[t]},{t,0,10},PlotRange->All]

hc[τ]/.τ->Log[((1+x0)/(2+(Sqrt[2])))]//Simplify

((3/2)-(Sqrt[2]))(((2+(Sqrt[2])))^2)//Simplify