Repository navigation
Benchmarks: indirect shooting on the Goddard problem — v0.8.25 vs v0.14.0-beta #348
ocots
started this conversation in
Show and tell
Replies: 5 comments
|
Stable version [deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d"
CTParser = "32681960-a1b1-40db-9bff-a1ca817385d1"
CTSolvers = "d3e8d392-8e4b-4d9b-8e92-d7d4e3650ef6"
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
[compat]
CTFlows = "0.8"#!/usr/bin/env julia
using Pkg
# ---------------------------------------------------------------------- #
# Add the project to the path
# ---------------------------------------------------------------------- #
Pkg.activate(@__DIR__)
# Pkg.develop(path=joinpath(@__DIR__, "..", "..")) # We use the stable CTFlows
# ---------------------------------------------------------------------- #
# Load packages
# ---------------------------------------------------------------------- #
using CTParser
using CTModels
using CTDirect
using CTSolvers
using NLPModelsIpopt
using CommonSolve
using Plots
using CTFlows
using OrdinaryDiffEq
using NonlinearSolve
using BenchmarkTools
# ---------------------------------------------------------------------- #
# Optimal control problem
# ---------------------------------------------------------------------- #
const t0 = 0 # initial time
const r0 = 1 # initial altitude
const v0 = 0 # initial speed
const m0 = 1 # initial mass
const vmax = 0.1 # maximal authorized speed
const mf = 0.6 # final mass to target
goddard = CTParser.@def begin # definition of the optimal control problem
tf ∈ R, variable
t ∈ [t0, tf], time
x = (r, v, m) ∈ R³, state
u ∈ R, control
x(t0) == [r0, v0, m0]
m(tf) == mf, (1)
0 ≤ u(t) ≤ 1
r(t) ≥ r0
0 ≤ v(t) ≤ vmax
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
-r(tf) → min
end
# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 2
F0(x) = begin
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
return [v, -D/m - 1/r^2, 0]
end
F1(x) = begin
r, v, m = x
return [0, Tmax/m, -b*Tmax]
end
# ---------------------------------------------------------------------- #
# Solve the problem
# ---------------------------------------------------------------------- #
init = CTModels.build_initial_guess(goddard, nothing)
discretizer = CTDirect.Collocation(grid_size=100, scheme=:trapeze)
docp = CTDirect.discretize(goddard, discretizer);
modeler = CTSolvers.ADNLP(backend=:optimized)
nlp = CTSolvers.nlp_model(docp, init, modeler)
solver = CTSolvers.Ipopt(print_level=5, tol=1e-8, mu_strategy="adaptive")
nlp_sol = CommonSolve.solve(nlp, solver; display=true)
direct_sol = CTSolvers.ocp_solution(docp, nlp_sol, modeler)
# ---------------------------------------------------------------------- #
# Plot the solution
# ---------------------------------------------------------------------- #
plt = plot(direct_sol, label="direct", size=(800, 800))
# ---------------------------------------------------------------------- #
# Structure of the solution
# ---------------------------------------------------------------------- #
t = CTModels.time_grid(direct_sol) # the time grid as a vector
x = CTModels.state(direct_sol) # the state as a function of time
u = CTModels.control(direct_sol) # the control as a function of time
p = CTModels.costate(direct_sol) # the costate as a function of time
H1 = CTFlows.Lift(F1) # H1(x, p) = p' * F1(x)
φ(t) = H1(x(t), p(t)) # switching function
g(x) = vmax - x[2] # state constraint v ≤ vmax
u_plot = plot(t, u, linewidth=2, label = "u(t)")
H1_plot = plot(t, φ, linewidth=2, label = "H₁(x(t), p(t))")
g_plot = plot(t, g ∘ x, linewidth=2, label = "g(x(t))")
plot(u_plot, H1_plot, g_plot, layout=(3,1), size=(600, 600))
# ---------------------------------------------------------------------- #
# Flows
# ---------------------------------------------------------------------- #
# Controls
const u0 = 0 # off control
const u1 = 1 # bang control
H0 = CTFlows.Lift(F0) # H0(x, p) = p' * F0(x)
H01 = CTFlows.@Lie {H0, H1}
H001 = CTFlows.@Lie {H0, H01}
H101 = CTFlows.@Lie {H1, H01}
us(x, p) = -H001(x, p) / H101(x, p) # singular control
⋅ = CTFlows.:(⋅)
ub(x) = -(F0⋅g)(x) / (F1⋅g)(x) # boundary control
μ(x, p) = H01(x, p) / (F1⋅g)(x) # multiplier associated to the state constraint g
# Flows
f0 = CTFlows.Flow(goddard, (x, p, tf) -> u0; abstol=1e-8, reltol=1e-8)
f1 = CTFlows.Flow(goddard, (x, p, tf) -> u1; abstol=1e-8, reltol=1e-8)
fs = CTFlows.Flow(goddard, (x, p, tf) -> us(x, p); abstol=1e-8, reltol=1e-8)
fb = CTFlows.Flow(goddard, (x, p, tf) -> ub(x), (x, u, tf) -> g(x), (x, p, tf) -> μ(x, p); abstol=1e-8, reltol=1e-8)
# ---------------------------------------------------------------------- #
# Shooting function
# ---------------------------------------------------------------------- #
x0 = [r0, v0, m0] # initial state
function shoot!(s, p0, t1, t2, t3, tf)
x1, p1 = f1(t0, x0, p0, t1)
x2, p2 = fs(t1, x1, p1, t2)
x3, p3 = fb(t2, x2, p2, t3)
xf, pf = f0(t3, x3, p3, tf)
s[1] = xf[3] - mf # final mass constraint
s[2:3] = pf[1:2] - [1, 0] # transversality conditions: r(tf) and v(tf) are free
# objective is -r(tf) → p_r(tf) = 1, and v(tf) free → p_v(tf) = 0
s[4] = H1(x1, p1) # H1 = H01 = 0
s[5] = H01(x1, p1) # at the entrance of the singular arc
s[6] = g(x2) # g = 0 when entering the boundary arc
s[7] = H0(xf, pf) # since tf is free
end
# ---------------------------------------------------------------------- #
# Initial guess
# ---------------------------------------------------------------------- #
η = 1e-3
t13 = t[ abs.(φ.(t)) .≤ η ]
t23 = t[ 0 .≤ (g ∘ x).(t) .≤ η ]
p0 = p(t0)
t1 = min(t13...)
t2 = min(t23...)
t3 = max(t23...)
tf = t[end]
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)
# Norm of the shooting function at solution
using LinearAlgebra: norm
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
# ---------------------------------------------------------------------- #
# Indirect shooting
# ---------------------------------------------------------------------- #
ξ_guess = [p0..., t1, t2, t3, tf] # initial guess
shoot!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
prob = NonlinearProblem(shoot!, ξ_guess)
# resolution of S(ξ) = 0
sciml_sol = solve(prob; show_trace=Val(true))
# we retrieve the costate solution together with the times
ξ = sciml_sol.u
p0, t1, t2, t3, tf = ξ[1:3], ξ[4:end]...
println("\nNonlinearSolve results:")
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)
# NonlinearSolve results:
# p0 = [3.945764658689699, 0.15039559623166285, 0.05371271293970429]
# t1 = 0.023509684041887746
# t2 = 0.05973738089985649
# t3 = 0.10157134842431054
# tf = 0.2020474405709961
# Norm of the shooting function at solution
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
@assert norm(s) < 1e-6 "Indirect shooting failed for NonlinearSolve"
# ---------------------------------------------------------------------- #
# Benchmark
# ---------------------------------------------------------------------- #
@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
# BenchmarkTools.Trial: 6 samples with 1 evaluation per sample.
# Range (min … max): 826.997 ms … 1.125 s ┊ GC (min … max): 17.46% … 39.65%
# Time (median): 946.809 ms ┊ GC (median): 17.81%
# Time (mean ± σ): 952.434 ms ± 103.767 ms ┊ GC (mean ± σ): 22.38% ± 8.87%
# ▁ ▁ █ ▁ ▁
# █▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
# 827 ms Histogram: frequency by time 1.13 s <
# Memory estimate: 3.31 GiB, allocs estimate: 2426563.
@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
# BenchmarkTools.Trial: 11 samples with 1 evaluation per sample.
# Range (min … max): 370.776 ms … 757.985 ms ┊ GC (min … max): 15.73% … 55.76%
# Time (median): 398.101 ms ┊ GC (median): 19.72%
# Time (mean ± σ): 463.804 ms ± 133.593 ms ┊ GC (mean ± σ): 30.63% ± 14.44%
# ▁█ █▁▁ ▁ ▁ ▁ ▁
# ██▁███▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁█ ▁
# 371 ms Histogram: frequency by time 758 ms <
# Memory estimate: 1.82 GiB, allocs estimate: 577198.
# ---------------------------------------------------------------------- #
# Integration from Hamiltonian
# ---------------------------------------------------------------------- #
# Pseudo-Hamiltonian
H(x, p, u) = (p' * F0(x) + u * p' * F1(x))
# Hamiltonians
Q0(x, p, tf) = H(x, p, u0)
Q1(x, p, tf) = H(x, p, u1)
Qs(x, p, tf) = H(x, p, us(x, p))
Qb(x, p, tf) = H(x, p, ub(x)) + μ(x,p)*g(x)
φ0 = CTFlows.Flow(CTFlows.Hamiltonian(Q0; variable=true); abstol=1e-8, reltol=1e-8)
φ1 = CTFlows.Flow(CTFlows.Hamiltonian(Q1; variable=true); abstol=1e-8, reltol=1e-8)
φs = CTFlows.Flow(CTFlows.Hamiltonian(Qs; variable=true); abstol=1e-8, reltol=1e-8)
φb = CTFlows.Flow(CTFlows.Hamiltonian(Qb; variable=true); abstol=1e-8, reltol=1e-8)
p0 = [3.945764658689699, 0.15039559623166285, 0.05371271293970429]
t1 = 0.023509684041887746
t2 = 0.05973738089985649
t3 = 0.10157134842431054
tf = 0.2020474405709961
φf = φ1 * (t1, φs) * (t2, φb) * (t3, φ0) # concatenation of the flows
sol = φf((t0, tf), x0, p0, tf) # compute the solution: state, costate, control...
t = sol.t
x = t -> sol(t)[1:3]
p = t -> sol(t)[4:6]
# Plot X and P in two columns, each component per row
x1_plot = plot(t, t -> x(t)[1], label="r(t)", xlabel="t", ylabel="r", linewidth=2)
x2_plot = plot(t, t -> x(t)[2], label="v(t)", xlabel="t", ylabel="v", linewidth=2)
x3_plot = plot(t, t -> x(t)[3], label="m(t)", xlabel="t", ylabel="m", linewidth=2)
p1_plot = plot(t, t -> p(t)[1], label="p_r(t)", xlabel="t", ylabel="p_r", linewidth=2)
p2_plot = plot(t, t -> p(t)[2], label="p_v(t)", xlabel="t", ylabel="p_v", linewidth=2)
p3_plot = plot(t, t -> p(t)[3], label="p_m(t)", xlabel="t", ylabel="p_m", linewidth=2)
plot(x1_plot, p1_plot, x2_plot, p2_plot, x3_plot, p3_plot, layout=(3,2), size=(800, 800))
# Benchmark
@benchmark φf((t0, tf), x0, p0, tf)
# BenchmarkTools.Trial: 138 samples with 1 evaluation per sample.
# Range (min … max): 33.492 ms … 41.657 ms ┊ GC (min … max): 13.17% … 29.92%
# Time (median): 36.273 ms ┊ GC (median): 19.43%
# Time (mean ± σ): 36.287 ms ± 1.324 ms ┊ GC (mean ± σ): 19.40% ± 2.32%
# ▃ ▂▃ ▂▆ ▅█▃ ▂▅
# ▄▇▁▅▅▄▄▁▁▅███████████████▅██▅▅▅▇▇▇█▅▇▄▁▇▄▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▄
# 33.5 ms Histogram: frequency by time 41 ms <
# Memory estimate: 139.22 MiB, allocs estimate: 198640.
# Shooting function
x0 = [r0, v0, m0] # initial state
function shoot_from_hamiltonian!(s, p0, t1, t2, t3, tf)
x1, p1 = φ1(t0, x0, p0, t1, tf)
x2, p2 = φs(t1, x1, p1, t2, tf)
x3, p3 = φb(t2, x2, p2, t3, tf)
xf, pf = φ0(t3, x3, p3, tf, tf)
s[1] = xf[3] - mf # final mass constraint
s[2:3] = pf[1:2] - [1, 0] # transversality conditions: r(tf) and v(tf) are free
# objective is -r(tf) → p_r(tf) = 1, and v(tf) free → p_v(tf) = 0
s[4] = H1(x1, p1) # H1 = H01 = 0
s[5] = H01(x1, p1) # at the entrance of the singular arc
s[6] = g(x2) # g = 0 when entering the boundary arc
s[7] = H0(xf, pf) # since tf is free
return nothing
end
p0 = [3.94, 0.15, 0.05]
t1 = 0.02
t2 = 0.05
t3 = 0.10
tf = 0.20
s = similar(p0, 7)
shoot_from_hamiltonian!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
#
ξ_guess = [p0..., t1, t2, t3, tf] # initial guess
shoot_from_hamiltonian!(s, ξ, λ) = shoot_from_hamiltonian!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
prob = NonlinearProblem(shoot_from_hamiltonian!, ξ_guess)
sciml_sol = solve(prob; show_trace=Val(true))
# benchamark
@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
# BenchmarkTools.Trial: 8 samples with 1 evaluation per sample.
# Range (min … max): 626.525 ms … 650.121 ms ┊ GC (min … max): 21.41% … 22.18%
# Time (median): 628.933 ms ┊ GC (median): 21.40%
# Time (mean ± σ): 634.982 ms ± 9.979 ms ┊ GC (mean ± σ): 21.51% ± 0.50%
# ▁ ▁▁ █ ▁ ▁ ▁
# █▁██▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁█ ▁
# 627 ms Histogram: frequency by time 650 ms <
# Memory estimate: 2.82 GiB, allocs estimate: 2151575.
@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
# BenchmarkTools.Trial: 13 samples with 1 evaluation per sample.
# Range (min … max): 341.528 ms … 703.757 ms ┊ GC (min … max): 15.25% … 54.92%
# Time (median): 360.111 ms ┊ GC (median): 18.50%
# Time (mean ± σ): 409.960 ms ± 122.170 ms ┊ GC (mean ± σ): 28.21% ± 13.99%
# ▃▃▃█
# ████▇▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▇ ▁
# 342 ms Histogram: frequency by time 704 ms <
# Memory estimate: 1.65 GiB, allocs estimate: 525496. |
0 replies
|
Latest version with [deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
CTLie = "6880e05b-3a7d-4cac-887c-30cb52c5fdde"
CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"#!/usr/bin/env julia
#
# For the moment, we cannot construct any flow from an optimal control proble
# We thus construct from Hamiltonians
#
using Revise
using Pkg
# ---------------------------------------------------------------------- #
# Add the project to the path
# ---------------------------------------------------------------------- #
Pkg.activate(@__DIR__)
Pkg.develop(path=joinpath(@__DIR__, "..", "..")) # We use the stable CTFlows
# ---------------------------------------------------------------------- #
# Load packages
# ---------------------------------------------------------------------- #
using CTFlows
using OrdinaryDiffEqTsit5
using Plots
using NonlinearSolve
using BenchmarkTools
using LinearAlgebra
using CTModels
using CTLie
using CTBase
# ---------------------------------------------------------------------- #
# Optimal control problem data
# ---------------------------------------------------------------------- #
const t0 = 0. # initial time
const r0 = 1. # initial altitude
const v0 = 0. # initial speed
const m0 = 1. # initial mass
const vmax = 0.1 # maximal authorized speed
const mf = 0.6 # final mass to target
const x0 = [r0, v0, m0] # initial state
# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 2
F0(x) = begin
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
return [v, -D/m - 1/r^2, 0]
end
F1(x) = begin
r, v, m = x
return [0, Tmax/m, -b*Tmax]
end
# Goddard OCP: ẋ = F0(x) + u F1(x), free tf (variable), path constraint v ≤ vmax.
# Mayer-only objective (irrelevant to the pseudo-Hamiltonian flow: no Lagrange term).
function _build_goddard_ocp()
pre = CTModels.Building.PreModel()
CTModels.Building.time_dependence!(pre; autonomous=true)
CTModels.Building.state!(pre, 3)
CTModels.Building.control!(pre, 1)
CTModels.Building.variable!(pre, 1)
CTModels.Building.time!(pre; t0=t0, indf=1) # free final time = variable[1]
CTModels.Building.dynamics!(
pre, (r, t, x, u, v) -> (r .= F0(x) .+ u[1] .* F1(x); nothing)
)
CTModels.Building.objective!(pre, :min; mayer=(x0, xf, v) -> -xf[3])
CTModels.Building.constraint!(
pre,
:path;
f=(r, t, x, u, v) -> (r[1] = x[2]; nothing),
lb=[-Inf],
ub=[vmax],
label=:vmax,
)
return CTModels.Building.build(pre)
end
goddard = _build_goddard_ocp()
# Controls
const u0 = 0 # off control
const u1 = 1 # bang control
H0 = CTLie.Lift(F0) # H0(x, p) = p' * F0(x)
H1 = CTLie.Lift(F1) # H1(x, p) = p' * F1(x)
H01 = CTLie.@Lie {H0, H1}
H001 = CTLie.@Lie {H0, H01}
H101 = CTLie.@Lie {H1, H01}
us(x, p) = -H001(x, p) / H101(x, p) # singular control
g(x) = vmax - x[2] # state constraint v ≤ vmax
ub(x) = -CTLie.ad(F0, g)(x) / CTLie.ad(F1, g)(x) # boundary control
μ(x, p) = H01(x, p) / (CTLie.ad(F1, g)(x)) # multiplier associated to the state constraint g
# # Pseudo-Hamiltonian
# H(x, p, u) = (p' * F0(x) + u * p' * F1(x))
# # Hamiltonians
# Q0(x, p, tf) = H(x, p, u0)
# Q1(x, p, tf) = H(x, p, u1)
# Qs(x, p, tf) = H(x, p, us(x, p))
# Qb(x, p, tf) = H(x, p, ub(x)) + μ(x,p)*g(x)
# φ0 = CTFlows.Flow(CTFlows.Hamiltonian(Q0; is_variable=true)) #; prepare_cache=false)
# φ1 = CTFlows.Flow(CTFlows.Hamiltonian(Q1; is_variable=true)) #; prepare_cache=false)
# φs = CTFlows.Flow(CTFlows.Hamiltonian(Qs; is_variable=true)) #; prepare_cache=false)
# φb = CTFlows.Flow(CTFlows.Hamiltonian(Qb; is_variable=true)) #; prepare_cache=false)
φ0 = CTFlows.Flow(goddard, (x, p, tf) -> u0; abstol=1e-8, reltol=1e-8)
φ1 = CTFlows.Flow(goddard, (x, p, tf) -> u1; abstol=1e-8, reltol=1e-8)
φs = CTFlows.Flow(goddard, (x, p, tf) -> us(x, p); abstol=1e-8, reltol=1e-8)
φb = CTFlows.Flow(goddard, (x, p, tf) -> ub(x); constraint=(x, u, tf) -> g(x), multiplier=(x, p, tf) -> μ(x, p), abstol=1e-8, reltol=1e-8)
p0 = [3.945764658689699, 0.15039559623166285, 0.05371271293970429]
t1 = 0.023509684041887746
t2 = 0.05973738089985649
t3 = 0.10157134842431054
tf = 0.2020474405709961
φf = φ1 * (t1, φs) * (t2, φb) * (t3, φ0) # concatenation of the flows
sol = φf((t0, tf), x0, p0; variable=tf) # compute the solution: state, costate, control...
# Benchmark
@benchmark φf((t0, tf), x0, p0; variable=tf)
##############################################################################
# BenchmarkTools.Trial: 1284 samples with 1 evaluation per sample.
# Range (min … max): 3.270 ms … 14.573 ms ┊ GC (min … max): 0.00% … 74.79%
# Time (median): 3.516 ms ┊ GC (median): 0.00%
# Time (mean ± σ): 3.896 ms ± 906.767 μs ┊ GC (mean ± σ): 9.43% ± 14.09%
# ▅█▅
# ▂▃▄███▇▄▃▃▂▂▂▁▂▂▁▂▁▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▃▃▃▃▃▃▃▃▃▂▂▂▂▂ ▃
# 3.27 ms Histogram: frequency by time 6.12 ms <
# Memory estimate: 8.71 MiB, allocs estimate: 74462.
##############################################################################
t = CTModels.time_grid(sol)
x = CTModels.state(sol) # callable : x(0.5) → vecteur état
p = CTModels.costate(sol) # callable : p(0.5) → vecteur co-état
# # Plot X and P in two columns, each component per row
# x1_plot = plot(t, t -> x(t)[1], label="r(t)", xlabel="t", ylabel="r", linewidth=2)
# x2_plot = plot(t, t -> x(t)[2], label="v(t)", xlabel="t", ylabel="v", linewidth=2)
# x3_plot = plot(t, t -> x(t)[3], label="m(t)", xlabel="t", ylabel="m", linewidth=2)
# p1_plot = plot(t, t -> p(t)[1], label="p_r(t)", xlabel="t", ylabel="p_r", linewidth=2)
# p2_plot = plot(t, t -> p(t)[2], label="p_v(t)", xlabel="t", ylabel="p_v", linewidth=2)
# p3_plot = plot(t, t -> p(t)[3], label="p_m(t)", xlabel="t", ylabel="p_m", linewidth=2)
# plot(x1_plot, p1_plot, x2_plot, p2_plot, x3_plot, p3_plot, layout=(3,2), size=(800, 800))
# Shooting function
function shoot!(s, p0, t1, t2, t3, tf)
x1, p1 = φ1(t0, x0, p0, t1; variable=tf)
x2, p2 = φs(t1, x1, p1, t2; variable=tf)
x3, p3 = φb(t2, x2, p2, t3; variable=tf)
xf, pf = φ0(t3, x3, p3, tf; variable=tf)
s[1] = xf[3] - mf # final mass constraint
s[2:3] = pf[1:2] - [1, 0] # transversality conditions: r(tf) and v(tf) are free
# objective is -r(tf) → p_r(tf) = 1, and v(tf) free → p_v(tf) = 0
s[4] = H1(x1, p1) # H1 = H01 = 0
s[5] = H01(x1, p1) # at the entrance of the singular arc
s[6] = g(x2) # g = 0 when entering the boundary arc
s[7] = H0(xf, pf) # since tf is free
return nothing
end
p0 = [3.94, 0.15, 0.05]
t1 = 0.02
t2 = 0.05
t3 = 0.10
tf = 0.20
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
#
ξ_guess = [p0..., t1, t2, t3, tf] # initial guess
shoot!(s, ξ, _) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
prob = NonlinearProblem(shoot!, ξ_guess)
sciml_sol = solve(prob; show_trace=Val(true))
# benchmark
@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
##############################################################################
# BenchmarkTools.Trial: 19 samples with 1 evaluation per sample.
# Range (min … max): 264.689 ms … 281.280 ms ┊ GC (min … max): 11.55% … 11.98%
# Time (median): 272.502 ms ┊ GC (median): 12.02%
# Time (mean ± σ): 272.345 ms ± 3.862 ms ┊ GC (mean ± σ): 11.90% ± 0.49%
# ▁ ▁ ▁ ▁ ▁█ █▁█▁ ▁█ ▁ ▁ ▁
# █▁█▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁█▁▁██▁████▁▁▁▁▁██▁█▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
# 265 ms Histogram: frequency by time 281 ms <
# Memory estimate: 638.54 MiB, allocs estimate: 3447173.
##############################################################################
@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
##############################################################################
# BenchmarkTools.Trial: 38 samples with 1 evaluation per sample.
# Range (min … max): 129.085 ms … 137.948 ms ┊ GC (min … max): 15.73% … 18.54%
# Time (median): 133.098 ms ┊ GC (median): 16.88%
# Time (mean ± σ): 133.519 ms ± 1.868 ms ┊ GC (mean ± σ): 17.12% ± 0.75%
# ▃ ▃ █▃ ▃ █ ▃▃ ▃ ▃
# ▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▇██▇█▇█▁▇▇▇▁██▇▁█▇▇▇▁▇▁▁▁▁▁▁▁▁▁▇▇█▁▁▇▁▁▁▁▇ ▁
# 129 ms Histogram: frequency by time 138 ms <
# Memory estimate: 348.05 MiB, allocs estimate: 873016.
############################################################################## |
0 replies
|
Latest with #!/usr/bin/env julia
#
# Benchmark for the Goddard problem using the :partial Hamiltonian vector field mode
# (AD takes partial derivatives of the pseudo-Hamiltonian at a fixed control value).
#
# For the moment, we cannot construct any flow from an optimal control proble
# We thus construct from Hamiltonians
#
using Revise
using Pkg
# ---------------------------------------------------------------------- #
# Add the project to the path
# ---------------------------------------------------------------------- #
Pkg.activate(@__DIR__)
Pkg.develop(path=joinpath(@__DIR__, "..", "..")) # We use the stable CTFlows
# ---------------------------------------------------------------------- #
# Load packages
# ---------------------------------------------------------------------- #
using CTFlows
using OrdinaryDiffEqTsit5
using Plots
using NonlinearSolve
using BenchmarkTools
using LinearAlgebra
using CTModels
using CTLie
using CTBase
# ---------------------------------------------------------------------- #
# Optimal control problem data
# ---------------------------------------------------------------------- #
const t0 = 0. # initial time
const r0 = 1. # initial altitude
const v0 = 0. # initial speed
const m0 = 1. # initial mass
const vmax = 0.1 # maximal authorized speed
const mf = 0.6 # final mass to target
const x0 = [r0, v0, m0] # initial state
# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 2
F0(x) = begin
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
return [v, -D/m - 1/r^2, 0]
end
F1(x) = begin
r, v, m = x
return [0, Tmax/m, -b*Tmax]
end
# Goddard OCP: ẋ = F0(x) + u F1(x), free tf (variable), path constraint v ≤ vmax.
# Mayer-only objective (irrelevant to the pseudo-Hamiltonian flow: no Lagrange term).
function _build_goddard_ocp()
pre = CTModels.Building.PreModel()
CTModels.Building.time_dependence!(pre; autonomous=true)
CTModels.Building.state!(pre, 3)
CTModels.Building.control!(pre, 1)
CTModels.Building.variable!(pre, 1)
CTModels.Building.time!(pre; t0=t0, indf=1) # free final time = variable[1]
CTModels.Building.dynamics!(
pre, (r, t, x, u, v) -> (r .= F0(x) .+ u[1] .* F1(x); nothing)
)
CTModels.Building.objective!(pre, :min; mayer=(x0, xf, v) -> -xf[3])
CTModels.Building.constraint!(
pre,
:path;
f=(r, t, x, u, v) -> (r[1] = x[2]; nothing),
lb=[-Inf],
ub=[vmax],
label=:vmax,
)
return CTModels.Building.build(pre)
end
goddard = _build_goddard_ocp()
# Controls
const u0 = 0 # off control
const u1 = 1 # bang control
H0 = CTLie.Lift(F0) # H0(x, p) = p' * F0(x)
H1 = CTLie.Lift(F1) # H1(x, p) = p' * F1(x)
H01 = CTLie.@Lie {H0, H1}
H001 = CTLie.@Lie {H0, H01}
H101 = CTLie.@Lie {H1, H01}
us(x, p) = -H001(x, p) / H101(x, p) # singular control
g(x) = vmax - x[2] # state constraint v ≤ vmax
ub(x) = -CTLie.ad(F0, g)(x) / CTLie.ad(F1, g)(x) # boundary control
μ(x, p) = H01(x, p) / (CTLie.ad(F1, g)(x)) # multiplier associated to the state constraint g
φ0 = CTFlows.Flow(goddard, (x, p, tf) -> u0; hamiltonian_type=:partial, abstol=1e-8, reltol=1e-8)
φ1 = CTFlows.Flow(goddard, (x, p, tf) -> u1; hamiltonian_type=:partial, abstol=1e-8, reltol=1e-8)
φs = CTFlows.Flow(goddard, (x, p, tf) -> us(x, p); hamiltonian_type=:partial, abstol=1e-8, reltol=1e-8)
φb = CTFlows.Flow(goddard, (x, p, tf) -> ub(x); constraint=(x, u, tf) -> g(x), multiplier=(x, p, tf) -> μ(x, p), hamiltonian_type=:partial, abstol=1e-8, reltol=1e-8)
p0 = [3.945764658689699, 0.15039559623166285, 0.05371271293970429]
t1 = 0.023509684041887746
t2 = 0.05973738089985649
t3 = 0.10157134842431054
tf = 0.2020474405709961
φf = φ1 * (t1, φs) * (t2, φb) * (t3, φ0) # concatenation of the flows
sol = φf((t0, tf), x0, p0; variable=tf) # compute the solution: state, costate, control...
# Benchmark
@benchmark φf((t0, tf), x0, p0; variable=tf)
##############################################################################
# BenchmarkTools.Trial: 2879 samples with 1 evaluation per sample.
# Range (min … max): 1.418 ms … 9.482 ms ┊ GC (min … max): 0.00% … 81.94%
# Time (median): 1.550 ms ┊ GC (median): 0.00%
# Time (mean ± σ): 1.737 ms ± 825.869 μs ┊ GC (mean ± σ): 10.17% ± 14.66%
# ▁█▇▃ ▁
# █████▅▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▆███▇ █
# 1.42 ms Histogram: log(frequency) by time 5.49 ms <
# Memory estimate: 2.85 MiB, allocs estimate: 44239.
##############################################################################
t = CTModels.time_grid(sol)
x = CTModels.state(sol) # callable : x(0.5) → vecteur état
p = CTModels.costate(sol) # callable : p(0.5) → vecteur co-état
# # Plot X and P in two columns, each component per row
# x1_plot = plot(t, t -> x(t)[1], label="r(t)", xlabel="t", ylabel="r", linewidth=2)
# x2_plot = plot(t, t -> x(t)[2], label="v(t)", xlabel="t", ylabel="v", linewidth=2)
# x3_plot = plot(t, t -> x(t)[3], label="m(t)", xlabel="t", ylabel="m", linewidth=2)
# p1_plot = plot(t, t -> p(t)[1], label="p_r(t)", xlabel="t", ylabel="p_r", linewidth=2)
# p2_plot = plot(t, t -> p(t)[2], label="p_v(t)", xlabel="t", ylabel="p_v", linewidth=2)
# p3_plot = plot(t, t -> p(t)[3], label="p_m(t)", xlabel="t", ylabel="p_m", linewidth=2)
# plot(x1_plot, p1_plot, x2_plot, p2_plot, x3_plot, p3_plot, layout=(3,2), size=(800, 800))
# Shooting function
function shoot!(s, p0, t1, t2, t3, tf)
x1, p1 = φ1(t0, x0, p0, t1; variable=tf)
x2, p2 = φs(t1, x1, p1, t2; variable=tf)
x3, p3 = φb(t2, x2, p2, t3; variable=tf)
xf, pf = φ0(t3, x3, p3, tf; variable=tf)
s[1] = xf[3] - mf # final mass constraint
s[2:3] = pf[1:2] - [1, 0] # transversality conditions: r(tf) and v(tf) are free
# objective is -r(tf) → p_r(tf) = 1, and v(tf) free → p_v(tf) = 0
s[4] = H1(x1, p1) # H1 = H01 = 0
s[5] = H01(x1, p1) # at the entrance of the singular arc
s[6] = g(x2) # g = 0 when entering the boundary arc
s[7] = H0(xf, pf) # since tf is free
return nothing
end
p0 = [3.94, 0.15, 0.05]
t1 = 0.02
t2 = 0.05
t3 = 0.10
tf = 0.20
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
#
ξ_guess = [p0..., t1, t2, t3, tf] # initial guess
shoot!(s, ξ, _) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
prob = NonlinearProblem(shoot!, ξ_guess)
sciml_sol = solve(prob; show_trace=Val(true))
# benchmark
@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
##############################################################################
# BenchmarkTools.Trial: 58 samples with 1 evaluation per sample.
# Range (min … max): 83.641 ms … 92.008 ms ┊ GC (min … max): 8.23% … 12.43%
# Time (median): 87.895 ms ┊ GC (median): 11.98%
# Time (mean ± σ): 86.911 ms ± 2.000 ms ┊ GC (mean ± σ): 10.62% ± 1.80%
# ▂▅ ▂█ ▂ ▅ ▂
# ▅▁▁▁▅▁█▁████▁▅█▁█▅▅█▁▁▅▁▅▁▁▁▁▁▁▅▁▁▁▁▁▁▁▁▅▁▅▁██▅▁█▁███▅██▁▁█ ▁
# 83.6 ms Histogram: frequency by time 89.4 ms <
# Memory estimate: 181.40 MiB, allocs estimate: 2013639.
##############################################################################
@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve
##############################################################################
# BenchmarkTools.Trial: 162 samples with 1 evaluation per sample.
# Range (min … max): 28.411 ms … 35.570 ms ┊ GC (min … max): 8.61% … 18.16%
# Time (median): 30.245 ms ┊ GC (median): 10.68%
# Time (mean ± σ): 30.910 ms ± 1.627 ms ┊ GC (mean ± σ): 14.02% ± 4.08%
# ▂▂ ▇▃▃█ ▂ ▂ ▃ ▅ ▃ ▂
# ▅▁▃▁▁▆▃██▇███████▇▁▅▃▁▁▁▃▃▁▁▁▁▅▁▁▁▁▁▁▁▃▅▃▆▁▆█▁███▆█████▇▃▁▅ ▃
# 28.4 ms Histogram: frequency by time 33.2 ms <
# Memory estimate: 86.16 MiB, allocs estimate: 499273.
############################################################################## |
0 replies
|
Prompt for analysis: # Deep analysis of CTFlows performance improvements
## Context
I have two versions of the Julia package **CTFlows.jl**.
- Stable version:
https://github.com/control-toolbox/CTFlows.jl/tree/v0.8.25
- Latest beta version:
https://github.com/control-toolbox/CTFlows.jl/tree/v0.14.0-beta
Benchmark results comparing both versions are available here:
https://github.com/control-toolbox/CTFlows.jl/discussions/348
The latest version is significantly faster and allocates much less memory.
I do **not** want a superficial explanation based only on reading the code.
I want a **deep forensic performance analysis** that explains, demonstrates and quantifies **where the gains come from**.
---
# Objective
Produce a rigorous technical report identifying every significant source of performance improvement between the two versions.
The report should rely on:
- source code inspection,
- compiler reasoning,
- Julia type inference,
- allocation analysis,
- generated LLVM/native code when relevant,
- and especially **minimal reproducible examples** that isolate each optimization.
The objective is to distinguish improvements that are:
- algorithmic,
- compiler-related,
- type-stability improvements,
- dispatch improvements,
- allocation reductions,
- cache friendliness,
- better data structures,
- inlining,
- constant propagation,
- specialization,
- escape analysis,
- etc.
---
# Important requirement
Whenever you identify a possible optimization,
**do not simply state it.**
Instead:
1. isolate the mechanism,
2. write a minimal Julia example reproducing it,
3. benchmark it,
4. explain why the compiler behaves differently,
5. quantify the impact.
For example, if replacing a mutable struct by an immutable struct improves performance, create a tiny self-contained example demonstrating exactly that effect.
The goal is to build a collection of "performance experiments" explaining the evolution of CTFlows.
---
# Required methodology
For every hypothesis:
## 1. Locate the code
Identify:
- commit(s),
- file(s),
- functions,
- changed data structures,
- API changes,
- internal implementation differences.
Explain why those modifications might influence performance.
---
## 2. Build a minimal reproduction
Construct the smallest possible Julia program reproducing the optimization.
For instance:
- replacing vectors by tuples
- concrete vs abstract fields
- Val dispatch
- StaticArrays
- function barriers
- removing Any
- avoiding splatting
- preallocation
- generated functions
- parametric structs
- immutable objects
- iterator redesign
- replacing Dict by tuple lookup
- avoiding closures
- removing captured variables
- etc.
The examples should be independent from CTFlows whenever possible.
---
## 3. Measure
Benchmark each example using BenchmarkTools.
Measure:
- runtime
- allocations
- allocated memory
Whenever useful:
- @code_warntype
- @code_typed
- @code_llvm
- @code_native
Explain what changes in generated code.
---
## 4. Explain
For every optimization explain:
- why allocations disappear
- why inference succeeds
- why dynamic dispatch disappears
- why LLVM optimizes better
- why SIMD/vectorization becomes possible
- why cache locality improves
- etc.
---
## 5. Estimate contribution
Estimate whether the optimization is likely
- negligible
- moderate
- major
- dominant
with respect to the benchmark improvements.
---
# Explore every possible source
Do not assume the gains come from one change.
Investigate all plausible causes, including:
## Type inference
- abstract fields
- unions
- Any
- parametric methods
- specialization
## Allocations
- temporary arrays
- tuples
- views
- slices
- comprehensions
- generators
- closures
- captured variables
## Dispatch
- dynamic dispatch
- multiple dispatch
- function barriers
## Memory layout
- immutable structs
- mutable structs
- StaticArrays
- tuples
- NamedTuples
- structs of arrays
- arrays of structs
## Algorithms
- asymptotic improvements
- fewer passes
- removed work
- caching
## Compiler friendliness
- inlining
- constant propagation
- dead code elimination
- escape analysis
- scalar replacement
- loop unrolling
## Data structures
- Dict
- Vector
- Tuple
- NTuple
- Set
- BitSet
- BitVector
## API redesign
Investigate whether API changes enabled stronger compiler specialization.
---
# Compare call graphs
When appropriate:
- inspect call graphs
- identify removed layers
- identify removed wrappers
- identify simplified abstractions
---
# Compare inference quality
Use
```
@code_warntype
```
to compare important functions between both versions.
Explain inference improvements.
---
# Compare allocations
Use
```
@allocated
```
and BenchmarkTools.
Identify where allocations disappeared.
---
# Compare generated code
Whenever interesting, compare
```
@code_typed
@code_llvm
@code_native
```
between both versions.
Explain how compiler output changed.
---
# Validate hypotheses
Do **not** speculate.
Every claim should be validated by either
- code inspection,
- compiler output,
- benchmark,
- or minimal reproduction.
If something cannot be demonstrated, explicitly state the uncertainty.
---
# Deliverables
Produce:
## 1. Executive summary
One page summarizing the main causes of the improvements.
---
## 2. Ranked list
Rank every identified optimization by estimated impact.
For example
| Optimization | Estimated impact | Confidence |
|--------------|-----------------|------------|
| Better type inference | High | High |
| Removed temporary allocations | Very High | High |
| Better specialization | Moderate | Medium |
| ... | ... | ... |
---
## 3. Detailed report
For every optimization include
- code location
- explanation
- minimal reproduction
- benchmark
- compiler analysis
- conclusion
---
## 4. Appendix
Include every minimal benchmark developed during the analysis.
Each benchmark should be executable independently.
---
# Goal
The final report should explain not only **what changed**, but **why Julia became faster**, in a way that could serve as documentation or even as material for a conference talk on Julia performance engineering. |
0 replies
!!! Caveat: the analysis below is with tol = 1e-10 for stable vs. 1e-8 for beta. Not consistant with the updated tests above that use tol = 1e-8 everywhereForensic performance analysis — part 1: reproduction & ranked causesAnalysis performed with Claude Code (Fable 5). Environment: Julia 1.12.6, 16× Intel Xeon Gold 5218 @ 2.3 GHz (cascadelake). Stage 1–2 of the analysis requested in the prompt above: setup/reproduction, source diff, and instrumentation. Minimal reproductions per mechanism (stage 3) to follow. Stage 1 — ReproductionTwo isolated environments: Reproduction quality: allocation counts match the posted tables almost exactly (e.g. stable from-Hamiltonian NewtonRaphson: 4,794,917 allocs — identical; beta
Stage 2 — Where the gains actually come fromThe accountingPer-arc integration at identical tolerance (1e-8), arcs solved separately (nf = RHS evaluations):
→ arc-level stepping is essentially identical across versions (±13%). The large step-count difference in the benchmarks comes from elsewhere:
Per-RHS-call cost (μs, BenchmarkTools minimum, representative point):
Cross-check: Σ(nf × cost/eval) reproduces the measured whole-solve times within ~20–30% (beta Ranked causes
Honest caveats
Interaction structuretime ≈ nf × cost/eval + GC(allocations). Items 1 & 4 act on nf; 2, 5–8 on cost/eval; 3 on allocations (hence GC). The shooting benchmarks multiply everything by the Newton loop with 7-partial Duals — which is why the memory ratios (43–50×) exceed the time ratios (20–33×). Next: stage 3 — one standalone minimal reproduction + benchmark per ranked item (switching-RHS vs sequential solves; gradient-through-control vs frozen-control partials; dense vs endpoint-only saving with Dual eltypes; Union-returning index helper inside an AD closure; 2n-dim vs split n-dim gradients; etc.). |
0 replies
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
This post documents the performance of CTFlows on the classic Goddard rocket problem (3-state, free final time, path constraint, 4-arc structure) solved by indirect shooting. All benchmarks were run with BenchmarkTools.jl on the same machine.
Versions compared
Problem recap
The optimal control problem and its structure are fully described in the tutorial. The indirect approach concatenates four Hamiltonian flows (bang–singular–boundary–bang) and solves the shooting equation
S(ξ) = 0(7 equations in 7 unknowns) via NonlinearSolve.jl.Flow construction: what changed
All flows are constructed with
abstol=1e-8, reltol=1e-8so that ODE tolerances are identical across all setups.In v0.8.25 (stable), flows are built from the OCP and a feedback control law:
or from a Hamiltonian:
In v0.14.0-beta (latest), the API remains almost the same:
A new keyword
hamiltonian_typecontrols how the Hamiltonian vector field is derived from the pseudo-HamiltonianH̃(x, p, u)::total(default) — substitutes the controlu(x, p, tf)first, then differentiates the full composition via AD.:partial— differentiatesH̃at a fixed control value (partial AD), then evaluates atu(x, p, tf). This skips AD through the control function, which can be expensive whenu(x, p)itself involves AD (e.g., singular control computed from Lie brackets).See the documentation for details.
Benchmark 1 — single flow evaluation
Raw results
Hamiltonian(stable baseline):total(default):partialGains
:totalvs stable:partialvs stable:partialvs:totalBenchmark 2 — indirect shooting, default solver (NewtonRaphson)
Raw results
Hamiltonian:total:partialGains (stable baseline = v0.8.25 from OCP)
:totalvs stable:partialvs stable:partialvs:totalBenchmark 3 — indirect shooting, SimpleNewtonRaphson
Raw results
Hamiltonian:total:partialGains (stable baseline = v0.8.25 from OCP)
:totalvs stable:partialvs stable:partialvs:totalSummary
The v0.14.0-beta is significantly faster and lighter than v0.8.25 at equal ODE tolerance (
abstol=reltol=1e-8). For a single flow evaluation, the:partialHamiltonian mode is 23× faster and uses 49× less memory than the stable baseline;:totalgives 10× faster and 16× less. For the full indirect shooting solve (SimpleNewtonRaphson),:partialis 13× faster and uses 22× less memory;:totalgives 3× faster and 5× less. Choosing:partialover:totaladds another 2–4× in time and 3–4× in memory at no cost in correctness.These results are specific to the Goddard problem; gains will vary depending on the complexity of the control law and the ODE solver settings.
Benchmark scripts are available below.
All reactions