6
$\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$
2
  • $\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$ – Alex Trounev May 9 at 17:11
  • $\begingroup$ @AlexTrounev I corrected this error. $\endgroup$ – dtn May 9 at 17:27
8
$\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$
8
  • $\begingroup$ The system loses performance if $x(0)> 1.217$. This is associated with the choice $\lambda(0)$. $\endgroup$ – dtn May 10 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 at 4:27
  • $\begingroup$ @dtn Could you explain what kind of solution your are looking for? $\endgroup$ – Alex Trounev May 10 at 15:23
  • 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$ – Tugrul Temel May 11 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$ – Alex Trounev May 11 at 17:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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