Problem with optimal control and Pontryagin's maximum principle

Question

For dynamic system:

$\dot{x}=\frac{df}{dx}+u$

where $f=e^{-x^2}$

It is necessary to develop optimal control, minimizing criterion:

$J= \int_{0}^{t_f} ((\frac{df}{dx})^2+u^2) \,dt $

Algorithm:

1. We write Hamiltonian: $H=((\frac{df}{dx})^2+u^2)+\lambda (\frac{df}{dx}+u)$

2. Write costate equation: $\dot{\lambda}=-\frac{dH}{d\lambda}$

3. Solve equation for control signal $u$: $\frac{dH}{du}=0$

4.Write resulting system of equation:$\begin{cases} \dot{x}=... \\ \dot{\lambda}=... \end{cases}$

5. Solve numerically:

I wrote this algorithm to Mathematica. There is my code:

(***)

Clear["Derivative"]

ClearAll["Global`*"]

f = Sech[(x[t] - 2)]

(***-Origin ODE)

eqn = -D[f, x[t]] + u

(***J)

J = Integrate[D[f, x[t]]^2 + u[t]^2, {t, 0, tf}]

(***Hamiltonian)

H = D[f, x[t]]^2 + 2 u^2 + \[Lambda][t] eqn

(***Costate-Equation)

cseqn = Derivative[1][\[Lambda]][t] == -D[H, x[t]]

(***Solution-For-Control-Signal)

Solve[D[H, u] == 0, u]

u = -\[Lambda][t]/4

(***Resulting-system-of-equation)

eqns = {x'[t] == D[f, x[t]] + u, cseqn}

sys = NDSolve[{eqns, 
   x[0] == 0, \[Lambda][0] == 0}, {x, \[Lambda]}, {t, 0, 150}]

Plot[{Evaluate[x[t] /. sys], 2}, {t, 0, 30}, PlotRange -> Full, 
 PlotPoints -> 100]

Plot[{Evaluate[D[f, x[t]] /. sys]}, {t, 0, 100}, PlotRange -> Full, 
 PlotPoints -> 100]

My problems:

1. System does not come to a state $\frac{df}{dx}=0$
2. I do not know how the minimized / maximized criterion $J$ is formed correctly (at an infinite time interval).

In other words, the system should achieve a state $\frac{df}{dx}=0$ as quickly as possible.

I would be grateful to help in awareness and correcting my mistakes.

Answer 1

Pontryagin's minimum principle means that we have to use Euler-Lagrange equations. Therefore code looks like this

ClearAll["Global`*"]

f = Exp[-x[t]^2];

(*Origin ODE *) eqn = D[f, x[t]] + u[t];
 J = Integrate[D[f, x[t]]^2 + u[t]^2, {t, 0, tf}];
H = D[f, x[t]]^2 + u[t]^2 + \[Lambda][t] eqn;
(*Costate-Equation*) 
cseqn = Derivative[1][\[Lambda]][t] == D[H, x[t]];
(*Solution-For-Control-Signal*) Solve[D[H, u[t]] == 0, u[t]] ;

u[t_] := -(\[Lambda][t]/2)

(*Resulting-system-of-equation*) eqns = {x'[t] == D[f, x[t]] + u[t],
   cseqn} ; sys = 
 NDSolve[{eqns, x[0] == 1, \[Lambda][0] == 0}, {x, \[Lambda]}, {t, 0, 
   10}]; {Plot[{Evaluate[x[t] /. sys]}, {t, 0, 10}, PlotRange -> Full,
   PlotPoints -> 100], 
 Plot[{Evaluate[(-\[Lambda][t]/2) /. sys]}, {t, 0, 10}, 
  PlotRange -> Full, PlotPoints -> 100]}

Update 1. In a case of additional constrains we can use FDM and NMinimize[] as a solver. For example, the problem with control considered above can be solved as optimization problem as follows

Clear["Global`*"]
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"]; \
Get["NumericalDifferentialEquationAnalysis`"];

L = 10; g = GaussianQuadratureWeights[100, 0, L];
ugrid = g[[All, 1]]; weights = g[[All, 2]]; tgrid = Join[{0}, ugrid];


fd = NDSolve`FiniteDifferenceDerivative[Derivative[1], tgrid]; m = 
 fd["DifferentiationMatrix"]; varu = 
 Table[u[i], {i, Length[tgrid]}]; varx = 
 Table[x[i], {i, Length[tgrid]}]; xt = m . varx; int = 
 Sum[weights[[i]] (u[i]^2 + (2 x[i] Exp[-x[i]^2])^2), {i, 
   Length[weights]}];

eqns = Table[
   xt[[i]] - (-2 x[i] Exp[-x[i]^2] + u[i]) == 0, {i, Length[xt]}];
ics = {u[1] == 0, x[1] == 1};

Solution

sol = NMinimize[{int, Join[eqns, ics]}, Join[varu, varx]];

Visualization

lst1 = Table[{tgrid[[i]], x[i] /. sol[[2]]}, {i, 
   Length[tgrid]}]; lst2 = 
 Table[{tgrid[[i]], u[i] /. sol[[2]]}, {i, Length[tgrid]}];

{ListLinePlot[lst1, AxesLabel -> {"t", "x"}], 
 ListLinePlot[lst2, AxesLabel -> {"t", "u"}]}

Note that Figure 2 looks different then Figure 1. Also numerical result for the first solution is $J=0.694438$, and for the second one $J=0.582044$.

Answer 2

Preamble: I think there are some formulation questions that I have about the original question, but I'd like to take crack and show my work. This way boundary conditions can be refined. I believe there's lots of vagueness in learning materials leading to trouble defining the terminal boundary conditions.

Also, I took some liberties from dtn's code to sketch out my solution, e.g. initial conditions:

Time parameters:

t0 = 0;
tf = 150;

Declare variables:

z = {
   x[t]
   };
zdt = {
   x'[t]
   };
U = {
   u[t]
   };

Derivation of the Lagrangian of the Hamiltonian:

f = {
   D[Exp[-x[t]^2], x[t]] + u[t]
   };
F = (D[Exp[-x[t]^2], x[t]])^2 + u[t]^2;
M = 0;
\[Lambda]z = Map[Indexed[\[Lambda], Head[#]][t] &, z];
\[Lambda]zdt = Map[Indexed[\[Lambda], Head[#]]'[t] &, z];
Fz = Transpose[Grad[{F}, z]];
fz = Transpose[Grad[f, z]];
H = F + Transpose[\[Lambda]z] . f + M;
L = H;

Dermine optimal control:

Lsol = Minimize[L, U];
Ustar = Last[Lsol];

Re-form state and costate equations from Lagranian with optimal U.

Lstar = L /. Ustar;
zeqn = Thread[zdt == Grad[Lstar, \[Lambda]z]];
\[Lambda]zeqn = Thread[\[Lambda]zdt == -Grad[Lstar, z]];

Declare boundary conditions and endpoint term:

bc = {
   x[t0] == 1
   };
tc = Thread[(\[Lambda]z /. t -> tf) == Grad[M, z]];

Equations of (co)state:

  Join[
   zeqn,
   \[Lambda]zeqn,
   bc,
   tc
   ] // FullSimplify // TableForm

Solve with numerical ODE solver:

ndsol = NDSolve[
  Join[
   zeqn,
   \[Lambda]zeqn,
   bc,
   tc
   ]
  ,
  Join[
   z,
   \[Lambda]z
   ],
  {t, t0, tf}
  ];

Plot results:

Plot[z /. ndsol, {t, t0, tf}]
Plot[\[Lambda]z /. ndsol, {t, t0, tf}]
Plot[U /. (Ustar /. ndsol), {t, t0, tf}]

or zooming in to (0,10):

Addendum: The procedure outlined by Alex uses NMinimize. However, this type of problem can be solved more quickly using IPOPT. The following code assumes you've run the code in Alex's answer.

Needs["IPOPTLink`"]
objfn = int;
vars = Join[varu, varx];
x0 = Join[
   ConstantArray[0, Length[varu]],
   Join[
    {1}, (* x[1] *)
    ConstantArray[0, Length[varx] - 1]
    ]
   ];
varbounds = Join[
   ConstantArray[{-Infinity, Infinity}, Length[varu]],
   ConstantArray[{-Infinity, Infinity}, Length[varx]]
   ];
constr = Join[
   eqns[[All,1]],
   {u[1]},
   {x[1]}
   ];
constrbounds = Join[
   ConstantArray[{0, 0}, Length[xt]],
   {{0, 0}}, (* u[1] == 0 *)
   {{1, 1}}  (* x[1] == 1 *)
   ];
solip = IPOPTMinimize[objfn, vars, x0, varbounds, constr, constrbounds];
ipoptsol = MapThread[Rule, {vars, IPOPTArgMin[solip]}];
lst1 = Table[{tgrid[[i]], x[i] /. ipoptsol}, {i, Length[tgrid]}];
lst2 = Table[{tgrid[[i]], u[i] /. ipoptsol}, {i, Length[tgrid]}];

{
 ListLinePlot[lst1, AxesLabel -> {"t", "x"}, Mesh -> All],
 ListLinePlot[lst2, AxesLabel -> {"t", "u"}, Mesh -> All]
 }