Energy minimisation

A wagon moves along a rail, and we control its acceleration with a force . Its state is  , where is the position and the velocity.

The mass is one and there is no friction, so the dynamics are

the double integrator. We transfer the wagon from rest at   , at  , to rest at  , at  , with the least energy,  . This is the problem of Your first problem: this page solves it by the direct and the indirect methods, and compares them with the exact solution.

We use OptimalControl to define the problem, NLPModelsIpopt to solve it by a direct method, OrdinaryDiffEqTsit5 and NonlinearSolve for the indirect method, and Plots to see the result.

using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

The problem

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

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

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

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

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

See Abstract syntax (@def) for every form these lines can take.

Direct method

The direct method^[1] discretises the problem in time and solves the resulting nonlinear program, here with Ipopt:

direct_sol = solve(ocp; display=false)
objective(direct_sol)

6.000096001536038

plt = plot(direct_sol; label="direct", size=(800, 600))

The exact solution

With the conventions of Notation and conventions (  ), the pseudo-Hamiltonian is

It is maximised by  . The adjoint equations,   and   , give   and   , so the control is affine,   . The final conditions,     and      , give   and  :

Indirect method

The indirect method, indirect simple shooting, computes the extremal from its initial costate. The flow of the Hamiltonian system, with the maximising control in feedback form, integrates from :

# maximising control, H(x, p, u) = p₁v + p₂u - u²/2
u_max(x, p) = p[2]

f = Flow(ocp, u_max)

The shooting function measures how far the flow from lands from the target:

S(p0) = f(t0, x0, p0, tf)[1] - xf

We solve   with a Newton method, starting from the costate of the direct solution:

nle!(s, p0, _) = (s .= S(p0))

p0_guess = costate(direct_sol)(t0)
prob = NonlinearProblem(nle!, p0_guess)
shooting_sol = NonlinearSolve.solve(prob; show_trace=Val(false))
p0_sol = shooting_sol.u

2-element Vector{Float64}:
 12.000000000000087
  6.000000000000021

The residual is zero, and the initial costate is the exact one, :

S(p0_sol)

2-element Vector{Float64}:
 3.438776657013338e-15
 7.32695851087319e-15

Comparison

The extremal integrated from p0_sol is a solution of the problem, with its cost:

indirect_sol = f((t0, tf), x0, p0_sol)
objective(direct_sol), objective(indirect_sol)

(6.000096001536038, 5.999999999999881)

For the plot, the flow is built with saveat, so that the trajectory has 101 points (see Plotting a flow):

f_plot = Flow(ocp, u_max; saveat=range(t0, tf, 101), dense=false)
plot!(plt, f_plot((t0, tf), x0, p0_sol); label="indirect", linestyle=:dash)

The indirect cost is the exact value, . The direct cost differs by the discretisation error of the default grid.

Notes

• MINPACK.jl can replace NonlinearSolve to solve the shooting equations.

• From an OCP explains the flow of a problem, and Shooting the shooting method.

• The Goddard tutorial uses the same steps on a problem with bang, singular and boundary arcs.

• State constraint adds a bound on the velocity or the position to this problem.

---

1. J. T. Betts, Practical methods for optimal control using nonlinear programming, SIAM, Philadelphia, PA, 2001. ↩︎