I continue to study the topic I started here: Problem with optimal control and Pontryagin's maximum principle

A simple ODE system $(1)$ is given:

$F=\begin{cases} \dot{x}=g \\ \dot{g}=-g+\frac{df}{dx}+u \\ \dot{h}=-h+\frac{d^2f}{dx^2} \end{cases} (1)$


$x,g,h$ - state-space variables


Purpose: compose cost function $$J= \int_{0}^{\infty} (?+g+u^2)dt$$ that minimize rate of transition to steady-state, using state-space variables $x,g,h,\dot{x},\dot{g},\dot{h}$.

A feature of the original system is the transition to a state in which $g=0$ and $h$ is negative.

I did it according to the same algorithm as last time, but in contrast to the previous case, as a method of acceleration, I chose the gain factor $k = 5$ in cost function:


  1. We write Hamiltonian: $H=-(k \cdot (g+u^2))+\boldsymbol{\lambda} \cdot F$ where $\boldsymbol{\lambda}=(\lambda_1,\lambda_2,\lambda_3)$

  2. Write costate equation: $\dot{\lambda}=\frac{dH}{d\boldsymbol{X}}$ where $\boldsymbol{X}=(x,g,h)$

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

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

  5. Solve numerically.

Just like in past, I am attaching my own version of the code in Mathematica:


x0 = {{1}, {0}, {0}};

\[Lambda]0 = {{0}, {0}, {1/4}};

f = Exp[-(x[t] + 1)^2];

J = Integrate[5 (g[t] + u[t])^2, {t, 0, Infinity}]

(*Origin ODE*)eqn = {{g[t]}, {-g[t] + D[f, x[t]] + u[t]}, {-h[t] + 
    D[f, {x[t], 2}]}}

L[t_] = 5 g[t] + 5 u[t]^2

lambda[t_] := {{\[Lambda]1[t]}, {\[Lambda]2[t]}, {\[Lambda]3[t]}}

X[t_] := {{x[t]}, {g[t]}, {h[t]}}

H[t_] = Flatten[-L[t] + lambda[t]\[Transpose].eqn][[1]]

uSol = First@Solve[0 == D[H[t], u[t]], u[t]]

 eqn1 = Table[
   D[lambda[t][[i, 1]], t] == D[H[t] /. uSol, X[t][[i, 1]]], {i, 1, 

 eqn2 = Table[
   D[X[t][[i, 1]], t] == D[H[t] /. uSol, lambda[t][[i, 1]]], {i, 1, 

bcx = Table[X[0][[i, 1]] == x0[[i, 1]], {i, 1, 3}]

bc\[Lambda] = Table[lambda[0][[i, 1]] == \[Lambda]0[[i, 1]], {i, 1, 3}]

Flatten[{eqn1, eqn2, bcx, bc\[Lambda]}]

sys = NDSolve[
  Flatten[{eqn1, eqn2, bcx, bc\[Lambda]}], {x[t], g[t], 
   h[t], \[Lambda]1[t], \[Lambda]2[t], \[Lambda]3[t]}, {t, 0, 500}]

Plot[{Evaluate[x[t] /. sys]}, {t, 0, 25}, PlotRange -> Full, 
 PlotPoints -> 100, ImageSize -> Small]

Plot[{Evaluate[g[t] /. sys]}, {t, 0, 20}, PlotRange -> Full, 
 PlotPoints -> 100, ImageSize -> Small]

Plot[{Evaluate[h[t] /. sys]}, {t, 0, 20}, PlotRange -> Full, 
 PlotPoints -> 100, ImageSize -> Small]

Plot[{Evaluate[L[t] /. uSol /. sys]}, {t, 0, 100}, PlotRange -> All, 
 ImageSize -> Small]

Plot[{Evaluate[u[t] /. uSol /. sys]}, {t, 0, 50}, PlotRange -> All, 
 ImageSize -> Small]


  1. How can a cost function $J$ be compiled to minimize the duration of transients over an infinite time interval.
  2. Should we use derivatives of state variables $\dot{x},\dot{g},\dot{h}$ for this? If so, than how?
  3. How does the initial state of the co-state vector $\boldsymbol{\lambda}$ affect this?

In general terms, the problem sounds like this: how the optimal control synthesized according to the Pontryagin maximum/minimum principle tuned?

I would be grateful to help.

  • 1
    $\begingroup$ We don't need Infinity limit for the numerical problem. Just put limit time in Integrate the same as in NDSolve. $\endgroup$ – Alex Trounev May 17 at 11:49
  • $\begingroup$ @AlexTrounev Thank you for your reply. The goal is to develop a control system and transfer it to simulink. That is, the structure of the control signal, subsystem of co-states and initial values. NDSolve in this case is used as an auxiliary tool for checking and assessing the effect of settings on transients. $\endgroup$ – dtn May 17 at 11:53
  • 1
    $\begingroup$ Your question #1 about duration can be reformulated for the interval used in NDSolve with using norm $|\bf X-X_f|$, where $\bf X_f$ is final state. $\endgroup$ – Alex Trounev May 17 at 12:30
  • $\begingroup$ @AlexTrounev and if the final state of variable $x$ is unknown. By the way, this is one of the features of the system, I forgot about it. $\endgroup$ – dtn May 17 at 12:32
  • 1
    $\begingroup$ We can use as a final state the solution of equation $\bf \dot {X}=0$ at $u=0$. Since we minimize $J=\int {(u^2+...)dt}$ we can suggest that u^2->0 at t->Infinity. $\endgroup$ – Alex Trounev May 17 at 17:51

The following comment is about the application of the PMP, not about Mathematica per se.

The Pontryagin maximum principle (PMP) is not designed for application on an infinite interval. The PMP provides a set of first-order necessary conditions for optimality in an optimal control problem, which in turn leads to a set of well-posed boundary value problems. For instance, if the optimal control problem is (and stackexchange is forcing me to treat the following LaTeX code as code):

& \text{minimize}_u && \int_0^T c(x(t), u(t))\, dt + c_F(x(T))\\
& \text{subject to} && 
\dot x(t) = f(x(t), u(t)),\\
x(0) = x_0,\; x(T) \text{ unspecified},

with standard assumptions on the instantaneous cost $c$, the system $f$, the terminal cost $c_F$, etc., then assuming that $(x_\star, u_\star)$ is an optimal solution of the problem, the PMP asserts that there exist

  • an abnormal multiplier $\eta\in\{0, 1\}$ and

  • an adjoint equation expressed in terms of the Hamiltonian $H^\eta(p, x, u) := \langle p, f(x, u) \rangle - \eta c(x, u)$ given by

    $$ \dot p(t) = - \frac{\partial H^\eta}{\partial x}(p(t), x_\star(t), u_\star(t)) $$

with boundary (final) condition $p(T) = -\frac{\partial c_F}{\partial x}(x_\star(T))$, together with the nontriviality condition $(\eta, p)\neq (0, 0)$.

There are cases where one first applies the PMP for finite $T$ and then argues via other means what happens to the optimal state-action trajectories if $T\to +\infty$; one such notable example is the linear quadratic regulator problem. Anyway, I digress.

The key takeaway from the PMP is that it leads to a two-point boundary value problem, which is then solved via, e.g., shooting via root-finding techniques, or collocation methods, etc.

In your problem, the first part of the "Purpose" seems to be suspect because you're looking at $[0, +\infty[$. Also, your cost has $g$ in it; what prevents the minimization from carrying $g$ to $-\infty$ (i.e., the problem being ill-posed)?

Sorry if I'm not being helpful with Mathematica.

  • $\begingroup$ I have only recently begun to study this topic and my knowledge in this area is very weak, but I seem to have caught what you are trying to say. I will give my comments. botik.ru/PSI/CPRC/sachkov/notes.pdf In this book, Chapter 13 "" has two interesting problems. 13.1 and 13.2 (fastest and cheapest stop of a train). $\endgroup$ – dtn May 18 at 7:41
  • $\begingroup$ If I understand correctly, these problems are formulated in such a way that the original dynamical system is transformed into a Hamiltonian system that operates on an infinite time interval. And it performs the functions inherent in it. I am guided by this. I also found such an article, which shows typical components of the cost function that perform one or another function (following a trajectory, minimizing a control signal, etc.) motion.cs.illinois.edu/RoboticSystems/OptimalControl.html Do you understand what I am talking about ? $\endgroup$ – dtn May 18 at 7:41
  • $\begingroup$ In my understanding, optimal control is an extension of the original dynamical system to the Hamiltonian one, which is a set of initial differential equations and co-state equations. Their interconnection leads to the structure of the control signal that minimizes the Hamiltonian. Moreover, this minimization should occur continuously from the initial values to the final ones (in my case, determined by the fundamental properties of the original system of differential equations). $\endgroup$ – dtn May 18 at 7:44
  • $\begingroup$ @dtn: It's not quite an accurate assessment to think that optimal control is an extension of the original controlled dynamical system to a Hamiltonian one; please see math.rutgers.edu/~sussmann/papers/brockettpaper.ps.gz for a detailed explanation of the differences, including between classical calculus of variations and optimal control. $\endgroup$ – dchatter May 18 at 12:41
  • 1
    $\begingroup$ @dtn: Including the square of the absolute value of the derivatives of the velocities, for instance, makes sense. I suggest taking a look at the LQR problem (Chapter 6 in Agrachev-Sachkov). $\endgroup$ – dchatter May 18 at 16:07

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.