|
| 1 | +# Discretisation methods |
| 2 | + |
| 3 | +## Discretisation formulas |
| 4 | +When calling `solve`, the option `disc_method=...` can be used to set the discretisation scheme. |
| 5 | +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. |
| 6 | +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. |
| 7 | + |
| 8 | +As an example we will use the [Goddard problem](@ref tutorial-goddard) |
| 9 | +```@example main |
| 10 | +using OptimalControl # to define the optimal control problem and more |
| 11 | +using NLPModelsIpopt # to solve the problem via a direct method |
| 12 | +using Plots # to plot the solution |
| 13 | +
|
| 14 | +t0 = 0 # initial time |
| 15 | +r0 = 1 # initial altitude |
| 16 | +v0 = 0 # initial speed |
| 17 | +m0 = 1 # initial mass |
| 18 | +vmax = 0.1 # maximal authorized speed |
| 19 | +mf = 0.6 # final mass to target |
| 20 | +
|
| 21 | +ocp = @def begin # definition of the optimal control problem |
| 22 | +
|
| 23 | + tf ∈ R, variable |
| 24 | + t ∈ [t0, tf], time |
| 25 | + x = (r, v, m) ∈ R³, state |
| 26 | + u ∈ R, control |
| 27 | +
|
| 28 | + x(t0) == [ r0, v0, m0 ] |
| 29 | + m(tf) == mf, (1) |
| 30 | + 0 ≤ u(t) ≤ 1 |
| 31 | + r(t) ≥ r0 |
| 32 | + 0 ≤ v(t) ≤ vmax |
| 33 | +
|
| 34 | + ẋ(t) == F0(x(t)) + u(t) * F1(x(t)) |
| 35 | +
|
| 36 | + r(tf) → max |
| 37 | +
|
| 38 | +end; |
| 39 | +
|
| 40 | +# Dynamics |
| 41 | +const Cd = 310 |
| 42 | +const Tmax = 3.5 |
| 43 | +const β = 500 |
| 44 | +const b = 2 |
| 45 | +
|
| 46 | +F0(x) = begin |
| 47 | + r, v, m = x |
| 48 | + D = Cd * v^2 * exp(-β*(r - 1)) # Drag force |
| 49 | + return [ v, -D/m - 1/r^2, 0 ] |
| 50 | +end |
| 51 | +
|
| 52 | +F1(x) = begin |
| 53 | + r, v, m = x |
| 54 | + return [ 0, Tmax/m, -b*Tmax ] |
| 55 | +end |
| 56 | +nothing # hide |
| 57 | +``` |
| 58 | +Now let us compare different discretisations |
| 59 | +```@example main |
| 60 | +sol_trapeze = solve(ocp; tol=1e-8) |
| 61 | +plot(sol_trapeze) |
| 62 | +
|
| 63 | +sol_midpoint = solve(ocp, disc_method=:midpoint; tol=1e-8) |
| 64 | +plot!(sol_midpoint) |
| 65 | +
|
| 66 | +sol_gl2 = solve(ocp, disc_method=:gauss_legendre_2; tol=1e-8) |
| 67 | +plot!(sol_gl2) |
| 68 | +
|
| 69 | +sol_gl3 = solve(ocp, disc_method=:gauss_legendre_3; tol=1e-8) |
| 70 | +plot!(sol_gl3) |
| 71 | +``` |
| 72 | + |
| 73 | +## Large problems and AD backend |
| 74 | +For some large problems, you may notice that solving spends a long time before the iterations actually begin. |
| 75 | +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. |
| 76 | +A possible alternative is to set the option `adnlp_backend=:manual`, which will use more basic sparsity patterns. |
| 77 | +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. |
| 78 | + |
| 79 | +```@example main |
| 80 | +solve(ocp, disc_method=:gauss_legendre_3, grid_size=1000, adnlp_backend=:manual) |
| 81 | +nothing # hide |
| 82 | +``` |
| 83 | + |
| 84 | +## Explicit time grid |
| 85 | +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. |
| 86 | +In the case of a free initial and/or final time, provide a normalised grid between 0 and 1. |
| 87 | +Note that `time_grid` will override `grid_size` if both are present. |
| 88 | + |
| 89 | +```@example main |
| 90 | +sol = solve(ocp, time_grid=[0, 0.1, 0.5, 0.9, 1], display=false) |
| 91 | +println(time_grid(sol)) |
| 92 | +``` |
0 commit comments