crude example of basic grid refinement

[PierreMartinon]

@joseph-gergaud To illustrate the possibility of a grid refinement algorithm, we do a discrete continuation on the discretization steps. We use limited iterations until we reach the desired steps, then do a final solve. Apparently this can be faster than than directlty solving the problem.

NB. In this example we only used uniform grids so we could just have passe grid_size instead of time_grid, but the point of grid refinement would to use non-uniform grids.

julia> include("test/test_grid.jl")
Test: grid options

Basic solve
steps  250, objective  1.012576, iterations   46
 15.297255 seconds (41.06 M allocations: 5.121 GiB, 1.38% gc time, 1.55% compilation time)

Step continuation
steps   10, objective  1.000258, iterations set to    5
steps   20, objective  1.004950, iterations set to    5
steps   40, objective  1.009693, iterations set to    5
steps   80, objective  1.012172, iterations set to    5
steps  160, objective  1.012389, iterations set to    5
steps  250, objective  1.011619, iterations set to    5
steps  250, objective  1.012576, iterations   14
  9.759167 seconds (21.74 M allocations: 2.728 GiB, 7.79% gc time)

[ocots]

I am about to add an application. I do a refinement to increase the speed of convergence. See here.

[PierreMartinon]

And a first try at optimizing the time steps themselves as variables (as in jump). Code is ugly, but hey, it gives something. Test case is basic double integrator min tf, the one with a bang-bang solution. With 10 steps to optimize, we obtain convergence, and the control looks ok, with the optimized time steps concentrated at t0, tf and the switching at mid-interval.

EXIT: Optimal Solution Found.
Optimized time steps [0.0, 0.009999994000050694, 0.9302150943487946, 0.9840614282533159, 0.9947620511971743, 1.005262939119436, 1.0248619765421794, 1.0508788049261277, 1.1209608643384075, 1.9900249959238299, 2.000024989923881]
And tf: 2.000024989923881

For the fearless, you can check the implementation in test_grid.jl, using variables for the time steps Dt, a fixed integer time grid of 'indices' for the Dt, and a manual rescaling of the dynamics (with trapeze rule in mind). The solution returned is interpolated against the fake time grid though, so the functions are incorrect and have to be rescaled manually. Just changing the time grid in the solution afterwards is not enough since the interpolations are already done.

Obviously this feature would require a specific option passed to solve, with an adjusted internal implementation. Might not be so hard though, since we can basically do it by cheating the ocp formulation :D

    N_vars = 10
    function f(t,x,u,v)
        if t == 0
            dt = v[1]
        elseif t == N_vars
            dt = v[end]
        else
            dt = 0.5 * (v[Int(t)] + v[Int(t)+1])
        end
        return [x[2], u] * dt
    end

    ocp = Model(variable=true, autonomous=false)
    state!(ocp, 2)
    control!(ocp, 1)
    variable!(ocp, N_vars)
    time!(ocp, t0=0, tf=N_vars)
    constraint!(ocp, :initial, lb=[0,0], ub=[0,0])
    constraint!(ocp, :final, lb=[1,0], ub=[1,0])
    constraint!(ocp, :control, lb=-1, ub=1)
    constraint!(ocp, :variable, lb=0.01*ones(N_vars), ub=10*ones(N_vars))
    dynamics!(ocp, f)
    objective!(ocp, :mayer, (x0, xf, v) -> sum(v))

    sol = solve(ocp, grid_size=N_vars)

[PierreMartinon]

Tests with goddard (free tf, max rf, v constraint) are not great, and the time grid found is very sensitive to the bounds on the Dt_i. Here are two sample controls.

Maybe the optimization finds 'bad' values for the time steps so it can cheat the ODE and get a better objective...
Also, there may be a derivatives problem for the external function dt.

This variable grid feature needs a proper implementation if we want to do more serious testing.

[ocots]

I agree.

[PierreMartinon]

For the syntax, maybe we can reuse the existing arguments with something like
solve(ocp, grid_size=N, time_grid=:optimized)