control-toolbox
[Bug] Goddard 3D does not compile
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
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 ]
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!
@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 https://github.com/control-toolbox/CTDirect.jl/issues/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?
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],))
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.
Thank you everyone, Using the smoothed sqrt as suggested by @ocots I get the same
EXIT: Restoration Failed!. By adding aninit=(state=x0,)), as suggested by @PierreMartinon , Ipopt progresses: 13 iterations andEXIT: Converged to a point of local infeasibility. Problem may be infeasible.Is there a way, ininit=...to provide a more precise warm start ? In my code for the Goddard problem withoutcontrol-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
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)
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!
@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()
@oameye closing in favour of https://github.com/control-toolbox/OptimalControl.jl/issues/384