[Bug] Goddard 3D does not compile

[horasio]

Dear CT community,

Describe the bug
If the F0() function for the dynamics just returns a vector of zeros, then it compiles.
But if I copy the velocity state into the time-deriv (x[4] in the line below marked as "BUG HERE") then the code does not compile, apparently in some automatic differentiation part.
To Reproduce
Execute the following:

module Goddard3D

using OptimalControl
#using NLPModelsIpopt
#using OrdinaryDiffEq  # to get the Flow function from OptimalControl
#using NonlinearSolve  # interface to NLE solvers
#using MINPACK         # NLE solver: use to solve the shooting equation
using Plots           # to plot the solution

t0 = 0      # initial time
r0 = [0.999949994, 0.0001, 0.01]      # initial altitude
v0 = [0,0,0]      # initial speed, v = 1 ==> 7910 m/s (orbital velocity at sea level)
m0 = 1      # initial mass
x0 = [ r0..., v0..., m0 ]
rf = [1.01,0,0] # r[1] = 1 ==> 6378 km

@def ocp begin # definition of the optimal control problem
    tf ∈ R, variable
    t ∈ [ t0, tf ], time
    x ∈ R⁷, state
    u ∈ R³, control

    r = x[1:3]
    v = x[4:6]
    m = x[7]

    x(t0) == x0
    r(tf) == rf,         (1)
    u(t)' * u(t) ≤ 1.0
    r(t)' * r(t) ≥ 1.0
    #0 ≤ v(t) ≤ vmax

    ẋ(t) == F0(x(t)) + F1(x(t), u(t))

    m(tf) → max

end;

# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 7

F0(x) = begin
    return [x[4],0,0,0,0,0,0] ###  <---- ************** BUG HERE
    #r = x[1:3]
    #v = x[4:6]
    #m = x[7]
    #rnorm2 = r' * r
    #rnorm = sqrt(rnorm2)
    #rnormalized = (1/rnorm) * r
    #vnorm = sqrt(v' * v)
    #D = ( Cd * vnorm * exp(-β*(rnorm - 1)) ) * v # Drag force
    #return [v..., (-D/m - (1/rnorm2)*rnormalized)..., 0 ]
end

F1(x,u) = begin
    m = x[7]
    unorm = sqrt(u' * u)
    t = (Tmax / m) * u
    return [ 0,0,0, t..., -b*Tmax*unorm ]
end

# DIRECT METHOD

function run_direct()
    direct_sol = solve(ocp, grid_size=100)
    plt = plot(direct_sol, size=(900, 900))
end

# INDIRECT SHOOTING METHOD

end # module

Goddard3D.run_direct()

Expected behavior
Should compile

Julia 1.10.5 on MacOS Ventura 13.7 up-to-date packages

[ocots]

Actually, even for this simple problem, there is a bug:

using OptimalControl
using NLPModelsIpopt

ocp = @def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    x(0) == [ -1, 0 ]
    x(1) == [ 0, 0 ]
    ẋ(t) == [ x₂(t), 0 ] + [ 0, u(t) ] # this compiles: ẋ(t) == [ x₂(t), u(t) ]
    ∫( 0.5u(t)^2 ) → min
end

sol = solve(ocp)

We get the error:

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

This compiles:

ẋ(t) == [ x₂(t), u(t) ]
# or 
ẋ(t) == [ x₂(t), 0 ] + [ x₂(t), u(t) ]
# or
ẋ(t) == [ x₂(t), 0 ] + [ x₁(t), u(t) ]
# or 
ẋ(t) == [ x₂(t), 0 ] + u(t) * [ 0, 1 ]

[ocots]

This for instance compiles but there is still an issue:

module Goddard3D

using OptimalControl
using NLPModelsIpopt
#using OrdinaryDiffEq  # to get the Flow function from OptimalControl
#using NonlinearSolve  # interface to NLE solvers
#using MINPACK         # NLE solver: use to solve the shooting equation
using Plots           # to plot the solution

t0 = 0                              # initial time
r0 = [0.999949994, 0.0001, 0.01]    # initial altitude
v0 = [0, 0, 0]                      # initial speed, v = 1 ==> 7910 m/s (orbital velocity at sea level)
m0 = 1                              # initial mass
x0 = [r0..., v0..., m0]
rf = [1.01, 0, 0]                   # r[1] = 1 ==> 6378 km

ocp = @def begin # definition of the optimal control problem

    tf ∈ R, variable
    t ∈ [ t0, tf ], time
    x ∈ R⁷, state
    u ∈ R³, control

    r = x[1:3]
    v = x[4:6]
    m = x[7]

    x(t0) == x0
    r(tf) == rf

    u(t)' * u(t) ≤ 1
    r(t)' * r(t) ≥ 1

    ẋ(t) == F(x(t), u(t))

    m(tf) → max

end;

# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 7

mysqrt(x) = sqrt(x + 1e-6)

F(x, u) = begin

    r = x[1:3]
    v = x[4:6]
    m = x[7]

    rnorm2 = r' * r
    rnorm = mysqrt(rnorm2)
    rnormalized = (1/rnorm) * r
    vnorm = mysqrt(v' * v)
    D = ( Cd * vnorm * exp(-β*(rnorm - 1)) ) * v # Drag force

    ṙ = v
    v̇ = -D/m - (1/rnorm2)*rnormalized + (Tmax / m) * u
    ṁ = -b*Tmax*mysqrt(u' * u)

    return [ṙ..., v̇..., ṁ]

end    

# DIRECT METHOD

function run_direct()
    direct_sol = solve(ocp; grid_size=100)
    plt = plot(direct_sol, size=(900, 900))
end

# INDIRECT SHOOTING METHOD

end # module

Goddard3D.run_direct()

returns

julia> include("goddard.jl");
WARNING: replacing module Goddard3D.
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     7510
Number of nonzeros in inequality constraint Jacobian.:      606
Number of nonzeros in Lagrangian Hessian.............:     4747

Total number of variables............................:     1011
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      710
Total number of inequality constraints...............:      202
        inequality constraints with only lower bounds:      101
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:      101

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -1.0000000e-01 7.98e+137 1.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 1
  Increasing icntl[13] from 1000 to 2000.
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 2
  Increasing icntl[13] from 2000 to 4000.
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 3
  Increasing icntl[13] from 4000 to 8000.
WARNING: Problem in step computation; switching to emergency mode.
   1r-1.0000000e-01 7.98e+137 9.99e+02 129.9 0.00e+00  20.0 0.00e+00 0.00e+00R  1
Restoration phase failed with unexpected solverreturn status 9

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:  -1.0000000000000001e-01   -1.0000000000000001e-01
Dual infeasibility......:   1.0000000000000000e+00    1.0000000000000000e+00
Constraint violation....:  7.9803725102132858e+129   7.9803725102132855e+137
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   9.7000000999999991e-01    9.7000000999999991e-01
Overall NLP error.......:  7.9803725102132858e+129   7.9803725102132855e+137


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 2
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 2
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 1.988

EXIT: Restoration Failed!

[ocots]

@jbcaillau, @PierreMartinon Any idea?

[jbcaillau]

@jbcaillau, @PierreMartinon Any idea?

Hi, thanks @horasio for the feedback and @ocots for the tests. Two things (at least 🤞🏾):

• there is a known issue with constants and (seemingly) AD, see Pb with zero in dyn rhs? CTDirect.jl#190 : we have not taken the time to delve into it but it is a good occasion (though there are simple workarounds, it is not acceptable)
• there seems to be a specific issue with this example (0 iteration, you must be kidding me 🥲)

@horasio norm (unsquared) is not differentiable, so it is always advisable to smooth it

@ocots I see, however, that changing to mysqrt is not enough; to be continued. have you tried to change the linear solver, replacing MUMPS by, e.g., MA27/57?

[PierreMartinon]

Hi everyone @horasio @jbcaillau @ocots !
I see JB already answered concerning the AD problem.
As for the restoration error, I strongly suspect this is due to the (lack of) initial guess in the solve call: default initial guess is constant 0.1, however for this problem this is not good since the altitude should always be greater than 1. Can you try

direct_sol = solve(ocp; grid_size=100, init=(state=[1.01,0,0,0,0,0,1],))

[horasio]

Thank you everyone,
Using the smoothed sqrt as suggested by @ocots I get the same EXIT: Restoration Failed!.
By adding an init=(state=x0,)), as suggested by @PierreMartinon , Ipopt progresses: 13 iterations and EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
Is there a way, in init=... to provide a more precise warm start ?
In my code for the Goddard problem without control-toolbox (JuMP+Ipopt) I do not use smoothed sqrt and still do not have divide-by-zero errors and I think it's thanks to the warm start, that moved the trajectory away from zero at all times.

[PierreMartinon]

Thank you everyone, Using the smoothed sqrt as suggested by @ocots I get the same EXIT: Restoration Failed!. By adding an init=(state=x0,)), as suggested by @PierreMartinon , Ipopt progresses: 13 iterations and EXIT: Converged to a point of local infeasibility. Problem may be infeasible. Is there a way, in init=... to provide a more precise warm start ? In my code for the Goddard problem without control-toolbox (JuMP+Ipopt) I do not use smoothed sqrt and still do not have divide-by-zero errors and I think it's thanks to the warm start, that moved the trajectory away from zero at all times.

Hi @horasio, yes, you can find more details on the initial guess here:
https://control-toolbox.org/OptimalControl.jl/stable/tutorial-initial-guess.html

[oameye]

A get a similar error message for

using OptimalControl, LinearAlgebra
using NLPModelsIpopt

T = 50  # Time horizon

# Define the vector field
f(u, v) = [u - u^3 - 10*u*v^2,  -(1 - u^2)*v]
f(x) = f(x...)

mysqrt(x) = sqrt(x + 1e-1)
@def action begin
    t ∈ [0, T], time
    x ∈ R², state
    u ∈ R², control
    x(0) == [-1, 0]    # Starting point (left well)
    x(T) == [1, 0]     # End point (right well)
    ẋ(t) == u(t)       # Path dynamics
    # ∫( sum((u(t) - f(x(t))).^2) ) → min  # Minimize deviation from deterministic flow
    ∫(mysqrt(dot(u(t),u(t))*dot(f(x(t)),f(x(t)))) - dot(u(t),f(x(t)))) → min=
end

# Linear interpolation for x₁
x1(t) = -(1 - t/T) + t/T
# Parabolic guess for x₂
x2(t) = 0.3(-x1(t)^2 + 1)
x(t) = [x1(t), x2(t)]
u(t) = f(x(t))
init = (state=x, control=u)

sol = solve(action; init=init, grid_size=50)

[ocots]

It seems it is because dot(f(x(t)),f(x(t))) is 0 at some points and there is a problem with automatic diff when dealing with constant!

[ocots]

@jbcaillau @PierreMartinon

[jbcaillau]

@horasio @ocots @PierreMartinon solved using zero(x) instead of 0 which allows overloading to propagate dual numbers when the AD backend is ForwardDiff.jl (note: one(x) also exists) 🙏🏽 @abavoil for the hint

NB. @horasio doesn't the convergence issue of the optimiser on this use case 🥲

Dear CT community,

Describe the bug If the F0() function for the dynamics just returns a vector of zeros, then it compiles. But if I copy the velocity state into the time-deriv (x[4] in the line below marked as "BUG HERE") then the code does not compile, apparently in some automatic differentiation part. To Reproduce Execute the following:

module Goddard3D

using OptimalControl
using NLPModelsIpopt
#using OrdinaryDiffEq  # to get the Flow function from OptimalControl
#using NonlinearSolve  # interface to NLE solvers
#using MINPACK         # NLE solver: use to solve the shooting equation
using Plots           # to plot the solution

t0 = 0      # initial time
r0 = [0.999949994, 0.0001, 0.01]      # initial altitude
v0 = [0,0,0]      # initial speed, v = 1 ==> 7910 m/s (orbital velocity at sea level)
m0 = 1      # initial mass
x0 = [ r0..., v0..., m0 ]
rf = [1.01,0,0] # r[1] = 1 ==> 6378 km

@def ocp begin # definition of the optimal control problem
    tf ∈ R, variable
    t ∈ [ t0, tf ], time
    x ∈ R⁷, state
    u ∈ R³, control

    r = x[1:3]
    v = x[4:6]
    m = x[7]

    x(t0) == x0
    r(tf) == rf,         (1)
    u(t)' * u(t) ≤ 1.0
    r(t)' * r(t) ≥ 1.0
    #0 ≤ v(t) ≤ vmax

    ẋ(t) == F0(x(t)) + F1(x(t), u(t))

    m(tf) → max

end;

# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 7

F0(x) = begin
    r = x[1:3]
    v = x[4:6]
    m = x[7]
    rnorm2 = r' * r
    rnorm = sqrt(rnorm2)
    rnormalized = (1/rnorm) * r
    vnorm = sqrt(v' * v)
    D = ( Cd * vnorm * exp(-β*(rnorm - 1)) ) * v # Drag force
    return [v..., (-D/m - (1/rnorm2)*rnormalized)..., zero(x)[1] ]
end

F1(x,u) = begin
    m = x[7]
    unorm = sqrt(u' * u)
    t = (Tmax / m) * u
    return [ zero(x)[1:3]; t...; -b*Tmax*unorm ]
end

# DIRECT METHOD

function run_direct()
    direct_sol = solve(ocp, grid_size=100)
    plt = plot(direct_sol, size=(900, 900))
end

# INDIRECT SHOOTING METHOD

end # module

Goddard3D.run_direct()

[jbcaillau]

@oameye closing in favour of #384