# Problem with optimal control and Pontryagin's maximum principle

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}$$

1. 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[\[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, \[Lambda] == 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.

• There is a typo in your code with J and H definition. It should be H = D[f, x[t]]^2 + u^2 + [Lambda][t] eqn May 9, 2021 at 17:11
• @AlexTrounev I corrected this error.
– dtn
May 9, 2021 at 17:27

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[\[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 == 1, \[Lambda] == 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["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"]; \
Get["NumericalDifferentialEquationAnalysis"];

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

fd = NDSolveFiniteDifferenceDerivative[Derivative, 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 == 0, x == 1};


Solution

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


Visualization

lst1 = Table[{tgrid[[i]], x[i] /. sol[]}, {i,
Length[tgrid]}]; lst2 =
Table[{tgrid[[i]], u[i] /. sol[]}, {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$$.

• The system loses performance if $x(0)> 1.217$. This is associated with the choice $\lambda(0)$.
– dtn
May 10, 2021 at 4:21
• How will the original Hamiltonian look like if it adds the criterion for the minimum convergence time of $\frac{df}{dx}$ to $0$ ?
– dtn
May 10, 2021 at 4:27
• @Alex Tournev: In the first part of your code, I get General::ivar: -(\[Lambda][t]/2) is not a valid variable. And in the second part of your code (Update 1), I receive NMinimize::ivar: -(\[Lambda]/2) is not a valid variable. and a bunch of errors relating to Lamda. Why do not I get the solution you have? Any idea? May 11, 2021 at 16:52
• @TugrulTemel Don't run code 1 and code 2 together. Put Clear["Global*"] before code 1 and code 2. May 11, 2021 at 17:56
• @Vajira Use in[t_] := Evaluate[(D[f, x[t]]^2 + u[t]^2) /. sys], and then evaluate NIntegrate[in[t], {t, 0, 10}] with result {0.694438}. Jan 21, 2022 at 7:38