Discretisation methods
In this tutorial, we will explore the different discretisation methods available in the OptimalControl.jl package. These methods are used to convert a continuous-time optimal control problem (OCP) into a discrete-time nonlinear programming (NLP) problem, which can then be solved using numerical optimization techniques.
Let us import the necessary packages and define the optimal control problem (Goddard problem) we will use as an example throughout this tutorial.
using BenchmarkTools # for benchmarking
using DataFrames # to store the results
using OptimalControl # to define the optimal control problem and more
using NLPModels # to retrieve information about the NLP problem
using NLPModelsIpopt # to solve the problem via a direct method
using Plots # to plot the solution
const t0 = 0 # initial time
const r0 = 1 # initial altitude
const v0 = 0 # initial speed
const m0 = 1 # initial mass
const vmax = 0.1 # maximal authorized speed
const mf = 0.6 # final mass to target
ocp = @def begin # definition of the optimal control problem
tf ∈ R, variable
t ∈ [t0, tf], time
x = (r, v, m) ∈ R³, state
u ∈ R, control
x(t0) == [r0, v0, m0]
m(tf) == mf
0 ≤ u(t) ≤ 1
r(t) ≥ r0
0 ≤ v(t) ≤ vmax
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
r(tf) → max
end;
# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 2
F0(x) = begin
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
return [ v, -D/m - 1/r^2, 0 ]
end
F1(x) = begin
r, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
endDiscretisation schemes
When calling solve, the option disc_method=... can be used to specify the discretization scheme. In addition to the default implicit :trapeze method (also known as Crank–Nicolson), other available options include the implicit :midpoint method, and Gauss–Legendre collocation schemes with 2 and 3 stages: :gauss_legendre_2 and :gauss_legendre_3, of order 4 and 6 respectively. Note that higher-order methods typically result in larger NLP problems for the same number of time steps, and their accuracy also depends on the smoothness of the problem.
Let us first solve the problem with the default :trapeze method and display the solution.
sol = solve(ocp; disc_method=:trapeze, display=false)
plot(sol; size=(800, 800))
Let us now compare different discretization schemes to evaluate their accuracy and performance.
# Solve the problem with different discretization methods
solutions = []
data = DataFrame(
Scheme=Symbol[],
Time=Float64[],
Objective=Float64[],
Iterations=Int[]
)
schemes = [
:euler,
:euler_implicit,
:trapeze,
:midpoint,
:gauss_legendre_2,
:gauss_legendre_3
]
for scheme in schemes
bt = @btimed solve($ocp; disc_method=$scheme, tol=1e-8, display=false)
sol = bt.value
push!(solutions, (scheme, sol))
push!(data, (
Scheme=scheme,
Time=bt.time,
Objective=objective(sol),
Iterations=iterations(sol)
))
end
println(data)┌ Warning: Assignment to `sol` in soft scope is ambiguous because a global variable by the same name exists: `sol` will be treated as a new local. Disambiguate by using `local sol` to suppress this warning or `global sol` to assign to the existing global variable.
└ @ tutorial-discretisation.md:93
6×4 DataFrame
Row │ Scheme Time Objective Iterations
│ Symbol Float64 Float64 Int64
─────┼───────────────────────────────────────────────────
1 │ euler 0.133998 1.01242 30
2 │ euler_implicit 0.134035 1.01273 29
3 │ trapeze 0.149832 1.01258 30
4 │ midpoint 0.209906 1.01258 31
5 │ gauss_legendre_2 0.596406 1.01258 21
6 │ gauss_legendre_3 1.28157 1.01258 21# Plot the results
x_style = (legend=:none,)
p_style = (legend=:none,)
styles = (state_style=x_style, costate_style=p_style)
scheme, sol = solutions[1]
plt = plot(sol; label=string(scheme), styles...)
for (scheme, sol) in solutions[2:end]
plt = plot!(sol; label=string(scheme), styles...)
end
plot(plt; size=(800, 800))
Large problems and AD backend
For some large problems, you may notice that the solver takes a long time before the iterations actually begin. This is due to the computation of sparse derivatives — specifically, the Jacobian of the constraints and the Hessian of the Lagrangian — which can be quite costly. One possible alternative is to set the option adnlp_backend=:manual, which uses simpler sparsity patterns. The resulting matrices are faster to compute but are also less sparse, so this represents a trade-off between automatic differentiation (AD) preparation time and the efficiency of the optimization itself.
solve(ocp; disc_method=:gauss_legendre_3, grid_size=1000, adnlp_backend=:manual)This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.
Number of nonzeros in equality constraint Jacobian...: 174028
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 104022
Total number of variables............................: 13004
variables with only lower bounds: 1001
variables with lower and upper bounds: 2001
variables with only upper bounds: 0
Total number of equality constraints.................: 12004
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.0100000e+00 2.22e+00 2.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0359714e+00 2.59e+00 2.26e+04 -1.0 6.29e+00 - 1.26e-01 9.81e-01f 1
2 -1.0421591e+00 2.58e+00 2.25e+04 0.5 1.61e+02 - 1.00e+00 5.62e-03f 1
3 -1.0581309e+00 2.37e+00 1.98e+05 1.6 2.29e+01 - 1.00e+00 7.89e-02f 1
4 -1.0297827e+00 1.37e+00 5.30e+05 2.3 1.01e+01 - 9.95e-01 6.06e-01h 1
5 -1.0168610e+00 7.22e-01 1.25e+06 2.3 4.82e+00 - 5.95e-01 9.84e-01h 1
6 -1.0101725e+00 8.73e-02 2.61e+07 2.3 2.06e+00 - 1.00e+00 8.67e-01h 1
7 -1.0088127e+00 5.28e-03 7.50e+05 1.6 6.07e-01 - 1.00e+00 9.97e-01h 1
8 -1.0088898e+00 7.84e-05 1.14e+03 0.2 9.19e-02 - 1.00e+00 1.00e+00f 1
9 -1.0088804e+00 2.32e-06 1.10e+01 -1.5 2.05e-02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -1.0088776e+00 2.71e-07 9.58e-02 -3.6 7.53e-03 - 1.00e+00 1.00e+00h 1
11 -1.0089150e+00 1.08e-05 1.50e+01 -5.5 1.88e-02 - 9.99e-01 1.00e+00h 1
12 -1.0109999e+00 4.39e-02 1.70e-03 -6.1 8.30e-01 - 1.00e+00 1.00e+00f 1
13 -1.0121802e+00 3.92e-02 1.06e+01 -6.6 9.51e-01 - 1.00e+00 7.90e-01h 1
14 -1.0124904e+00 9.56e-03 4.00e-05 -7.0 7.89e-01 - 1.00e+00 1.00e+00h 1
15 -1.0125500e+00 3.91e-03 3.80e-01 -7.5 6.77e-01 - 1.00e+00 9.43e-01h 1
16 -1.0125686e+00 2.11e-03 4.64e-06 -8.0 6.09e-01 - 1.00e+00 1.00e+00h 1
17 -1.0125744e+00 9.51e-04 8.30e-03 -8.6 5.47e-01 - 1.00e+00 9.89e-01h 1
18 -1.0125756e+00 4.32e-04 2.82e-02 -9.2 6.93e-01 - 1.00e+00 8.38e-01h 1
19 -1.0125759e+00 1.83e-04 6.00e-05 -9.3 6.45e-01 - 9.94e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 -1.0125762e+00 1.07e-04 3.43e-08 -10.0 3.84e-01 - 1.00e+00 1.00e+00h 1
21 -1.0125763e+00 9.04e-05 1.61e-08 -11.0 6.51e-01 - 1.00e+00 1.00e+00h 1
22 -1.0125763e+00 2.05e-05 2.14e-09 -11.0 2.98e-01 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 22
(scaled) (unscaled)
Objective...............: -1.0125763121308899e+00 -1.0125763121308899e+00
Dual infeasibility......: 2.1389119035387901e-09 2.1389119035387901e-09
Constraint violation....: 4.8819803755151270e-10 2.0494566410900195e-05
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 1.6640648931629154e-11 1.6640648931629154e-11
Overall NLP error.......: 2.1389119035387901e-09 2.0494566410900195e-05
Number of objective function evaluations = 23
Number of objective gradient evaluations = 23
Number of equality constraint evaluations = 23
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 23
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 22
Total seconds in IPOPT = 7.847
EXIT: Optimal Solution Found.Let us now compare the performance of the two backends. We will use the @btimed macro from the BenchmarkTools package to measure the time taken for both the preparation of the NLP problem and the execution of the solver. Besides, we will also collect the number of non-zero elements in the Jacobian and Hessian matrices, which can be useful to understand the sparsity of the problem, thanks to the functions get_nnzo, get_nnzj, and get_nnzj from the NLPModels package. The problem is first discretized with the direct_transcription method and then solved with the ipopt solver, see the tutorial on direct transcription for more details.
# DataFrame to store the results
data = DataFrame(
Backend=Symbol[],
NNZO=Int[],
NNZJ=Int[],
NNZH=Int[],
PrepTime=Float64[],
ExecTime=Float64[],
Objective=Float64[],
Iterations=Int[]
)
# The different AD backends to test
backends = [:optimized, :manual]
# Loop over the backends
for adnlp_backend in backends
# Discretize the problem with a large grid size and Gauss-Legendre method
bt = @btimed direct_transcription($ocp;
disc_method=:gauss_legendre_3,
grid_size=10,
adnlp_backend=$adnlp_backend,
)
docp, nlp = bt.value
prepa_time = bt.time
# Get the number of non-zero elements
nnzo = get_nnzo(nlp) # Gradient of the Objective
nnzj = get_nnzj(nlp) # Jacobian of the constraints
nnzh = get_nnzh(nlp) # Hessian of the Lagrangian
# Solve the problem
bt = @btimed ipopt($nlp;
print_level=0,
mu_strategy="adaptive",
tol=1e-8,
sb="yes",
)
nlp_sol = bt.value
exec_time = bt.time
# Store the results in the DataFrame
push!(data,
(
Backend=adnlp_backend,
NNZO=nnzo,
NNZJ=nnzj,
NNZH=nnzh,
PrepTime=prepa_time,
ExecTime=exec_time,
Objective=-nlp_sol.objective,
Iterations=nlp_sol.iter
)
)
end
println(data)2×8 DataFrame
Row │ Backend NNZO NNZJ NNZH PrepTime ExecTime Objective Iterations
│ Symbol Int64 Int64 Int64 Float64 Float64 Float64 Int64
─────┼──────────────────────────────────────────────────────────────────────────────
1 │ optimized 134 874 951 0.00612104 0.0323602 1.01255 16
2 │ manual 134 1768 1062 0.00190795 0.0376257 1.01255 16Explicit time grid
The option time_grid=... allows you to provide the full time grid vector t0, t1, ..., tf, which is especially useful if a non-uniform grid is desired. In the case of a free initial and/or final time, you should provide a normalized grid ranging from 0 to 1. Note that time_grid overrides grid_size if both options are specified.
sol = solve(ocp; time_grid=[0, 0.1, 0.5, 0.9, 1], display=false)
println(time_grid(sol))[0.0, 0.016495556123488363, 0.08247778061744182, 0.1484600051113953, 0.16495556123488364]