control-toolbox / CTDirect.jl

Direct transcription of an optimal control problem and resolution
MIT License
4 stars 5 forks source link

Performance test (optimisation) #61

Open jbcaillau opened 1 year ago

jbcaillau commented 1 year ago

@ocots @PierreMartinon @gergaud Some preliminary tests with last package versions (CTDirect 0.4.1)

JuMP [^1]

IMG_2353 IMG_2354

ADNLP [^2]

IMG_2355 IMG_2356

[^1]: with this file:

using JuMP, Ipopt, Plots, MINPACK, BenchmarkTools

#JuMP model, Ipopt solver
sys = Model(optimizer_with_attributes(Ipopt.Optimizer,"print_level"=> 5))
set_optimizer_attribute(sys,"tol",1e-8)
set_optimizer_attribute(sys,"constr_viol_tol",1e-6)
set_optimizer_attribute(sys,"max_iter",1000)

# Parameters
Cd = 310.0
Tmax = 3.5
β = 500.0
b = 2.0
N = 100
t0 = 0.0
r0 = 1.0
v0 = 0.0
vmax = 0.1
m0 = 1.0
x0 = [ r0, v0, m0 ]
mf = 0.6

# Variables (some state constraints have been added to ease convergence)
@variables(sys, begin
    0.0 ≤ Δt                      # time step (unknown as tf is free)
    r[1:N+1] ≥ r0                 # r
    0 ≤ v[1:N+1] ≤ vmax           # v
    mf ≤ m[1:N+1] ≤ m0            # m
    0.0 ≤ u[1:N+1] ≤ 1.0          # u
end)

# Objective
@objective(sys, Max, r[N+1])

# Boundary constraints
@constraints(sys, begin
    con_r0, r[1] == r0
    con_v0, v[1] == v0
    con_m0, m[1] == m0
end)

# Dynamics
@NLexpressions(sys, begin
    # D = Cd v^2 exp(-β(r-1))
    D[i = 1:N+1], Cd * v[i]^2 * exp(-β * (r[i] - 1.0))
    # r'= v
    dr[i = 1:N+1], v[i]
    # v' = (Tmax.u-D)/m - 1/r^2
    dv[i = 1:N+1], (Tmax*u[i]-D[i])/m[i] - 1/r[i]^2
    # m' = -b.Tmax.u
    dm[i = 1:N+1], -b*Tmax*u[i]
end)

# Crank-Nicolson scheme
@NLconstraints(sys, begin
    con_dr[i = 1:N], r[i+1] == r[i] + Δt * (dr[i] + dr[i+1])/2.0
    con_dv[i = 1:N], v[i+1] == v[i] + Δt * (dv[i] + dv[i+1])/2.0
    con_dm[i = 1:N], m[i+1] == m[i] + Δt * (dm[i] + dm[i+1])/2.0
end);

# Solves for the control and state
println("Solving...")
@btime status = optimize!(sys)

[^2]: with this file

ocots commented 1 year ago

We can see more nonzeros for jump.

jbcaillau commented 1 year ago

We can see more nonzeros for jump.

Yes, this is in favour of ADNLP and should provide some increased performance 🤞🏾

PierreMartinon commented 3 months ago

Ok, this one needs some updating.