OptimalControl.jl icon indicating copy to clipboard operation
OptimalControl.jl copied to clipboard

[Bug] Goddard 3D does not compile

Open horasio opened this issue 4 months ago • 5 comments

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

horasio avatar Oct 01 '24 15:10 horasio