OptimalControl.jl
OptimalControl.jl copied to clipboard
[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