10
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$ May 9, 2021 at 17:11
  • $\begingroup$ @AlexTrounev I corrected this error. $\endgroup$
    – dtn
    May 9, 2021 at 17:27
  • $\begingroup$ 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? $\endgroup$
    – Eric Brown
    Apr 18, 2023 at 15:07
  • $\begingroup$ @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/… $\endgroup$
    – dtn
    Apr 19, 2023 at 6:50
  • $\begingroup$ @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. $\endgroup$
    – Eric Brown
    Jun 18, 2023 at 0:41

2 Answers 2

12
$\begingroup$

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

Figure 1

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"}]}

Figure 2

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

$\endgroup$
18
  • $\begingroup$ The system loses performance if $x(0)> 1.217$. This is associated with the choice $\lambda(0)$. $\endgroup$
    – dtn
    May 10, 2021 at 4:21
  • $\begingroup$ How will the original Hamiltonian look like if it adds the criterion for the minimum convergence time of $\frac{df}{dx}$ to $0$ ? $\endgroup$
    – dtn
    May 10, 2021 at 4:27
  • 1
    $\begingroup$ @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? $\endgroup$ May 11, 2021 at 16:52
  • 2
    $\begingroup$ @TugrulTemel Don't run code 1 and code 2 together. Put Clear["Global`*"] before code 1 and code 2. $\endgroup$ May 11, 2021 at 17:56
  • 2
    $\begingroup$ @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}. $\endgroup$ Jan 21, 2022 at 7:38
5
$\begingroup$

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

ustar

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

equations of (co)state and boundary conditions

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

x(t) lambda_x(t) ustar

or zooming in to (0,10): x(t) zoomed lambda_x(t) zoomed ustar zoomed

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]
 }
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 13, 2023 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.