search

control-toolbox / CTDirect.jl

Direct transcription of an optimal control problem and resolution
http://control-toolbox.org/CTDirect.jl/
MIT License
9 stars 6 forks source link

Added implicit midpoint discretization #228

Closed PierreMartinon closed 2 months ago

PierreMartinon commented 2 months ago

@ocots @jbcaillau @joseph-gergaud Hi everyone, We now have the implicit midpoint method in addition to the default explicit trapeze. The choice is made by setting the discretization="midpoint" argument in the direct_solve call.

direct_transcription takes the argument as well, I need to add some tests on this.

Default trapeze method

direct_solve(ocp)
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:     3004
Number of nonzeros in inequality constraint Jacobian.:        2
Number of nonzeros in Lagrangian Hessian.............:     1004

Total number of variables............................:      755
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      253
                     variables with only upper bounds:        0
Total number of equality constraints.................:      504
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   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.0999999e-01 9.00e-01 1.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.0009999e-01 8.96e-01 1.95e+00  -6.0 4.80e+00    -  2.06e-03 4.04e-03h  1
   2 -1.2380702e-01 8.92e-01 1.25e+03   0.8 5.12e+01    -  1.97e-02 5.32e-03f  3
   3 -1.9816285e-01 8.85e-01 2.50e+03   1.2 3.85e+01    -  1.28e-02 7.50e-03f  3
   4 -2.1932237e-01 8.81e-01 2.49e+03   0.9 5.33e+01    -  2.38e-02 4.08e-03h  3
   5 -2.5709205e-01 8.75e-01 2.48e+03   0.9 4.29e+01    -  1.01e-01 7.64e-03h  2
   6 -1.7933148e-01 8.56e-01 2.37e+03   0.3 8.73e+00    -  1.38e-01 2.16e-02h  2
   7 -2.7929078e-01 7.77e-01 1.48e+03   0.2 6.33e+00    -  3.59e-01 9.24e-02h  2
   8 -4.9360839e-01 2.77e-01 1.10e+03  -0.4 2.82e+00    -  2.88e-01 6.43e-01h  1
   9 -7.3403809e+00 4.52e-02 2.30e+02   0.7 8.32e+00    -  9.36e-01 8.37e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -7.3882300e+00 4.49e-07 8.68e-01  -1.6 4.78e-02    -  9.89e-01 1.00e+00h  1
  11 -7.4582323e+00 5.60e-06 5.81e-03  -2.5 7.00e-02    -  9.94e-01 1.00e+00f  1
  12 -7.7592754e+00 8.50e-04 3.64e-02  -3.2 7.14e-01    -  1.00e+00 1.00e+00f  1
  13 -7.9376306e+00 5.75e-04 1.68e-02  -3.7 8.22e-01    -  9.86e-01 1.00e+00h  1
  14 -7.9860241e+00 8.97e-05 1.74e-03  -4.4 4.66e-01    -  9.92e-01 1.00e+00h  1
  15 -7.9973201e+00 2.38e-05 8.85e-05  -5.0 5.44e-01    -  9.82e-01 1.00e+00h  1
  16 -7.9997665e+00 2.18e-06 4.94e-06  -6.2 3.08e-01    -  1.00e+00 9.88e-01h  1
  17 -7.9999664e+00 3.20e-08 1.14e-07 -11.0 4.44e-02    -  9.56e-01 9.92e-01h  1
  18 -7.9999681e+00 5.36e-12 3.66e-12 -11.0 1.37e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 18

                                   (scaled)                 (unscaled)
Objective...............:  -7.9999681066635722e+00   -7.9999681066635722e+00
Dual infeasibility......:   3.6628478028433165e-12    3.6628478028433165e-12
Constraint violation....:   5.3607118744025684e-12    5.3607118744025684e-12
Variable bound violation:   9.9989831525704176e-08    9.9989831525704176e-08
Complementarity.........:   2.8843111065815013e-11    2.8843111065815013e-11
Overall NLP error.......:   2.8843111065815013e-11    2.8843111065815013e-11

Number of objective function evaluations             = 36
Number of objective gradient evaluations             = 19
Number of equality constraint evaluations            = 36
Number of inequality constraint evaluations          = 36
Number of equality constraint Jacobian evaluations   = 19
Number of inequality constraint Jacobian evaluations = 19
Number of Lagrangian Hessian evaluations             = 18
Total seconds in IPOPT                               = 0.145

EXIT: Optimal Solution Found.

Implicit midpoint method (note the larger NLP size due to the k_i variables)

direct_solve(ocp, discretization = "midpoint")
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:     3754
Number of nonzeros in inequality constraint Jacobian.:        2
Number of nonzeros in Lagrangian Hessian.............:     1000

Total number of variables............................:     1254
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      252
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1004
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   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.0999999e-01 9.00e-01 5.00e-01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.0009999e-01 8.96e-01 9.54e-01  -6.0 4.80e+00    -  2.06e-03 4.04e-03h  1
   2 -1.2386037e-01 8.92e-01 1.26e+03   0.8 5.10e+01    -  1.97e-02 5.33e-03f  3
   3 -1.9794118e-01 8.85e-01 2.49e+03   1.2 3.86e+01    -  1.28e-02 7.49e-03f  3
   4 -2.1896928e-01 8.81e-01 2.49e+03   0.9 5.34e+01    -  2.37e-02 4.07e-03h  3
   5 -2.5653279e-01 8.75e-01 2.47e+03   0.9 4.30e+01    -  1.02e-01 7.63e-03h  2
   6 -1.7904905e-01 8.56e-01 2.37e+03   0.3 8.71e+00    -  1.38e-01 2.13e-02h  2
   7 -2.2959007e-01 8.16e-01 2.09e+03   0.2 6.37e+00    -  3.59e-01 4.62e-02h  3
   8 -2.2631164e+00 5.32e-02 8.62e+03  -0.2 4.28e+00    -  1.77e-01 9.35e-01H  1
   9 -4.2444467e+00 1.95e-05 9.57e+02  -0.0 1.98e+00    -  6.91e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -4.2498587e+00 1.67e-09 8.11e-01  -1.6 5.41e-03   2.0 1.00e+00 1.00e+00h  1
  11 -7.5466761e+00 1.77e-09 7.96e-01  -7.6 2.01e+02    -  2.00e-02 1.64e-02f  1
  12 -7.5799492e+00 3.17e-09 7.83e-01  -2.6 1.80e+00    -  1.00e+00 1.85e-02f  1
  13 -7.8912006e+00 8.99e-04 1.69e-01  -3.1 7.21e-01    -  1.00e+00 1.00e+00f  1
  14 -7.9314294e+00 4.10e-05 1.80e-02  -4.0 2.59e-01    -  9.99e-01 1.00e+00h  1
  15 -7.9921542e+00 1.39e-04 1.41e-03  -4.6 5.75e-01    -  1.00e+00 1.00e+00h  1
  16 -7.9974291e+00 2.10e-05 2.02e-04  -5.2 6.97e-01    -  9.52e-01 8.75e-01h  1
  17 -7.9997761e+00 6.38e-06 3.79e-06  -6.1 6.82e-01    -  9.52e-01 1.00e+00h  1
  18 -7.9999999e+00 3.69e-08 8.46e-08 -11.0 4.14e-02    -  9.80e-01 1.00e+00h  1
  19 -8.0000001e+00 3.79e-13 2.07e-14 -11.0 5.53e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 19

                                   (scaled)                 (unscaled)
Objective...............:  -8.0000001074856755e+00   -8.0000001074856755e+00
Dual infeasibility......:   2.0650148258027912e-14    2.0650148258027912e-14
Constraint violation....:   3.7936320751441599e-13    3.7936320751441599e-13
Variable bound violation:   9.9990000279603919e-08    9.9990000279603919e-08
Complementarity.........:   1.4993415768214275e-11    1.4993415768214275e-11
Overall NLP error.......:   1.4993415768214275e-11    1.4993415768214275e-11

Number of objective function evaluations             = 39
Number of objective gradient evaluations             = 20
Number of equality constraint evaluations            = 39
Number of inequality constraint evaluations          = 39
Number of equality constraint Jacobian evaluations   = 20
Number of inequality constraint Jacobian evaluations = 20
Number of Lagrangian Hessian evaluations             = 19
Total seconds in IPOPT                               = 0.484

Now that all the plumbing is in place, it should be quite easy to add new discretization schemes. I'll add a 2-stage gauss method to check, then back to allocation / performance aspects.

github-actions[bot] commented 2 months ago
Package name latest stable
OptimalControl.jl compat: v0.12.0 compat: v0.12.0
ocots commented 2 months ago

Could you give the possibility to provide a String or a Symbol? 

discretization = "midpoint"
# or 
discretization = :midpoint

You can simply add:

function direct_solve(
    ocp::OptimalControlModel,
    description::Symbol...;
    init = CTBase.__ocp_init(),
    grid_size::Int = CTDirect.__grid_size(),
    time_grid = CTDirect.__time_grid(),
    discretization::Union{String, Symbol} = __discretization(),
    kwargs...,
)
discretization = string(discretization)
...
end

# Idem for direct_transcription

[!NOTE] It would be nice to have this for all options given by a Symbol.

PierreMartinon commented 2 months ago

Could you give the possibility to provide a String or a Symbol?

discretization = "midpoint"
# or 
discretization = :midpoint

You can simply add:

function direct_solve(
    ocp::OptimalControlModel,
    description::Symbol...;
    init = CTBase.__ocp_init(),
    grid_size::Int = CTDirect.__grid_size(),
    time_grid = CTDirect.__time_grid(),
    discretization::Union{String, Symbol} = __discretization(),
    kwargs...,
)
discretization = string(discretization)
...
end

# Idem for direct_transcription

Note

It would be nice to have this for all options given by a Symbol.

Sure. By the way, yes, we should adopt the same rule for all arguments. Maybe always allow for String/Symbol ? If we want to keep things compact we could write a basic equality function in CTBase that would work for both types. Also define a type for the Union, something like TextOption ?

On the other hand, just adding the string conversion as you wrote is not that hard.

PierreMartinon commented 2 months ago

Note: the aforementioned plumbing caused some additional allocations and slower execution than the hardcoded trapeze. I'll investigate this before merging.

ocots commented 2 months ago

The conversion for all the options in kwargs...:

function replace_symbols_by_strings(; kwargs...)
    return replace(p -> p.second isa Symbol ? p.first => string(p.second) : p, kwargs)
end

and the extension of equality:

import Base: :(==)
Base.:(==)(a::String, b::Symbol) = a == string(b)
Base.:(==)(a::Symbol, b::String) = Base.:(==)(b, a)
jbcaillau commented 2 months ago

@PierreMartinon nice! yes, symbols are more or less the standard for julia options

ocots commented 2 months ago

@PierreMartinon You can create if you want a file src/ctbase.jl where you put everything you think that should be placed in CTBase.jl.

ocots commented 2 months ago

@PierreMartinon Je vois que tu es devenu fan des "Tag" :-)

PierreMartinon commented 2 months ago

@PierreMartinon Je vois que tu es devenu fan des "Tag" :-)

Actually I went back a bit from the separate tag, now the whole docp is parametrized by the discretization ie DOCP{Trapeze} or DOCP{Midpoint}. It seemed to make things a bit more explicit for julia and removed the extra allocations I first had when introducing multiple dispatch (still not completely clear why this happened).

Merge and release probably next week, I'd like to add the Gauss Legendre 4th order scheme as well. At this point we should be close to the generic implicit RK as well.

jbcaillau commented 2 months ago

@PierreMartinon seen that you have added quite a lot of things 👍🏽!

ocots commented 2 months ago

Such a big PR!