-
Notifications
You must be signed in to change notification settings - Fork 6
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
Ocp, docp, nlp models... #115
Comments
What do you mean ? |
@PierreMartinon what I meant 🙂 Use case 1using 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 2using 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 using Plots # plot is available (extension of CTBase.jl)
sol = build(solution(nlp_sol), docp) # sol::OCPSolution, more or less OCPSolutionFromDOCP...
plot(sol) |
@jbcaillau Regarding the Use Case 2, using current functions it would look like
So I would need to clean up a bit and combine the part
into something like
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) |
More precisely, we can already split the calls as in
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...
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.
|
Also linked to the older #74 |
@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 Footnotes
|
@jbcaillau @ocots Yes this is what I went for. 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 ? |
@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 there might indeed be some variability in the way constraints are treated / multipliers are ordered by the NLP solver. |
See control-toolbox/OptimalControl.jl#207 and #154 Have a nice weeekend ! |
Can we close? |
@ocots @PierreMartinon about that, @PierreMartinon are you done with changes in
|
The only change in the API is that the
while (re)setting the initial guess in the
The change is only visible when handling the NLP directly and has no impact on the standard The only updates needed in
(as usual, releases will need to be done together for |
If the api has changed the new release of CTDirect should have a different Y. |
Yes, that would be 0.10.0 I guess. |
@ocots @PierreMartinon @0Yassine0 @frapac first try 🙂, Goddard problem julia> @def ocp begin # definition of the optimal control problem
tf ∈ R, variable
t ∈ [ t0, tf ], time
x = (r, v, m) ∈ R³, state
u ∈ R, control
x(t0) == [ r0, v0, m0 ]
m(tf) == mf, (1)
0 ≤ u(t) ≤ 1
r(t) ≥ r0
0 ≤ v(t) ≤ vmax
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
r(tf) → max
end
julia> docp = direct_transcription(ocp; grid_size=5000)
julia> nlp = get_nlp(docp)
ADNLPModel - Model with automatic differentiation backend ADModelBackend{
ReverseDiffADGradient,
ReverseDiffADHvprod,
ForwardDiffADJprod,
ReverseDiffADJtprod,
SparseADJacobian,
SparseReverseADHessian,
ForwardDiffADGHjvprod,
}
Problem name: Generic
All variables: ████████████████████ 20005 All constraints: ████████████████████ 15004
free: ██████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 5002 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: █████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 5001 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
low/upp: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 10002 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ████████████████████ 15004
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: ( 99.97% sparsity) 55011 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ████████████████████ 15004
nnzj: ( 99.97% sparsity) 105004
Counters:
obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
julia> @time ipopt(nlp)
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.
Number of nonzeros in equality constraint Jacobian...: 105004
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 55011
Total number of variables............................: 20005
variables with only lower bounds: 5001
variables with lower and upper bounds: 10002
variables with only upper bounds: 0
Total number of equality constraints.................: 15004
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
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.0100000e+00 9.00e-01 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0096633e+00 8.66e-01 8.39e+02 -1.0 2.60e+00 - 3.18e-03 3.81e-02f 1
2 -1.0107757e+00 7.47e-01 2.04e+03 -1.7 1.03e+00 - 3.22e-02 1.38e-01f 1
3 -1.0113617e+00 5.12e-01 1.01e+03 -1.7 7.47e-01 - 1.31e-01 3.14e-01f 1
4 -1.0093886e+00 3.01e-01 2.26e+03 -1.7 5.12e-01 - 4.02e-01 4.12e-01f 1
5 -1.0066535e+00 3.01e-03 4.85e+03 -1.7 3.01e-01 - 4.98e-03 9.90e-01h 1
6 -1.0079763e+00 3.01e-05 4.52e+02 -1.7 2.89e-01 - 2.33e-01 9.90e-01f 1
7 -1.0083752e+00 4.10e-06 1.43e+04 -1.7 2.44e-01 - 7.47e-01 9.93e-01f 1
8 -1.0083340e+00 4.45e-08 3.08e+00 -2.5 3.46e-02 - 1.00e+00 1.00e+00h 1
9 -1.0083454e+00 8.05e-09 2.72e-01 -3.8 7.88e-03 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -1.0083655e+00 8.40e-10 1.38e-02 -5.7 2.11e-03 - 1.00e+00 1.00e+00h 1
11 -1.0096160e+00 3.12e-07 3.57e+01 -8.6 5.98e-02 - 8.06e-01 1.00e+00h 1
12 -1.0119286e+00 2.37e-06 1.16e+01 -8.6 1.60e-01 - 6.74e-01 8.14e-01h 1
13 -1.0123971e+00 9.85e-07 1.42e+00 -8.6 1.05e-01 - 8.78e-01 6.88e-01h 1
14 -1.0125381e+00 4.57e-07 2.58e-01 -8.6 1.26e-01 - 8.18e-01 8.14e-01h 1
15 -1.0125630e+00 1.51e-07 3.05e-02 -8.6 1.04e-01 - 8.82e-01 8.73e-01f 1
16 -1.0125660e+00 5.50e-08 1.02e-03 -8.6 7.60e-02 - 9.67e-01 9.69e-01f 1
17 -1.0125660e+00 2.82e-08 1.27e-04 -8.6 7.61e-02 - 1.00e+00 5.00e-01f 2
18 -1.0125726e+00 2.57e-08 3.54e-07 -9.0 5.19e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 18
(scaled) (unscaled)
Objective...............: -1.0125726087110427e+00 -1.0125726087110427e+00
Dual infeasibility......: 3.5366498671414553e-07 3.5366498671414553e-07
Constraint violation....: 7.6224774270272633e-09 2.5654078028569671e-08
Variable bound violation: 1.1132245658748600e-31 1.1132245658748600e-31
Complementarity.........: 1.2766492761943900e-09 1.2766492761943900e-09
Overall NLP error.......: 7.6224774270272633e-09 3.5366498671414553e-07
Number of objective function evaluations = 20
Number of objective gradient evaluations = 19
Number of equality constraint evaluations = 20
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 19
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 18
Total seconds in IPOPT = 11.305
EXIT: Optimal Solution Found.
11.312735 seconds (35.63 M allocations: 1.828 GiB, 19.71% gc time)
"Execution stats: first-order stationary"
julia> @time madnlp(nlp)
This is MadNLP version v0.8.4, running with umfpack
Number of nonzeros in constraint Jacobian............: 105004
Number of nonzeros in Lagrangian Hessian.............: 55011
Total number of variables............................: 20005
variables with only lower bounds: 5001
variables with lower and upper bounds: 10002
variables with only upper bounds: 0
Total number of equality constraints.................: 15004
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
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.0100000e+00 9.00e-01 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0096633e+00 8.66e-01 4.25e+01 -1.0 2.60e+00 - 3.18e-03 3.81e-02h 1
2 -1.0100423e+00 7.03e-01 6.82e+01 -1.0 8.66e-01 - 4.11e-02 1.88e-01h 1
3 -1.0083987e+00 5.14e-01 8.02e+00 -1.0 7.03e-01 - 4.71e-01 2.69e-01h 1
4 -1.0061304e+00 5.14e-03 1.01e+02 -1.0 5.14e-01 - 2.20e-02 9.90e-01h 1
5 -1.0081414e+00 1.52e-03 2.59e+01 -1.0 4.28e-01 - 8.87e-01 7.04e-01h 1
6 -1.0084646e+00 1.51e-05 1.38e+00 -1.0 1.43e-01 - 1.00e+00 9.90e-01h 1
7 -1.0085001e+00 1.91e-07 1.35e+01 -1.0 3.29e-02 - 1.00e+00 9.96e-01h 1
8 -1.0084940e+00 6.88e-10 1.82e-03 -1.0 1.34e-02 - 1.00e+00 1.00e+00h 1
9 -1.0084940e+00 3.33e-16 1.25e+04 -3.8 2.50e-06 - 9.99e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -1.0085019e+00 1.44e-11 3.37e-06 -3.8 7.52e-04 - 1.00e+00 1.00e+00h 1
11 -1.0085172e+00 1.08e-10 3.94e+01 -8.6 1.59e-03 - 9.97e-01 1.00e+00h 1
12 -1.0114077e+00 3.68e-06 1.52e+01 -8.6 2.22e-01 - 6.13e-01 8.64e-01h 1
13 -1.0121777e+00 1.94e-06 2.69e+00 -8.6 1.35e-01 - 8.24e-01 6.12e-01h 1
14 -1.0124634e+00 8.75e-07 3.72e-01 -8.6 1.38e-01 - 8.62e-01 7.13e-01h 1
15 -1.0125548e+00 3.65e-07 4.57e-02 -8.6 1.21e-01 - 8.87e-01 8.54e-01h 1
16 -1.0125648e+00 1.16e-07 6.41e-03 -8.6 9.21e-02 - 8.65e-01 8.52e-01h 1
17 -1.0125660e+00 2.89e-08 5.57e-04 -8.6 1.08e-01 - 9.14e-01 9.11e-01h 1
18 -1.0125660e+00 4.16e-09 6.65e-09 -8.6 8.27e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 18
(scaled) (unscaled)
Objective...............: -1.0125660371054928e+00 -1.0125660371054928e+00
Dual infeasibility......: 6.6492546317747740e-09 6.6492546317747740e-09
Constraint violation....: 4.1625690672120186e-09 4.1625690672120186e-09
Complementarity.........: 3.1425345666755113e-09 3.1425345666755113e-09
Overall NLP error.......: 6.6492546317747740e-09 6.6492546317747740e-09
Number of objective function evaluations = 19
Number of objective gradient evaluations = 19
Number of constraint evaluations = 19
Number of constraint Jacobian evaluations = 19
Number of Lagrangian Hessian evaluations = 18
Total wall-clock secs in solver (w/o fun. eval./lin. alg.) = 1.733
Total wall-clock secs in linear solver = 4.475
Total wall-clock secs in NLP function evaluations = 6.742
Total wall-clock secs = 12.950
EXIT: Optimal Solution Found (tol = 1.0e-08).
12.951147 seconds (35.64 M allocations: 39.980 GiB, 23.45% gc time)
"Execution stats: Optimal Solution Found (tol = 1.0e-08)." |
@PierreMartinon for the record: on the current release of CTDirect (v0.10) and associated syntax: to what extent do you need the |
@PierreMartinon please answer 🥹 and close! |
Not sure I understand your question: you only need the solution for plot. You do need the docp to build the ocp functional solution from the discrete docp solution. For instance from https://github.com/control-toolbox/CTDirect.jl/blob/bolza/test/suite/test_misc.jl
You can then do Or is it something else ? Ps. The code above tests some internal functions for the NLP handling, normally you would just run |
@PierreMartinon OK. Fine! Closing. |
@PierreMartinon still some things to be discussed in the process
ocp -> docp -> sol -> plot
. Check this old issueThe text was updated successfully, but these errors were encountered: