Ocp, docp, nlp models...

[jbcaillau]

@PierreMartinon still some things to be discussed in the process ocp -> docp -> sol -> plot. Check this old issue

[PierreMartinon]

What do you mean ?

[jbcaillau]

@PierreMartinon what I meant 🙂

Use case 1

using OptimalControl
using NLPModelsIpopt # solve (with Ipopt) is available (extension of CTDirect.jl)
using Plots # plot is available (extension of CTBase.jl)

@def ocp ...
sol = solve(ocp) # sol::OCPSolution
plot(sol)

Use case 2

using OptimalControl
using MyNlpSolver # has its own solve for an NLP (e.g., an ADNLPModel)

@def ocp ...
docp = DirectTranscription(ocp; model=:adnlp) # docp::ADNLPModel
nlp_sol = solve(docp)

Use case 2 (cont'ed)

Enforcing some traits (API...) for the solution return by an NLP solver, namely enforcing that the returned type has solution getter (returning a vector), one can then write:

using Plots # plot is available (extension of CTBase.jl)

sol = build(solution(nlp_sol), docp) # sol::OCPSolution, more or less OCPSolutionFromDOCP...
plot(sol)

[PierreMartinon]

@jbcaillau Regarding the Use Case 2, using current functions it would look like

@def ocp...
docp = directTranscription(ocp)
nlp_sol = solve(docp) 
# here I assume nlp_sol contains at least the vector of the primal variables at convergence, and maybe multipliers (see next message)
lots_of_stuff = parse_DOCP_solution(docp, nlp_sol)
sol = OCPSolutionFromDOCP_raw(docp, lots_of_stuff)
plot(sol)

So I would need to clean up a bit and combine the part

lots_of_stuff = parse_DOCP_solution(docp, nlp_sol)
sol = OCPSolutionFromDOCP_raw(docp, lots_of_stuff)

into something like

sol = OCPSolutionFromNLPSolution(docp, nlp_sol)

All in all, I think we're pretty close. Should not be too hard to finalize. The one point may be the exact format of nlp_sol: it's easy enough if we consider only the primal solution, but if we want multipliers too we'll need to choose a format (see message below)

[PierreMartinon]

More precisely, we can already split the calls as in

docp = directTranscription(ocp)
dsol = solve(docp, print_level=0, tol=1e-12)
sol1 = OCPSolutionFromDOCP(docp, dsol)

For a custom NLP solver we would typically get something a bit different from the dsol above (this one is an ipopt solution from NLPModelIpopt). We currently have the more low level subfunction actually used to build the OCP solution using only vectors from the NLP solution, mainly T,X,U,v and P. The function still requires the docp, but we have it in this use case.

And yes, I really need to combine the arguments somehow...

function OCPSolutionFromDOCP_raw(docp, T, X, U, v, P;
    objective=0, iterations=0, constraints_violation=0,
    message="No msg", stopping=nothing, success=nothing,
    sol_control_constraints=nothing, sol_state_constraints=nothing, sol_mixed_constraints=nothing, sol_variable_constraints=nothing, mult_control_constraints=nothing, mult_state_constraints=nothing, mult_mixed_constraints=nothing, mult_variable_constraints=nothing, mult_state_box_lower=nothing, 
    mult_state_box_upper=nothing, mult_control_box_lower=nothing, mult_control_box_upper=nothing, mult_variable_box_lower=nothing, mult_variable_box_upper=nothing)

To parse the NLP solution into this awful bunch of different variables and multipliers we have the following function, in which 'solution' is the NLP solution vector (primal variables). I think I could make the multipliers and constraints optional for users who are just interested in the optimal trajectory.

Note that taking the constraints and multipliers into account is a bit more involved since each solver could have its own way of handling the constraints (eg separate between boxes/ranges, linear, general nonlinear...). The current format follows ipopt, with box and nonlinear constraints.

parse_DOCP_solution(docp, solution, multipliers_constraints, multipliers_LB, multipliers_UB, constraints)

[PierreMartinon]

Also linked to the older https://github.com/control-toolbox/CTDirect.jl/issues/74

[jbcaillau]

@jbcaillau Regarding the Use Case 2, using current functions it would look like

@def ocp...
docp = directTranscription(ocp)
nlp_sol = solve(docp) 
# here I assume nlp_sol contains at least the vector of the primal variables at convergence, and maybe multipliers (see next message)
lots_of_stuff = parse_DOCP_solution(docp, nlp_sol)
sol = OCPSolutionFromDOCP_raw(docp, lots_of_stuff)
plot(sol)

So I would need to clean up a bit and combine the part [...]

@PierreMartinon Regarding what is returned by an arbitrary NLP solver, it is reasonable to expect that the user is able to access the primal vars (kind of basic trait that we want any solver to have) and, possibly, to the dual (= multipliers - if not, we just don't plot them). One way to enforce this can be:[^1]

using OptimalControl
using MyNLPSolver

@def ocp ...
docp = direct_transcription(ocp; model=:my_model)
nlp_sol = solve(docp) # Could be solve!(docp) depending on the solver
sol = build(docp; primal=nlp_sol.x, dual=nlp_sol.λ) # A simple version of parse_DOCP_solution
print(sol)

where build (a simplified version of OCPSolutionFromDOCP_raw, with a shorter name 🤭) only assumes primal and dual solutions are passed in the same order as what is defined in docp.

[^1]: Julia Blue style recommends to limit upper camel-case convention to modules and types (and lowercase and underscores for functions)

[PierreMartinon]

@jbcaillau @ocots Yes this is what I went for. Currently only the primal, ie sol=ocp_solution_from_nlp(docp, solution) but the extension to sol=ocp_solution_from_nlp(docp, solution, multipliers) is almost ready.

The part that will be a bit more involved is the range constraints, since they are handled separately and may be less standardized among NLP solvers. An overhaul will be needed here anyway, since this is the section with the gazillion variables. I should have a look at other solver interfaces before I rewrite this, for better portability.

Ps. the contest for a better name is open :D Maybe just build_solution ?

[jbcaillau]

@PierreMartinon well, I'd say that it's definitely the user's job to retrieve the multipliers properly. But you are right: although we know exactly the structure of what we pass, i.e. an NLP structured as

F(X) \rightarrow \min\\
L_G \leq G(X) \leq U_G\\
L \leq X \leq U

there might indeed be some variability in the way constraints are treated / multipliers are ordered by the NLP solver.

[PierreMartinon]

See https://github.com/control-toolbox/OptimalControl.jl/pull/207 and https://github.com/control-toolbox/CTDirect.jl/pull/154

Have a nice weeekend !