Skip to content

[Bug] Reverse over forward AD issues with ADNLP #481

Open
Open
[Bug] Reverse over forward AD issues with ADNLP#481
Assignees
ocots
Labels
bugSomething isn't working

Description

@jbcaillau
jbcaillau
opened on Apr 27, 2025

@ocots @PierreMartinon Same kind of error with AD as in #384 tough here the problem seems to be directly related to using a constant in the computation: this fails (Cannot determine ordering of Dual tags ForwardDiff.Tag)

using OptimalControl
using NLPModelsIpopt

o = @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x(0) - [-1, v] == [0, 0] # OK: the boundary contraint may involve the variable
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫( 0.5u(t)^2 ) → min
       end

s = solve(o)

this also fails (same error)

o = @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x(0) - [-one(v), v] == [0, 0] # OK: the boundary contraint may involve the variable
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫( 0.5u(t)^2 ) → min
       end

while this is OK

o = @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x(0) - [-1 + 0v, v] == [0, 0] # OK: the boundary contraint may involve the variable
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫( 0.5u(t)^2 ) → min
       end

@PierreMartinon I thought that using one as above (second version; see also zero) solved the problem, but this is not the case. any clue?

Activity

  1. jbcaillau
    added
    bugSomething isn't working
    on Apr 27, 2025
  2. jbcaillau
    assigned
    ocots
    on Apr 27, 2025

  3. jbcaillau commented on Apr 27, 2025

    @jbcaillau
    jbcaillau
    on Apr 27, 2025
    MemberAuthor

    @amontoison @gdalle input welcome 🤞🏽

  4. amontoison commented on Apr 27, 2025

    @amontoison
    amontoison
    on Apr 27, 2025

    @jbcaillau If you can provide the full error, it will help.

  5. jbcaillau commented on Apr 28, 2025

    @jbcaillau
    jbcaillau
    on Apr 28, 2025
    MemberAuthor

    @amontoison full log attached (copy paste of the first cell above)

    @jbcaillau If you can provide the full error, it will help.

    julia> using OptimalControl


julia> using NLPModelsIpopt

julia> o = @def begin
                  v ∈ R, variable
                  t ∈ [0, 1], time
                  x ∈ R², state
                  u ∈ R, control
                  -1 ≤ v ≤ 1
                  x(0) - [-1, v] == [0, 0] # OK: the boundary contraint may involve the variable
                  x(1) == [0, 0]
                  ẋ(t) == [x₂(t), u(t)]
                  ∫( 0.5u(t)^2 ) → min
              end
    v ∈ R, variable
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    -1 ≤ v ≤ 1
    x(0) - [-1, v] == [0, 0]
    x(1) == [0, 0]
    ẋ(t) == [x₂(t), u(t)]
    ∫(0.5 * u(t) ^ 2) → min

The optimal control problem is of the form:

    minimize  J(x, u, v) = ∫ f⁰(t, x(t), u(t), v) dt, over [0, 1]

    subject to

        ẋ(t) = f(t, x(t), u(t), v), t in [0, 1] a.e.,

        ϕ₋ ≤ ϕ(x(0), x(1), v) ≤ ϕ₊, 
        v₋ ≤ v ≤ v₊, 

    where x(t) ∈ R², u(t) ∈ R and v ∈ R.

Declarations (* required):

╭────────┬────────┬──────────┬──────────┬───────────┬────────────┬─────────────╮
│ times* │ state* │ control* │ variable │ dynamics* │ objective* │ constraints │
├────────┼────────┼──────────┼──────────┼───────────┼────────────┼─────────────┤
│   V    │   V    │    V     │    V     │     V     │     V      │      V      │
╰────────┴────────┴──────────┴──────────┴───────────┴────────────┴─────────────╯


julia> s = solve(o)
This is Ipopt version 3.14.17, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     3006
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      251

Total number of variables............................:     1005
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:      755
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

ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{ReverseDiff.var"#130#133"{typeof(-), ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}, Float64, 1}}, ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}, Float64, 1}} and ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}
Stacktrace:
  [1] value
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/dual.jl:98 [inlined]
  [2] (::ForwardDiff.var"#76#77"{ForwardDiff.Tag{…}})(d::ForwardDiff.Dual{ForwardDiff.Tag{…}, ForwardDiff.Dual{…}, 1})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/apiutils.jl:6
  [3] value!(f::ForwardDiff.var"#76#77"{…}, r::DiffResults.ImmutableDiffResult{…}, x::ForwardDiff.Dual{…})
    @ DiffResults ~/.julia/packages/DiffResults/YLo25/src/DiffResults.jl:168
  [4] extract_value!
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/apiutils.jl:5 [inlined]
  [5] derivative!
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/derivative.jl:49 [inlined]
  [6] binary_scalar_forward_exec!(f::typeof(-), output::ReverseDiff.TrackedReal{…}, input::Tuple{…}, cache::Base.RefValue{…})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/rKZaG/src/derivatives/scalars.jl:114
  [7] scalar_forward_exec!(instruction::ReverseDiff.ScalarInstruction{typeof(-), Tuple{…}, ReverseDiff.TrackedReal{…}, Base.RefValue{…}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/rKZaG/src/derivatives/scalars.jl:88
  [8] forward_exec!(instruction::ReverseDiff.ScalarInstruction{typeof(-), Tuple{…}, ReverseDiff.TrackedReal{…}, Base.RefValue{…}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/rKZaG/src/tape.jl:82
  [9] ForwardExecutor
    @ ~/.julia/packages/ReverseDiff/rKZaG/src/api/tape.jl:76 [inlined]
 [10] (::FunctionWrappers.CallWrapper{Nothing})(f::ReverseDiff.ForwardExecutor{ReverseDiff.ScalarInstruction{…}})
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
 [11] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
 [12] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
 [13] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
 [14] forward_pass!(compiled_tape::ReverseDiff.CompiledTape{ReverseDiff.GradientTape{ADNLPModels.var"#ψ#72"{…}, Tuple{…}, ReverseDiff.TrackedReal{…}}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/rKZaG/src/api/tape.jl:124
 [15] seeded_forward_pass!
    @ ~/.julia/packages/ReverseDiff/rKZaG/src/api/tape.jl:42 [inlined]
 [16] gradient!(result::Tuple{Vector{…}, Vector{…}}, tape::ReverseDiff.CompiledTape{ReverseDiff.GradientTape{…}}, input::Tuple{Vector{…}, Vector{…}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/rKZaG/src/api/gradients.jl:79
 [17] ∇l!
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:215 [inlined]
 [18] ∇l!
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:214 [inlined]
 [19] sparse_hess_coord!(b::ADNLPModels.SparseReverseADHessian{…}, x::Vector{…}, obj_weight::Float64, y::SubArray{…}, vals::Vector{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:351
 [20] hess_coord!(b::ADNLPModels.SparseReverseADHessian{…}, nlp::ADNLPModels.ADNLPModel{…}, x::Vector{…}, y::SubArray{…}, obj_weight::Float64, vals::Vector{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:375
 [21] hess_coord!(nlp::ADNLPModels.ADNLPModel{Float64, Vector{…}, Vector{…}}, x::Vector{Float64}, y::Vector{Float64}, vals::Vector{Float64}; obj_weight::Float64)
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nlp.jl:706
 [22] hess_coord!
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/nlp.jl:695 [inlined]
 [23] (::NLPModelsIpopt.var"#eval_h#5"{…})(x::Vector{…}, rows::Vector{…}, cols::Vector{…}, σ::Float64, λ::Vector{…}, values::Vector{…})
    @ NLPModelsIpopt ~/.julia/packages/NLPModelsIpopt/OGzSv/src/NLPModelsIpopt.jl:132
 [24] _Eval_H_CB(n::Int32, x_ptr::Ptr{…}, ::Int32, obj_factor::Float64, m::Int32, lambda_ptr::Ptr{…}, ::Int32, nele_hess::Int32, iRow::Ptr{…}, jCol::Ptr{…}, values_ptr::Ptr{…}, user_data::Ptr{…})
    @ Ipopt ~/.julia/packages/Ipopt/3Ojbq/src/C_wrapper.jl:129
 [25] #5
    @ ~/.julia/packages/Ipopt/3Ojbq/src/C_wrapper.jl:407 [inlined]
 [26] disable_sigint(f::Ipopt.var"#5#6"{Ipopt.IpoptProblem, Base.RefValue{Float64}})
    @ Base ./c.jl:167
 [27] IpoptSolve
    @ ~/.julia/packages/Ipopt/3Ojbq/src/C_wrapper.jl:406 [inlined]
 [28] solve!(solver::IpoptSolver, nlp::ADNLPModels.ADNLPModel{…}, stats::SolverCore.GenericExecutionStats{…}; callback::Function, kwargs::@Kwargs{…})
    @ NLPModelsIpopt ~/.julia/packages/NLPModelsIpopt/OGzSv/src/NLPModelsIpopt.jl:263
 [29] solve!
    @ ~/.julia/packages/NLPModelsIpopt/OGzSv/src/NLPModelsIpopt.jl:184 [inlined]
 [30] #solve!#16
    @ ~/.julia/packages/SolverCore/vTpnV/src/solver.jl:40 [inlined]
 [31] solve!
    @ ~/.julia/packages/SolverCore/vTpnV/src/solver.jl:38 [inlined]
 [32] solve_docp(solver_backend::CTDirect.IpoptBackend, docp::CTDirect.DOCP{…}, nlp::ADNLPModels.ADNLPModel{…}; display::Bool, max_iter::Int64, tol::Float64, print_level::Int64, mu_strategy::String, linear_solver::String, kwargs::@Kwargs{})
    @ CTSolveExtIpopt ~/.julia/packages/CTDirect/rvI7x/ext/CTSolveExtIpopt.jl:55
 [33] solve(::Model{…}; display::Bool, grid_size::Int64, disc_method::Symbol, time_grid::Nothing, init::Nothing, adnlp_backend::Symbol, kwargs::@Kwargs{})
    @ CTDirect ~/.julia/packages/CTDirect/rvI7x/src/solve.jl:75
 [34] solve(::Model{…})
    @ CTDirect ~/.julia/packages/CTDirect/rvI7x/src/solve.jl:38
 [35] solve(::Model{…}; kwargs::@Kwargs{})
    @ OptimalControl ~/.julia/packages/OptimalControl/SiL76/src/solve.jl:96
 [36] solve(::Model{…})
    @ OptimalControl ~/.julia/packages/OptimalControl/SiL76/src/solve.jl:87
 [37] top-level scope
    @ REPL[4]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia>

  6. amontoison commented on Apr 28, 2025

    @amontoison
    amontoison
    on Apr 28, 2025

    Hum it is the forward over reverse mode for the Hessian that fails. If you don't use backend=:optimized, it should work because it's forward over forward by default.

    I should investigate but maybe better to spend time on the support of DI.jl in ADNLPModels.jl. Guillaume is maybe doing something smarter than us with the tape.

  7. PierreMartinon commented on May 2, 2025

    @PierreMartinon
    PierreMartinon
    on May 2, 2025
    Member

    @jbcaillau Do you still have the error with option adnlp_backend=:default added to the solve call ?

  8. ocots commented on May 2, 2025

    @ocots
    ocots
    on May 2, 2025
    Member

    @jbcaillau Do you still have the error with option adnlp_backend=:default added to the solve call ?

    It works (but as you said it is slower so cannot be the default for us).

  9. ocots commented on May 2, 2025

    @ocots
    ocots
    on May 2, 2025
    Member

    Should be explained on the doc.

  10. jbcaillau commented on May 3, 2025

    @jbcaillau
    jbcaillau
    on May 3, 2025
    MemberAuthor

    @PierreMartinon @ocots @amontoison defaulting adnlp_backend=:default indeed works. as pointed by @ocots 6 times slower on that example (o2 is equivalent but raises no bug as not single constant is used, check -1 + 0 * v instead of merely -1 😬).

    @jbcaillau Do you still have the error with option adnlp_backend=:default added to the solve call ?

    It works (but as you said it is slower so cannot be the default for us).

    julia> o1 = @def begin
                         v ∈ R, variable
                         t ∈ [0, 1], time
                         x ∈ R², state
                         u ∈ R, control
                         -1 ≤ v ≤ 1
                         x(0) - [-1, v] == [0, 0] # OK: the boundary contraint may involve the variable
                         x(1) == [0, 0]
                         ẋ(t) == [x₂(t), u(t)]
                         ∫( 0.5u(t)^2 ) → min
                     end

julia> @btime solve(o1; adnlp_backend=:default) # no more exception, indeed
159.537 ms (502305 allocations: 237.32 MiB)

julia> o2 = @def begin
                         v ∈ R, variable
                         t ∈ [0, 1], time
                         x ∈ R², state
                         u ∈ R, control
                         -1 ≤ v ≤ 1
                         x(0) - [-1 + 0 * v, v] == [0, 0] # OK: the boundary contraint may involve the variable
                         x(1) == [0, 0]
                         ẋ(t) == [x₂(t), u(t)]
                         ∫( 0.5u(t)^2 ) → min
                     end

julia> @btime solve(o2)
27.333 ms (113391 allocations: 15.55 MiB) # and same objective reached up to

  11. jbcaillau commented on May 3, 2025

    @jbcaillau
    jbcaillau
    on May 3, 2025
    MemberAuthor

    @ocots adding a Know issues section (pointing to issues) at the end of the "define a problem" tutorial #480

    Should be explained on the doc.

  12. jbcaillau commented on May 3, 2025

    @jbcaillau
    jbcaillau
    on May 3, 2025
    MemberAuthor

    @PierreMartinon @amontoison any of you has tried to use ADNLP on a problem with a constraint function involving a vectorial function having one constant component to replicate the issue?

  13. amontoison commented on May 25, 2025

    @amontoison
    amontoison
    on May 25, 2025

    @jbcaillau
    I was cleaning my emails and checked again this issue. Why do you use some Vector{Real} ?
    It is not a concrete type for Julia.

    ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{ReverseDiff.var"#130#133"{typeof(-), ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}, Float64, 1}}, ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}, Float64, 1}} and ForwardDiff.Tag{ADNLPModels.var"#ψ#72"{CTDirect.var"#14#16"{CTDirect.DOCP{CTDirect.Trapeze, Model{CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}, CTModels.StateModel, CTModels.ControlModel, CTModels.VariableModel, var"##248#3", CTModels.LagrangeObjectiveModel{var"##252#4"}, CTModels.ConstraintsModel{Tuple{Vector{Real}, CTModels.var"#path_cons_nl!#59"{Int64, Vector{Int64}, Tuple{}}, Vector{Real}}, Tuple{Vector{Real}, CTModels.var"#boundary_cons_nl!#61"{Int64, Vector{Int64}, Tuple{var"##242#2", var"##237#1"}}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Tuple{Vector{Real}, Vector{Int64}, Vector{Real}}, Dict{Symbol, Tuple{Symbol, Union{Function, OrdinalRange{<:Int64}}, AbstractVector{<:Real}, AbstractVector{<:Real}}}}}}}}, Float64}

  14. jbcaillau commented on May 25, 2025

    @jbcaillau
    jbcaillau
    on May 25, 2025
    MemberAuthor

    i think everything is eventually set to Float64 and type stable. @ocots and @PierreMartinon do you confirm?

  15. ocots commented on May 25, 2025

    @ocots
    ocots
    on May 25, 2025
    Member

    Where do you that we use Vector{Real}?

  16. 18 remaining items

  17. jbcaillau commented on Feb 25, 2026

    @jbcaillau
    jbcaillau
    on Feb 25, 2026
    MemberAuthor

    @PierreMartinon @ocots the fix above is now merged into main at ADNLModels (see JuliaSmoothOptimizers/ADNLPModels.jl#378). the three problems below solve with it and backend=:enzyme:

    • what about integrating these tests into those of OC.jl / CTDirect.jl?
    • on top of that, would be nice to test more intensively backend=:enzyme to potentially use it as the default (instead of :optimized, currently, that also does forward over reverse but with ReverseDiff.jl) 👀
    using OptimalControl
using NLPModelsIpopt
using Enzyme

backend = :enzyme # vs. :default that does forward over reverse with ReverseDiff and fails

# use case no. 1

x0 = [ 0.
       1. ]
A = [  0. 1.
      -1. 0. ]
B = [ 0.
      1. ]
Q = [ 1. 0.
      0. 1. ]
R = 1.
tf = 3.

o = @def begin
    t ∈ [0, tf], time
    x ∈ R², state
    u ∈ R, control
    x(0) == x0
    ẋ(t) == A * x(t) + B * u(t)
    0.5∫( x(t)' * Q * x(t) + u(t)' * R * u(t) ) → min
end

sol = solve(o; backend=backend)

# use case no. 2

o = @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x(0) - [-1, v] == [0, 0] # OK: the boundary contraint may involve the variable
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫( 0.5u(t)^2 ) → min
       end

sol = solve(o; backend=backend)

# use case no. 3

t0 = 0.
vf = 20.
x0 = [0., 0., 0.1, 0.1]
p = [1., 1., 1., 1., 2., 1.5]

# Functions
e(x) = p[1] + x[3] + x[4]
g(x) = p[2]
f₁₁(x) = p[3]
f₁₂(x) = - p[4]*x[3]
f₂₁(x) = p[5]
f₂₂(x) = - p[6]*x[4]

# Vector fields
Ff(x) = [e(x),  g(x), f₁₁(x), f₂₁(x)]
Fb(x) = [e(x), -g(x), f₁₂(x), f₂₂(x)]
F0(x) = Ff(x) + Fb(x)
F1(x) = Ff(x) - Fb(x)

o = @def begin
    tf ∈ R, variable 
    t ∈ [t0, tf], time
    x = (e, v, r1, r2) ∈ R⁴, state
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(t0) == x0
    v(tf) == vf
    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
    e(tf) → min
end

sol = solve(o; backend=backend)
  18. jbcaillau
    linked a pull request that will close this issue481 bug reverse over forward ad issues with adnlp #706on Feb 25, 2026

  19. jbcaillau commented on Feb 27, 2026

    @jbcaillau
    jbcaillau
    on Feb 27, 2026
    MemberAuthor

  20. amontoison commented on Feb 27, 2026

    @amontoison
    amontoison
    on Feb 27, 2026

    @jbcaillau https://github.com/JuliaSmoothOptimizers/ADNLPModels.jl/blob/main/src/predefined_backend.jl#L61-L66

  21. jbcaillau commented on Feb 27, 2026

    @jbcaillau
    jbcaillau
    on Feb 27, 2026
    MemberAuthor

    @amontoison

    using OptimalControl
using NLPModelsIpopt

backend = :optimized # errors while it is now ok with :enzyme (see above)

# use case no. 2

o = @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x(0) - [-1, v] == [0, 0] # FTR: x(0) - [-1 + 0v, v] == [0, 0] works
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫( 0.5u(t)^2 ) → min
       end

sol = solve(o; backend=backend)

    raises the error documented at the beginning of the issue:

    ▫ OptimalControl v1.2.1-beta solving with: collocation → adnlp → ipopt

  📦 Configuration:
   ├─ Discretizer: collocation (no user options)
   ├─ Modeler: adnlp (backend = optimized)
   └─ Solver: ipopt (no user options)

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:     1755
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      250

Total number of variables............................:      753
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        1
                     variables with only upper bounds:        0
Total number of equality constraints.................:      504
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

ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{ReverseDiff.var"...

    everything comes from the generated function below:

    julia> @macroexpand o = @def begin
                  v ∈ R, variable
                  t ∈ [0, 1], time
                  x ∈ R², state
                  u ∈ R, control
                  -1 ≤ v ≤ 1
                  x(0) - [-1, v] == [0, 0] # OK: the boundary contraint may involve the variable
                  x(1) == [0, 0]
                  ẋ(t) == [x₂(t), u(t)]
                  ∫( 0.5u(t)^2 ) → min
              end
              [...]
              #= REPL[33]:7 =#
              begin
                  #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:186 =#
                  try
                      #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:187 =#
                      begin
                          #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:653 =#
                          function var"fun##108335"(var"r##108338", var"x0##108336", var"xf##108337", var"v##108329")
                              #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:653 =#
                              #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:654 =#
                              var"r##108338"[:] .= var"x0##108336" - [-1, (Base.maybeview)(var"v##108329", 1)]
                              #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:655 =#
                              return nothing
                          end
                          #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:657 =#
                          CTModels.constraint!(var"p_ocp##108320", :boundary; f = var"fun##108335", lb = [0, 0], ub = [0, 0], label = Symbol("label##108333"))
                      end
                  catch
                      #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:189 =#
                      println("Line ", 6, ": ", "x(0) - [-1, v] == [0, 0]")
                      #= /net/home/c/caillau/.julia/packages/CTParser/AmHLr/src/onepass.jl:190 =#
                      rethrow()
                  end
              end

    the function that is accepted replaces x(0) - [-1, v] == [0, 0] with x(0) - [-1 + 0v, v] == [0, 0] so, in the generated code, the line

    var"r##108338"[:] .= var"x0##108336" - [-1, (Base.maybeview)(var"v##108329", 1)]

    is replaced with

    var"r##108423"[:] .= var"x0##108421" - [-1 + 0 * (Base.maybeview)(var"v##108414", 1), (Base.maybeview)(var"v##108414", 1)]

  22. jbcaillau commented on Feb 27, 2026

    @jbcaillau
    jbcaillau
    on Feb 27, 2026
    MemberAuthor
    •  @jbcaillau build a mock up optimisation using the generated function above to define constraints

  23. jbcaillau commented on Apr 27, 2026

    @jbcaillau
    jbcaillau
    on Apr 27, 2026
    MemberAuthor

    cc @abavoil

  24. jbcaillau commented on Apr 27, 2026

    @jbcaillau
    jbcaillau
    on Apr 27, 2026
    MemberAuthor

    @PierreMartinon can you please confirm to @abavoil that the option to use the sparsity structure (that you computed) is :manual, as shown below (excerpt from https://control-toolbox.org/OptimalControl.jl/stable/manual-solve.html)?

    Image

  25. ocots commented on Apr 27, 2026

    @ocots
    ocots
    on Apr 27, 2026
    Member

    Yep see https://control-toolbox.org/Tutorials.jl/stable/tutorial-discretisation.html#Large-problems-and-AD-backend

  26. ocots
    added 2 commits that reference this issue on Oct 9, 2026

    docs(abstract-syntax): structure section, ASCII table, labels, verifi…

    d4331e3

    docs(probes): check the claims the docs build cannot execute

    1a3db34

  27. ocots commented on Oct 9, 2026

    @ocots
    ocots
    on Oct 9, 2026
    Member

    ADNLPModels 0.8.14 fixed the boundary-constraint cases above: x(0) - [-1, v] == [0, 0] and x(0) - [-one(v), v] == [0, 0] now solve with :adnlp. Another case of the same error remains, with ADNLPModels 0.8.14 as with 0.8.13: a dynamics component that is a constant.

    using OptimalControl, NLPModelsIpopt
ocp = @def begin
    t ∈ [0, 1], time
    x ∈ R³, state
    u ∈ R, control
    x(0) == [-1, 0, 2]
    x₁(1) == 0
    x₂(1) == 0
    ẋ(t) == [x₂(t), u(t), 0]
    0.5∫(u(t)^2) → min
end
solve(ocp)          # ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{ReverseDiff…}
solve(ocp, :exa)    # works, J = 6
    Third component of ẋ(t) :adnlp
    0 ✗ Dual tags error
    0.0 ✗
    zero(u(t)) ✗
    1 ✗
    0 * u(t) ✓ J = 6

    So any component that does not depend on the inputs fails, whatever its value, and 0 * u(t) is the workaround (it was already the workaround of CTDirect#190, dtf = 0, closed in January 2025). The same happens when the constant state then appears in the cost or in the dynamics (ẋ = [x₂ u, 0]). This is the usual way to turn a free final time into a state (ṫ_f = 0), so users hit it.

    The tags in the error (ForwardDiff duals inside a ReverseDiff tape, through ADNLPModels' ψ) point to the second-order AD backend, as in the original report; I have not narrowed it down further. Versions: OptimalControl 2.2.1-beta, CTDirect 1, ADNLPModels 0.8.14 (also 0.8.13), Julia 1.12.6.

  28. amontoison commented on Oct 10, 2026

    @amontoison
    amontoison
    on Oct 10, 2026

    @jbcaillau @ocots @PierreMartinon New release 0.8.15 of ADNLPModels.jl.
    It should fix the last issue. 🤞

  29. ocots commented on Oct 10, 2026

    @ocots
    ocots
    on Oct 10, 2026
    Member

    Thank you @amontoison for 0.8.15! I checked it: with ADNLPModels 0.8.15, the last case of this issue now solves with the default :adnlp modeler, a constant component of the dynamics:

    ocp = @def begin
    t ∈ [0, 1], time
    x ∈ R³, state
    u ∈ R, control
    x(0) == [-1, 0, 2]
    x₁(1) == 0
    x₂(1) == 0
    ẋ(t) == [x₂(t), u(t), 0]
    0.5∫(u(t)^2) → min
end
solve(ocp)   # J = 6.000096, same as with `0 * u(t)` and with :exa

    "Cannot determine ordering of Dual tags" is gone. A constant 1.0 component works too. The warning on the Abstract syntax page of the documentation is removed in #971, and the probe that tracked this bug now checks that it stays fixed. Thanks again for the quick fixes, 0.8.14 and now 0.8.15.

  30. ocots
    added a commit that references this issue on Oct 10, 2026

    docs(changelog): constant component of the dynamics, with the example (…

    9656013
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Type

No type

Projects

No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions