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]

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

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)

([2.134348517083835e-15, 2.3650010008626292e-14], [12.0, -5.999999999999998])

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)

([2.134348517083835e-15, 2.3650010008626292e-14], [12.0, -5.999999999999998])

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]

([2.152733210148547e-15, 2.6983767684360683e-14], [12.0, -5.999999999999992])

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)

2-element Vector{Float64}:
 -8.881784197001252e-16
 -4.440892098500626e-16

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

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