Export _OptimalControlSolution

[paraynaud]

Hello, I recently gave a student lab in which I used the great work done in #73 to extract an ADLNPModel:

docp = DirectTranscription(ocp, grid_size=n)
nlp = getNLP(docp)

I asked my students to solve the problem using:

• sol = solveDOCP(docp) which uses Ipopt;
• percival(nlp) a solver from JuliaSmoothOptimizers/Percival.jl;
  and to print the solution, which is straightforward once sol is obtained, i.e. plot(sol).

However, it is less direct from the solution returned by percival.
I found the function _OptimalControlSolution(...) to build a proper solution from the GenericExecutionStats that percival returns:

ges = percival(nlp)
sol = CTDirect._OptimalControlSolution(ges, docp)

which can be printed with plot(sol) later.

As the ADNLPModel is not only embedded in docp anymore, It would be quite interesting for users wanting to test their own solvers to highlight the feature that _OptimalControlSolution provides.

Currently, it needs a full GenericExecutionStats, but I believe the resulting vector should be enough, as the rest of the information could be inferred from docp.
That way, a user could implement its own solver without having to fully understand what a GenericExecutionStats is.

Thanks for your work !

[PierreMartinon]

Hi @paraynaud , great timing :D
I had just renamed it and exported as OCPSolutionFromDOCP (still on the export branch for now, merge pending). Indeed the aim is to make it as robust as possible regarding the NLP solver. I'll start to see what we can use from the DOCP instead of the GenericExecutionStats.

Don't hesitate to put your feedback here !

[paraynaud]

Thanks @PierreMartinon for adding those features.
I gave an extensive feedback to Jean Baptiste Caillau, but briefly, OptimalControl makes optimal control problem really simple to model and solve ; while the plot() method makes it easy for any user to have the graphical solutions !
Practically, my students grasped pretty quickly the functionalities of OptimalControl and cleared easily the lab's instructions.

[PierreMartinon]

Doc pending.

[jbcaillau]

Hi @paraynaud ; just gave a try to your problem: solves fine 👍🏽. Regarding the dynamics, or expressions in @def in general, you can use whatever julia code you want (provided it is defined at call time 😅). Check for instance the function $\gamma$ below (as it depends on time, the package detects that the problem is NonAutonomous).

# invest.jl

using OptimalControl

T = 1
β = 0.2
γ(t, x) = 3 * exp(-β * t) # yield
x0 = 0.1 # initial capital

@def ocp begin
  t ∈ [ 0, T ], time
  x ∈ R, state
  u ∈ R, control
  0 ≤ u(t) ≤ 1
  
  x(0) == x0
  ẋ(t) == γ(t, x(t)) * u(t) * x(t) # capital dynamics

  ∫( (1 - u(t)) * x(t) ) → max
end

sol = solve(ocp)

plot(sol)

[jbcaillau]

@paraynaud please consider adding some comments / motivation / context for your example and including it in OptimalControl.jl tutos

[PierreMartinon]

@paraynaud I did a more 'raw' version of the export, is this what you had in mind ?

function OCPSolutionFromDOCP_raw(docp, solution; objective=nothing, constraints_violation=nothing, iterations=0, multipliers_con=nothing, multipliers_L=nothing, multipliers_U=nothing, message=nothing)

[jbcaillau]

@PierreMartinon @gergaud @ocots Regarding optimisation solvers: as for plots (WIP on recipes), it would be nicer (and seems possible) to keep the solver outside our packages

• to have something more flexible / less depending on one particular solver (this is good for benchmarking)
• to have lighter dependencies

For instance, in the spirit of OCPSolutionFromDOCP_raw in the previous post:

using OptimalControl
using ADNLPModels

@def ocp begin
    ...
end

docp = discretise(ocp, model=:adnlp)
dsol = solve(docp)
sol = foo(ocp, dsol) # needs a better name 🥲
plot(sol)

NB. This is strongly related to having a proper interface defining what traits a "discrete optimal control problem" should have.

[PierreMartinon]

It is certainly possible, the current way is quite close.

docp = directTranscription(ocp, grid_size=100)
sol = solveDOCP(docp, print_level=5, tol=1e-12)

Your foo is actually called OCPSolutionFromDOCP(docp, docp_solution), and we could move the call outside of solveDOCP

function solveDOCP(docp::DOCP; init=nothing, display::Bool=__display(),print_level::Integer=__print_level_ipopt(),mu_strategy::String=__mu_strategy_ipopt(),kwargs...)

    # solve DOCP with NLP solver
    ... IPOPT call with or without starting guess
    ... omitted for clarity

    # return solution for original OCP
    return OCPSolutionFromDOCP(docp, docp_solution)
end

I also keep a solveDirect(ocp, ...) that does everything in one call and returns the OCP solution.

[jbcaillau]

@PierreMartinon I think we need to keep track (= pass as an arg) the original ocp to build a proper optimal control solution (e.g. to know that the variable, state and control dimensions, their names for further plotting, etc.) hence the signature

sol = OCPSolutionFromDOCP(ocp, docp)

I am assuming that docp is strictly suited for the chosen optimal solver and contains no additional information. It might possible to enrich the structure while keeping it compatible with what the optimisation solver requires; not a crucial point, though.

Regarding solveDOCP, I think we should not provide it:

• providing a docp (with directTranscription parametrised by the targeted optimisation solver) is enough
• then solving the docp itself is done outside OptimalControl.jl
• might still be nice to keep a generic solve function in our package with a default choice for optimisation model and solver (ADNLP+ iPopt...) to provide the def + solve + plot basic chain of operations (instead of def + ocp2docp + solve + docp2ocpsol + plot is not bad and much more flexible)
• what is gained wrt. having the choice of optimiser as an option of our solve is that, for a given type of model, there are already many solvers available (and there will be more)

NB. The name is much better but too long 🥲. Should check what the standard for conversion functions in Julia is (typeA2typeB, typeBfromtypeA, typeBoftypeA...)

[PierreMartinon]

A few precisions regarding DOCP:

• it contains a copy of the original OCP, used both for the construction of the discretized problem and the generation of the continuous solution after optimization
• it contains a NLPModel that can be passed to any compatible solver. Currently the problem is built by ADNLPModel but there could be other options. Similarly, the current solver is Ipopt but other choices could be possible.

As for the solve, the question of the interface is still open yes. I guess more use cases are needed to better see what is the most practical. Although I think we definitely need to keep at least some form of a solve() function. Users can still use their own solving method on the NLP from DOCP, and call the solution generation afterwards.

For development I will keep a solve somewhere in CTDirect since I need to be able to solve a problem using just CTBase and CTDirect.

[jbcaillau]

Hi @paraynaud ; any feedback on #74 (comment) above?

[paraynaud]

Hi @jbcaillau, I've been quite busy by personal stuffs lately.
This is unrelated to research, but it takes most of my time.
I did not forget the tuto, but I missed the time and the bandwidth to make it happen for now.

[jbcaillau]

@paraynaud no worry! courage. quite a few things in order (check control-toolbox/CTBase.jl#173), we keep in touch.