How to compute flows

In this tutorial, we explain the Flow function from OptimalControl.jl package.

Basic usage

Les us define a basic optimal control problem.

using OptimalControl

t0 = 0
tf = 1
x0 = [-1, 0]

ocp = @def begin

    t ∈ [ t0, tf ], time
    x = (q, v) ∈ R², state
    u ∈ R, control

    x(t0) == x0
    x(tf) == [ 0, 0 ]

    ẋ(t)  == [ v(t), u(t) ]

    ∫( 0.5u(t)^2 ) → min

end

The pseudo-Hamiltonian of this problem is

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

where $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}\]

The Flow function aims to compute $(x, p)$ from the optimal control problem ocp and the control in feedback form u(x, p).

Let us try to get the associated flow:

julia> f = Flow(ocp, u)ERROR: ExtensionError. Please make: julia> using OrdinaryDiffEq

As you can see, an error occured since we need the package OrdinaryDiffEq.jl. This package provides numerical integrators to compute solutions of the ordinary differential equation $\dot{z}(t) = \vec{\mathbf{H}}(z(t))$.

using OrdinaryDiffEq

f = Flow(ocp, u)

Now we have the flow of the associated Hamiltonian vector field, we can use it. Some simple calculations shows that the initial covector $p(0)$ solution of the Pontryagin maximum principle is $[12, 6]$. Let us check that integrating the flow from $(t_0, x_0) = (0, [-1, 0])$ to the final time $t_f$ we reach the target $x_f = [0, 0]$.

p0 = [12, 6]
xf, pf = f(t0, x0, p0, tf)

([2.425225564003491e-15, 2.4789995105012755e-14], [12.0, -5.999999999999998])

If you prefer to get the state, costate and control trajectories at any time, you can call the flow:

sol = f((t0, tf), x0, p0)

In this case, you obtain a data that you can plot exactly like when solving the optimal control problem with the function solve. See for instance the basic example or the plot tutorial.

using Plots

plot(sol)

You can notice from the graph of v that the integrator has made very few steps:

time_grid(sol)

6-element Vector{Float64}:
 0.0
 0.0027652656319165627
 0.01991270125908861
 0.10962307175239358
 0.42412906329190825
 1.0

To have a better visualisation (the accuracy won't change), you can provide a fine grid.

sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 100))
plot(sol)

The argument saveat is an option from OrdinaryDiffEq.jl. Please check the list of common options. For instance, one can change the integrator with the keyword argument alg or the absolute tolerance with abstol. Note that you can set an option when declaring the flow or set an option in a particular call of the flow. In the following example, the integrator will be BS5() and the absolute tolerance will be abstol=1e-8.

f = Flow(ocp, u; alg=BS5(), abstol=1)   # alg=BS5(), abstol=1
xf, pf = f(t0, x0, p0, tf; abstol=1e-8) # alg=BS5(), abstol=1e-8

([9.187823365203511e-17, 1.6903858162793364e-16], [12.0, -6.0])

Extremals and trajectories

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))
xf, pf = z(t0, x0, p0, tf)

([2.425225564003491e-15, 2.4789995105012755e-14], [12.0, -5.999999999999998])

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.425225564003491e-15, 2.4789995105012755e-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.936701702520115e-15, 4.115723153754471e-14], [12.0, -5.999999999999967])

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
  1.5543122344752192e-15

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

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

Variable

Let us consider an optimal control problem with a (decision / optimisation) variable.

t0 = 0
x0 = 0

ocp = @def begin

    tf ∈ R, variable             # the optimisation variable is tf
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control

    x(t0) == x0
    x(tf) == 1

    ẋ(t) == tf * u(t)

    tf + 0.5∫(u(t)^2) → min

end

As you can see, the variable is the final time tf. Note that the dynamics depends on tf. From the Pontryagin maximum principle, the solution is given by:

tf = (3/2)^(1/4)
p0 = 2tf/3

The input arguments of the maximising control are now the state x, the costate p and the variable tf.

u(x, p, tf) = tf * p

Let us check that the final condition x(tf) = 1 is satisfied.

f = Flow(ocp, u)
xf, pf = f(t0, x0, p0, tf, tf)

(1.0000000000000004, 0.7377879464668812)

The usage of the flow f is the following: f(t0, x0, p0, tf, v) where v is the variable. If one wants to compute the state at time t1 = 0.5, then, one must write:

t1 = 0.5
x1, p1 = f(t0, x0, p0, t1, tf)

(0.45180100180492255, 0.7377879464668812)

xf, pf = f(t0, x0, p0, tf)

(1.0000000000000004, 0.7377879464668812)

Since the variable is the final time, we can make the time-reparameterisation $t = t_f \, s$ to normalise the time $s$ in $[0, 1]$.

ocp = @def begin

    tf ∈ R, variable
    s ∈ [0, 1], time
    x ∈ R, state
    u ∈ R, control

    x(0) == 0
    x(1) == 1

    ẋ(s) == tf^2 * u(s)

    tf + (0.5*tf)*∫(u(s)^2) → min

end

f = Flow(ocp, u)
xf, pf = f(0, x0, p0, 1, tf)

(1.0000000000000002, 0.7377879464668812)

Another possibility is to add a new state variable $t_f(s)$. The problem has no variable anymore.

ocp = @def begin

    s ∈ [0, 1], time
    y = (x, tf) ∈ R², state
    u ∈ R, control

    x(0) == 0
    x(1) == 1

    dx = tf(s)^2 * u(s)
    dtf = 0 * u(s) # 0
    ẏ(s) == [dx, dtf]

    tf(1) + 0.5∫(tf(s) * u(s)^2) → min

end

u(y, q) = y[2] * q[1]

f = Flow(ocp, u)
yf, pf = f(0, [x0, tf], [p0, 0], 1)

([1.0000000000000002, 1.1066819197003217], [0.7377879464668812, -1.0000000000000004])

Concatenation of arcs

In this part, we present how to concatenate several flows. Let us consider the following problem.

t0 =  0
tf =  1
x0 = -1
xf =  0

@def ocp begin

    t ∈ [ t0, tf ], time
    x ∈ R, state
    u ∈ R, control

    x(t0) == x0
    x(tf) == xf

    -1 ≤ u(t) ≤ 1

    ẋ(t) == -x(t) + u(t)

    ∫( abs(u(t)) ) → min

end

From the Pontryagin maximum principle, the optimal control is a concatenation of an off arc ($u=0$) followed by a positive bang arc ($u=1$). The initial costate is

\[p_0 = \frac{1}{x_0 - (x_f-1) e^{t_f}}\]

and the switching time is $t_1 = -\ln(p_0)$.

p0 = 1/( x0 - (xf-1) * exp(tf) )
t1 = -log(p0)

Let us define the two flows and the concatenation. Note that the concatenation of two flows is a flow.

f0 = Flow(ocp, (x, p) -> 0)     # off arc: u = 0
f1 = Flow(ocp, (x, p) -> 1)     # positive bang arc: u = 1

f = f0 * (t1, f1)               # f0 followed by f1 whenever t ≥ t1

Now, we can check that the state reach the target.

sol = f((t0, tf), x0, p0)
plot(sol)

For the moment, this concatenation is not equivalent to an exact concatenation.

f = Flow(x ->  x)
g = Flow(x -> -x)

x0 = 1
φ(t) = (f * (t/2, g))(0, x0, t)
ψ(t) = g(t/2, f(0, x0, t/2), t)

println("φ(t) = ", abs(φ(1)-x0))
println("ψ(t) = ", abs(ψ(1)-x0))

t = range(1, 5e2, 201)

plt = plot(yaxis=:log, legend=:bottomright, title="Comparison of concatenations", xlabel="t")
plot!(plt, t, t->abs(φ(t)-x0), label="OptimalControl")
plot!(plt, t, t->abs(ψ(t)-x0), label="Classical")

Callbacks

You can use any callback from OrdinaryDiffEq.jl. For instance, we reproduce the bouncing ball example.

function condition(x, t, integrator) # Event when condition(u, t, integrator) == 0
    x[1]
end

function affect!(integrator)
    integrator.u[2] = -integrator.u[2] # the state is called u
end

cb = ContinuousCallback(condition, affect!)

g = 9.81
V(x) = [x[2], -g]
f = Flow(V; callback = cb)

t0 = 0
tf = 15
x0 = [50, 0]
sol = f((t0, tf), x0)
plot(sol)