How to compute Hamiltonian flows and trajectories

In this tutorial, we explain the Flow function, in particular to compute flows from a Hamiltonian vector fields, but also from general vector fields.

Introduction

Consider the simple optimal control problem from the basic example page. The pseudo-Hamiltonian is

\[ H(x, p, u) = p_q\, v + p_v\, u + p^0 u^2 /2,\]

where $x=(q,v)$, $p=(p_q,p_v)$, $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum principle, the maximising control is given in feedback form by

\[u(x, p) = p_v\]

since $\partial^2_{uu} H = p^0 = - 1 < 0$.

u(x, p) = p[2]
nothing # hide

<< @example-block not executed in draft mode >>

Actually, if $(x, u)$ is a solution of the optimal control problem, then, the Pontryagin maximum principle tells us that there exists a costate $p$ such that $u(t) = u(x(t), p(t))$ and such that the pair $(x, p)$ satisfies:

\[\begin{array}{l} \dot{x}(t) = \displaystyle\phantom{-}\nabla_p H(x(t), p(t), u(x(t), p(t))), \\[0.5em] \dot{p}(t) = \displaystyle - \nabla_x H(x(t), p(t), u(x(t), p(t))). \end{array}\]

Let us import the necessary packages.

using OptimalControl
using OrdinaryDiffEq

<< @example-block not executed in draft mode >>

The package OrdinaryDiffEq.jl provides numerical integrators to compute solutions of ordinary differential equations.

Extremals from the Hamiltonian

The pairs $(x, p)$ solution of the Hamitonian vector field are called extremals. We can compute some constructing the flow from the optimal control problem and the control in feedback form. Another way to compute extremals is to define explicitely the Hamiltonian.

H(x, p, u) = p[1] * x[2] + p[2] * u - 0.5 * u^2     # pseudo-Hamiltonian
H(x, p) = H(x, p, u(x, p))                          # Hamiltonian

z = Flow(Hamiltonian(H))

t0 = 0
tf = 1
x0 = [-1, 0]
p0 = [12, 6]
xf, pf = z(t0, x0, p0, tf)

<< @example-block not executed in draft mode >>

Extremals from the Hamiltonian vector field

You can also provide the Hamiltonian vector field.

Hv(x, p) = [x[2], p[2]], [0.0, -p[1]]     # Hamiltonian vector field

z = Flow(HamiltonianVectorField(Hv))
xf, pf = z(t0, x0, p0, tf)

<< @example-block not executed in draft mode >>

Note that if you call the flow on tspan=(t0, tf), then you obtain the output solution from OrdinaryDiffEq.jl.

sol = z((t0, tf), x0, p0)
xf, pf = sol(tf)[1:2], sol(tf)[3:4]

<< @example-block not executed in draft mode >>

Trajectories

You can also compute trajectories from the control dynamics $(x, u) \mapsto (v, u)$ and a control law $t \mapsto u(t)$.

u(t) = 6-12t
x = Flow((t, x) -> [x[2], u(t)]; autonomous=false) # the vector field depends on t
x(t0, x0, tf)

<< @example-block not executed in draft mode >>

Again, giving a tspan you get an output solution from OrdinaryDiffEq.jl.

using Plots
sol = x((t0, tf), x0)
plot(sol)

<< @example-block not executed in draft mode >>