Discretisation methods
Discretisation formulas
When calling solve, the option disc_method=... can be used to set the discretisation scheme. In addition to the default implicit :trapeze method (aka Crank-Nicolson), other choices are available, namely implicit :midpoint and the Gauss-Legendre collocations with 2 and stages, :gauss_legendre_2 and :gauss_legendre_3, of order 4 and 6 respectively. Note that higher order methods will typically lead to larger NLP problems for the same number of time steps, and that accuracy will also depend on the smoothness of the problem.
As an example we will use the Goddard problem
using OptimalControl # to define the optimal control problem and more
using NLPModelsIpopt # to solve the problem via a direct method
using Plots # to plot the solution
t0 = 0 # initial time
r0 = 1 # initial altitude
v0 = 0 # initial speed
m0 = 1 # initial mass
vmax = 0.1 # maximal authorized speed
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, (1)
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 ]
endNow let us compare different discretisations
sol_trapeze = solve(ocp; tol=1e-8)
plot(sol_trapeze)
sol_midpoint = solve(ocp, disc_method=:midpoint; tol=1e-8)
plot!(sol_midpoint)
sol_euler = solve(ocp, disc_method=:euler; tol=1e-8)
plot!(sol_euler)
sol_euler_imp = solve(ocp, disc_method=:euler_implicit; tol=1e-8)
plot!(sol_euler_imp)
sol_gl2 = solve(ocp, disc_method=:gauss_legendre_2; tol=1e-8)
plot!(sol_gl2)
sol_gl3 = solve(ocp, disc_method=:gauss_legendre_3; tol=1e-8)
plot!(sol_gl3)
Large problems and AD backend
For some large problems, you may notice that solving spends a long time before the iterations actually begin. This is due to the computing of the sparse derivatives, namely the Jacobian of the constraints and the Hessian of the Lagrangian, that can become quite costly. A possible alternative is to set the option adnlp_backend=:manual, which will use more basic sparsity patterns. The resulting matrices are faster to compute but are also less sparse, so this is a trade-off bewteen the AD preparation and 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.7.3.
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
Number of Iterations....: 21
(scaled) (unscaled)
Objective...............: -1.0125763121283586e+00 -1.0125763121283586e+00
Dual infeasibility......: 1.6147371975487365e-08 1.6147371975487365e-08
Constraint violation....: 3.8848892991438788e-09 9.0393234664976063e-05
Variable bound violation: 2.7626403077445802e-31 2.7626403077445802e-31
Complementarity.........: 6.3724399828030864e-11 6.3724399828030864e-11
Overall NLP error.......: 1.6147371975487365e-08 9.0393234664976063e-05
Number of objective function evaluations = 22
Number of objective gradient evaluations = 22
Number of equality constraint evaluations = 22
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 22
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 21
Total seconds in IPOPT = 8.506
EXIT: Optimal Solution Found.Explicit time grid
The option time_grid=... allows to pass the complete time grid vector t0, t1, ..., tf, which is typically useful if one wants a non uniform grid. In the case of a free initial and/or final time, provide a normalised grid between 0 and 1. Note that time_grid will override grid_size if both are present.
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]