Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added implicit midpoint discretization #228

Merged
merged 42 commits into from
Sep 2, 2024
Merged

Added implicit midpoint discretization #228

merged 42 commits into from
Sep 2, 2024

Conversation

PierreMartinon
Copy link
Member

@PierreMartinon PierreMartinon commented Aug 29, 2024

@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.

@PierreMartinon PierreMartinon marked this pull request as draft August 29, 2024 10:17
@PierreMartinon PierreMartinon marked this pull request as ready for review August 29, 2024 10:22
@PierreMartinon PierreMartinon marked this pull request as draft August 29, 2024 10:31
Copy link
Contributor

github-actions bot commented Aug 29, 2024

Package name latest stable
OptimalControl.jl compat: v0.12.0 compat: v0.12.0

@ocots
Copy link
Member

ocots commented Aug 29, 2024

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
Copy link
Member Author

PierreMartinon commented Aug 29, 2024

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
Copy link
Member Author

PierreMartinon commented Aug 29, 2024

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

@ocots
Copy link
Member

ocots commented Aug 29, 2024

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
Copy link
Member

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

@ocots
Copy link
Member

ocots commented Aug 30, 2024

@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
Copy link
Member

ocots commented Aug 30, 2024

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

@PierreMartinon
Copy link
Member Author

PierreMartinon commented Aug 30, 2024

@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.

@PierreMartinon PierreMartinon marked this pull request as ready for review September 2, 2024 15:38
@jbcaillau
Copy link
Member

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

@PierreMartinon PierreMartinon merged commit a050125 into main Sep 2, 2024
5 checks passed
@PierreMartinon PierreMartinon deleted the midpoint branch September 2, 2024 18:08
@ocots
Copy link
Member

ocots commented Sep 2, 2024

Such a big PR!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants