Pb with zero in dyn rhs?

[jbcaillau]

@PierreMartinon

@def ocp2 begin
    s ∈ [0, 1], time
    y = (x, tf) ∈ R², state
    u ∈ R, control
    dx = tf(s)^2 * u(s)
    dtf = 0 # while 0 * u(s) is fine 😱 (1)
    ẏ(s) == [dx, dtf]
    x(0) == 0
    x(1) == 1
    tf(1) + 0.5∫(tf(s) * u(s)^2) → min
end

(1) check #189

[jbcaillau]

@PierreMartinon Pb with ForwardDiff when evaluating to solve:

julia> @def ocp2 begin
           s ∈ [0, 1], time
           y = (x, tf) ∈ R², state
           u ∈ R, control
           dx = tf(s)^2 * u(s)
           dtf = 0
           ẏ(s) == [dx, dtf]
           x(0) == 0
           x(1) == 1
           tf(1) + 0.5∫(tf(s) * u(s)^2) → min
       end

The (autonomous) optimal control problem is given by:

    s ∈ [0, 1], time
    y = ((x, tf) ∈ R², state)
    u ∈ R, control
    dx = tf(s) ^ 2 * u(s)
    dtf = 0
    ẏ(s) == [dx, dtf]
    x(0) == 0
    x(1) == 1
    tf(1) + 0.5 * ∫(tf(s) * u(s) ^ 2) → min

The (autonomous) optimal control problem is of the form:

    minimize  J(y, u) = g(y(0), y(1)) + ∫ f⁰(y(s), u(s)) ds, over [0, 1]

    subject to

        ẏ(s) = f(y(s), u(s)), s in [0, 1] a.e.,

        ϕl ≤ ϕ(y(0), y(1)) ≤ ϕu, 

    where y(s) = (x(s), tf(s)) ∈ R² and u(s) ∈ R.

Declarations (* required):

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


julia> solve(ocp2)
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:     1403
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      303

Total number of variables............................:      404
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      303
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"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}, Float64, 1}}, ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}, Float64, 1}} and ForwardDiff.Tag{ADNLPModels.var"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}
Stacktrace:
  [1] value
    @ ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:94 [inlined]
  [2] (::ForwardDiff.var"#76#77"{ForwardDiff.Tag{…}})(d::ForwardDiff.Dual{ForwardDiff.Tag{…}, ForwardDiff.Dual{…}, 1})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:6

[PierreMartinon]

More generally, same error if you put a constant instead of 0. Maybe this is AD related ?

Not really a problem for us since we would typically use a variable instead of a constant state, but it is a bit annoying.

[jbcaillau]

More generally, same error if you put a constant instead of 0. Maybe this is AD related ?

Not really a problem for us since we would typically use a variable instead of a constant state, but it is a bit annoying.

@PierreMartinon yes, the error seems to come from ForwardDiff. (Actually wasted my early morning on it 🙁.) Try to investigate a bit, maybe worth sending an issue to them once the problem is clearly located.

[jbcaillau]

@amontoison does the error below ring a bell? Seems to happen when calling ForwardDiff on functions involving constants, as in

f(x, u) = [x[1]^2 * u, 0] # crashes while [x[1]^2 * u, 0 * u] is fine

ERROR: Cannot determine ordering of Dual tags ForwardDiff.Tag{ReverseDiff.var"#130#133"{typeof(+), ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}, Float64, 1}}, ForwardDiff.Dual{ForwardDiff.Tag{ADNLPModels.var"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}, Float64, 1}} and ForwardDiff.Tag{ADNLPModels.var"#ψ#73"{CTDirect.var"#37#39"{CTDirect.DOCP}}, Float64}

@gdalle I see that you've had a look at JuliaDiff/ForwardDiff.jl#320 : any hint?

[amontoison]

The issue is the combination of ForwardDiff.jl with ReverseDiff.jl for the sparse Hessian.
Can you try to create the ADNLPModel without backend = :optimized to confirm it?
By default, ADNLPModels.jl will do ForwardDiff.jl over ForwardDiff.jl for the Hessian (slower but more robust).

I would like to switch to DI.jl and only rely on Enzyme.jl in the future. Enzyme.jl supports both reverse and forward.

[gdalle]

I assume either ADNLPModels or CTBase must be creating some custom tags somewhere which cause this issue?
If you want to keep abreast of current events, the switch to DI in ADNLPModels is being discussed in JuliaSmoothOptimizers/ADNLPModels.jl#302, but it's gonna be painful

[jbcaillau]

@PierreMartinon @ocots @amontoison @gdalle one has to use zero to propagate the dual number (🙏🏽 @abavoil for the hint) :

ocp2 = @def begin
    s ∈ [0, 1], time
    y = (x, tf) ∈ R², state
    u ∈ R, control
    dx = tf(s)^2 * u(s)
    dtf = zero(y(s))[1] # will propagate dual number
    ẏ(s) == [dx, dtf]
    x(0) == 0
    x(1) == 1
    tf(1) + 0.5∫(tf(s) * u(s)^2) → min
end

solve(ocp2)

[jbcaillau]

NB. Check #189 (comment) for the complete example (and converging solve!)