# 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[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.

• 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 Commented May 9, 2021 at 17:11
• @AlexTrounev I corrected this error.
– ayr
Commented May 9, 2021 at 17:27
• The more I study this problem, the more I think there's something missing, such as the value of Lambda[Tfinal=150], i.e. Termination conditions. Aren't PMP problems two-point BVP's? Commented Apr 18, 2023 at 15:07
• @EricBrown As far as I know, classical optimal control was really considered on finite time intervals. The problem where the end time is not explicitly set is called "optimal control on an infinite horizon" math.stackexchange.com/questions/2178582/…
– ayr
Commented Apr 19, 2023 at 6:50
• @dtn thanks for the clarification. Indeed your function and its derivative decay to 0 pretty quickly. Also my understanding is that lambda(infinity) = 0 in infinite horizon problems. Question I have is: since u is not constrained, is it possible to choose u so that the function practically never decays? q.v. my answer below. Commented Jun 18, 2023 at 0:41

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["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[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$$.

• The system loses performance if $x(0)> 1.217$. This is associated with the choice $\lambda(0)$.
– ayr
Commented 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$ ?
– ayr
Commented 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][1]/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? Commented 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. Commented 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}. Commented Jan 21, 2022 at 7:38

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];
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;


Declare boundary conditions and endpoint term:

bc = {
x[t0] == 1
};


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];

• Thank you for IPOPT method implementation to solve this problem (+1). What is a difference to optimize the problem with NDSolve and NMinimize? In the first case we are looking for stationary points to functional J using Euler-Lagrange equation. As well known it could be extremum (minimum or maximum). In the second case we are looking for minimum int using NMinimize. This is why minimum J= 0.582044 while extremum J=0.694438`. Commented Apr 13, 2023 at 15:17