GitHub

Double integrator: time minimisation

The problem consists in minimising the final time $t_f$ for the double integrator system

\[ \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t), \quad u(t) \in [-1,1],\]

and the limit conditions

\[ x(0) = (1,2), \quad x(t_f) = (0,0).\]

This problem can be interpretated as a simple model for a wagon with constant mass moving along a line without fricton.

First, we need to import the OptimalControl.jl package to define the optimal control problem and NLPModelsIpopt.jl to solve it. We also need to import the Plots.jl package to plot the solution.

using OptimalControl
using NLPModelsIpopt
using Plots

Optimal control problem

Let us define the problem:

ocp = @def begin

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

    -1 ≤ u(t) ≤ 1

    q(0)  == -1
    v(0)  == 0
    q(tf) == 0
    v(tf) == 0

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

    tf → min

end
Nota bene

For a comprehensive introduction to the syntax used above to define the optimal control problem, see this abstract syntax tutorial. In particular, non-Unicode alternatives are available for derivatives, integrals, etc.

Solve and plot

Direct method

Let us solve it with a direct method (we set the number of time steps to 200):

sol = solve(ocp; grid_size=200)
▫ This is OptimalControl version v1.1.1 running with: direct, adnlp, ipopt.

▫ The optimal control problem is solved with CTDirect version v0.16.4.

   ┌─ The NLP is modelled with ADNLPModels and solved with NLPModelsIpopt.
   │
   ├─ Number of time steps⋅: 200
   └─ Discretisation scheme: trapeze

▫ This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.1.

Number of nonzeros in equality constraint Jacobian...:     2004
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      402

Total number of variables............................:      604
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      201
                     variables with only upper bounds:        0
Total number of equality constraints.................:      404
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 1.10e+00 2.51e-04   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.2780101e-01 1.09e+00 6.69e+01  -6.0 1.06e+01    -  7.72e-02 1.28e-02h  4
   2  7.5366951e-01 1.03e+00 1.70e+02   0.4 9.79e+00    -  2.37e-01 5.37e-02h  3
   3  1.4473632e+00 8.70e-01 9.88e+01   0.2 5.71e+00    -  1.00e+00 1.53e-01H  1
   4  3.3885017e+00 4.33e-03 1.26e+02  -0.6 1.94e+00    -  1.00e+00 1.00e+00h  1
   5  2.2206258e+00 1.28e-03 7.75e+00  -1.5 1.17e+00   0.0 1.00e+00 1.00e+00h  1
   6  2.3300741e+00 2.17e-05 1.11e+00  -3.2 1.09e-01    -  1.00e+00 1.00e+00h  1
   7  2.0599681e+00 9.30e-04 2.59e-01  -3.6 6.91e-01    -  1.00e+00 1.00e+00f  1
   8  2.0440329e+00 2.77e-05 1.22e-02  -4.5 3.49e-01    -  1.00e+00 1.00e+00h  1
   9  2.0027584e+00 1.30e-04 8.64e-04  -5.2 8.17e-01    -  1.00e+00 9.83e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.0006050e+00 1.15e-06 4.22e-06  -5.6 2.33e-01    -  9.80e-01 1.00e+00h  1
  11  2.0000587e+00 9.27e-08 2.16e-07  -7.6 4.22e-02    -  9.87e-01 9.93e-01h  1
  12  2.0000500e+00 5.17e-11 6.45e-10 -11.0 2.26e-03    -  9.98e-01 1.00e+00h  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:   2.0000499958370912e+00    2.0000499958370912e+00
Dual infeasibility......:   6.4540800308321433e-10    6.4540800308321433e-10
Constraint violation....:   5.1717741200718592e-11    5.1717741200718592e-11
Variable bound violation:   9.9715664614308253e-09    9.9715664614308253e-09
Complementarity.........:   2.5370603077805507e-10    2.5370603077805507e-10
Overall NLP error.......:   6.4540800308321433e-10    6.4540800308321433e-10


Number of objective function evaluations             = 21
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 21
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 12
Total seconds in IPOPT                               = 1.942

EXIT: Optimal Solution Found.

and plot the solution:

plt = plot(sol; label="Direct", size=(800, 600))
Example block output
Nota bene

The solve function has options, see the solve tutorial. You can customise the plot, see the plot tutorial.

Indirect method

We now turn to the indirect method, which relies on Pontryagin’s Maximum Principle. The pseudo-Hamiltonian is given by

\[H(x, p, u) = p_1 v + p_2 u - 1,\]

where $p = (p_1, p_2)$ is the costate vector. The optimal control is of bang–bang type:

\[u(t) = \mathrm{sign}(p_2(t)),\]

with one switch from $u=+1$ to $u=-1$ at one single time denoted $t_1$. Let us implement this approach. First, we import the necessary packages:

using OrdinaryDiffEq
using NonlinearSolve

Define the bang–bang control and Hamiltonian flow:

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

# bang–bang control
u_max = +1
u_min = -1

# Hamiltonian flow
f_max = Flow(ocp, (x, p, tf) -> u_max)
f_min = Flow(ocp, (x, p, tf) -> u_min)

The shooting function enforces the conditions:

t0 = 0
x0 = [-1, 0]
xf = [ 0, 0]
function shoot!(s, p0, t1, tf)
    x_t0, p_t0 = x0, p0
    x_t1, p_t1 = f_max(t0, x_t0, p_t0, t1)
    x_tf, p_tf = f_min(t1, x_t1, p_t1, tf)
    s[1:2] = x_tf - xf                          # target conditions
    s[3] = p_t1[2]                              # switching condition
    s[4] = H(x_tf, p_tf, -1)                    # free final time
end

We are now ready to solve the shooting equations:

# in-place shooting function
nle!(s, ξ, λ) = shoot!(s, ξ[1:2], ξ[3], ξ[4])

# initial guess: costate and final time
ξ_guess = [0.1, 0.1, 0.5, 1]

# NLE problem
prob = NonlinearProblem(nle!, ξ_guess)

# resolution of the shooting equations
sol = solve(prob; show_trace=Val(true))
p0, t1, tf = sol.u[1:2], sol.u[3], sol.u[4]

# print the solution
println("\np0 = ", p0, "\nt1 = ", t1, "\ntf = ", tf)

Algorithm: NewtonRaphson(
    descent = NewtonDescent(),
    autodiff = AutoForwardDiff(),
    vjp_autodiff = AutoReverseDiff(
        compile = false
    ),
    jvp_autodiff = AutoForwardDiff(),
    concrete_jac = Val{false}()
)

----    	-------------       	-----------
Iter    	f(u) inf-norm       	Step 2-norm
----    	-------------       	-----------
0       	1.00000000e+00      	0.00000000e+00
1       	2.62500000e+00      	2.58553669e+00
2       	3.22650000e-01      	8.75907529e-01
3       	6.42532991e-03      	1.41327491e-01
4       	1.77171853e-06      	2.98492940e-03
5       	7.83817455e-14      	8.45802733e-07
Final   	7.83817455e-14
----------------------

p0 = [1.000000000000038, 1.0000000000000002]
t1 = 1.0000000000000016
tf = 2.000000000000003

Finally, we reconstruct and plot the solution obtained by the indirect method:

# concatenation of the flows
φ = f_max * (t1, f_min)

# compute the solution: state, costate, control...
flow_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 200))

# plot the solution on the previous plot
plot!(plt, flow_sol; label="Indirect", color=2, linestyle=:dash)
Example block output
Note