Skip to content

Commit

Permalink
Examples for initial guess options (#78)
Browse files Browse the repository at this point in the history
* Examples for constant initial guess (vector)

* Examples for warm start and passing initial guess at different levels

* Examples for functional initial guess and mixed functional/constant; bugfix in warm start generation with free final time

* updated CI test for initial guess options

* renamed solveDOCP and solveDirect to solve
  • Loading branch information
PierreMartinon authored Apr 17, 2024
1 parent b1dd97e commit 4e9fff4
Show file tree
Hide file tree
Showing 13 changed files with 359 additions and 156 deletions.
10 changes: 5 additions & 5 deletions src/CTDirect.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,18 +22,18 @@ include("solution.jl")
include("solve.jl")

# export functions only for user
export available_methods
export is_solvable
export solve
export OptimalControlInit
export directTranscription
export getNLP
export setDOCPInit
export solveDOCP
export solveDirect
export OCPSolutionFromDOCP
export initial_guess
export ipopt_objective
export ipopt_constraint
export available_methods
export DOCP_objective
export DOCP_constraints!


# CTBase reexports
export @def
Expand Down
6 changes: 3 additions & 3 deletions src/solve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ $(TYPEDSIGNATURES)
Solve a discretized optimal control problem
"""
function solveDOCP(docp::DOCP;
function solve(docp::DOCP;
init=nothing,
display::Bool=__display(),
print_level::Integer=__print_level_ipopt(),
Expand All @@ -76,7 +76,7 @@ function solveDOCP(docp::DOCP;
end


function solveDirect(ocp::OptimalControlModel,
function solve(ocp::OptimalControlModel,
description...;
init::OptimalControlInit=OptimalControlInit(),
grid_size::Integer=__grid_size_direct(),
Expand All @@ -88,7 +88,7 @@ function solveDirect(ocp::OptimalControlModel,
# build discretized OCP
docp = directTranscription(ocp, description, init=init, grid_size=grid_size)
# solve DOCP and retrieve OCP solution
ocp_solution = solveDOCP(docp; display=display, print_level=print_level, mu_strategy=mu_strategy, kwargs...)
ocp_solution = solve(docp; display=display, print_level=print_level, mu_strategy=mu_strategy, kwargs...)

return ocp_solution
end
13 changes: 7 additions & 6 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,13 @@ function initial_guess(docp, init::OptimalControlInit=OptimalControlInit())
# note: internal variables (lagrange cost, k_i for RK schemes) will keep these default values
xuv = 0.1 * ones(docp.dim_NLP_variables)

#init = docp.NLP_init
# set variables if provided
# (needed first in case of free times !)
if !isnothing(init.variable_init)
set_variable!(xuv, init.variable_init, docp)
end

# get time grid for state / control variables
N = docp.dim_NLP_steps
t0 = get_initial_time(xuv, docp)
tf = get_final_time(xuv, docp)
Expand All @@ -142,11 +148,6 @@ function initial_guess(docp, init::OptimalControlInit=OptimalControlInit())
if !isnothing(init.control_init(ti))
set_control_at_time_step!(xuv, init.control_init(ti), docp, i)
end

# set variables if provided
if !isnothing(init.variable_init)
set_variable!(xuv, init.variable_init, docp)
end
end

return xuv
Expand Down
6 changes: 3 additions & 3 deletions test/suite/abstract_ocp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ println("Test: abstract OCP definition")
end

@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin
sol1 = solveDirect(ocp1, grid_size=100, print_level=0, tol=1e-12)
sol1 = solve(ocp1, grid_size=100, print_level=0, tol=1e-12)
@test sol1.objective 2.0 rtol=1e-2
end

Expand All @@ -42,7 +42,7 @@ end
tf min
end

@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constraints" begin
sol2 = solveDirect(ocp2, grid_size=100, print_level=0, tol=1e-12)
@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constr" begin
sol2 = solve(ocp2, grid_size=100, print_level=0, tol=1e-12)
@test sol2.objective 5.46 rtol=1e-2
end
80 changes: 6 additions & 74 deletions test/suite/all_constraints_types.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
using CTDirect

println("Test: all constraint types")

# goddard max final altitue (all constraint types formulation)
# +++ todo: redo this one with different constraint formulations, use default initial guess only
# goddard max final altitude (all constraint types formulation)
ocp = Model(variable=true)
Cd = 310
Tmax = 3.5
Expand Down Expand Up @@ -39,86 +39,18 @@ constraint!(ocp, :control, Index(1), 0, Inf, :control_box_lb)
# variable box
constraint!(ocp, :variable, Index(1), 0.01, Inf, :variable_box_tfmin)
objective!(ocp, :mayer, (x0, xf, v) -> xf[1], :max)
function F0(x)
function FF0(x)
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1))
return [ v, -D/m - 1/r^2, 0 ]
end
function F1(x)
function FF1(x)
r, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
end
dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
dynamics!(ocp, (x, u, v) -> FF0(x) + u*FF1(x) )

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints" begin
sol1 = solveDirect(ocp, grid_size=100, print_level=0, tol=1e-8)
sol1 = solve(ocp, grid_size=100, print_level=0, tol=1e-8)
@test sol1.objective 1.0125 rtol=1e-2
end

# with constant initial guess (vector / function formats)
x_init = [1.05, 0.1, 0.8]
u_init = 0.5
v_init = 0.1

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant" begin
init_constant = OptimalControlInit(x_init=x_init, u_init=u_init, v_init=v_init)
sol2 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=init_constant)
@test sol2.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_function_x" begin
init_function_x = OptimalControlInit(x_init=t->x_init, u_init=u_init, v_init=v_init)
sol = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=init_function_x)
@test sol.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_function_u" begin
init_function_u = OptimalControlInit(x_init=x_init, u_init=t->u_init, v_init=v_init)
sol = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=init_function_u)
@test sol.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_function_xu" begin
init_function_xu = OptimalControlInit(x_init=t->x_init, u_init=t->u_init, v_init=v_init)
sol = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=init_function_xu)
@test sol.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (x)" begin
sol3 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(x_init=x_init))
@test sol3.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (u)" begin
sol4 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(u_init=u_init))
@test sol4.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (v)" begin
sol5 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(v_init=v_init))
@test sol5.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (x,u)" begin
sol6 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(x_init=x_init, u_init=u_init))
@test sol6.objective 1.0125 rtol=1e-2
end

@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (x,v)" begin
sol7 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(x_init=x_init, v_init=v_init))
@test sol7.objective 1.0125 rtol=1e-2
end


@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_constant (u,v)" begin
sol8 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=OptimalControlInit(u_init=u_init, v_init=v_init))
@test sol8.objective 1.0125 rtol=1e-2
end

# with initial guess from solution
sol = solveDirect(ocp, grid_size=100, print_level=0, tol=1e-8)
@testset verbose = true showtiming = true ":goddard :max_rf :all_constraints :init_sol" begin
init_sol = OptimalControlInit(sol)
sol9 = solveDirect(ocp, grid_size=30, print_level=0, tol=1e-8, init=init_sol)
@test sol9.objective 1.0125 rtol=1e-2
end
10 changes: 5 additions & 5 deletions test/suite/double_integrator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ dynamics!(ocp1, (x, u, v) -> [x[2], u])
objective!(ocp1, :mayer, (x0, xf, v) -> v)

@testset verbose = true showtiming = true ":double_integrator :min_tf" begin
sol1 = solveDirect(ocp1, grid_size=100, print_level=0, tol=1e-12)
sol1 = solve(ocp1, grid_size=100, print_level=0, tol=1e-12)
@test sol1.objective 2.0 rtol=1e-2
end

Expand All @@ -35,7 +35,7 @@ dynamics!(ocp2, (x, u, v) -> [x[2], u])
objective!(ocp2, :lagrange, (x, u, v) -> 1)

@testset verbose = true showtiming = true ":double_integrator :min_tf :lagrange" begin
sol2 = solveDirect(ocp2, grid_size=100, print_level=0, tol=1e-12)
sol2 = solve(ocp2, grid_size=100, print_level=0, tol=1e-12)
@test sol2.objective 2.0 rtol=1e-2
end

Expand All @@ -54,7 +54,7 @@ dynamics!(ocp3, (x, u, v) -> [x[2], u[1]])
objective!(ocp3, :mayer, (x0, xf, v) -> v[1])

@testset verbose = true showtiming = true ":double_integrator :min_tf :vectorial" begin
sol3 = solveDirect(ocp3, grid_size=100, print_level=0, tol=1e-12)
sol3 = solve(ocp3, grid_size=100, print_level=0, tol=1e-12)
@test sol3.objective 2.0 rtol=1e-2
end

Expand All @@ -73,7 +73,7 @@ dynamics!(ocp4, (x, u, v) -> [x[2], u])
objective!(ocp4, :mayer, (x0, xf, v) -> v[1], :max)

@testset verbose = true showtiming = true ":double_integrator :max_t0" begin
sol4 = solveDirect(ocp4, grid_size=100, print_level=0, tol=1e-12)
sol4 = solve(ocp4, grid_size=100, print_level=0, tol=1e-12)
@test sol4.objective 8.0 rtol=1e-2
end

Expand All @@ -89,6 +89,6 @@ dynamics!(ocp5, (x, u) -> [x[2], -u[1] + u[2]])
objective!(ocp5, :lagrange, (x, u) -> u[1]*u[1] + u[2]*u[2])

@testset verbose = true showtiming = true ":double_integrator :min_energy" begin
sol5 = solveDirect(ocp5, grid_size=50, print_level=0, tol=1e-12)
sol5 = solve(ocp5, grid_size=50, print_level=0, tol=1e-12)
@test sol5.objective 9.6e-2 rtol=1e-2
end
150 changes: 150 additions & 0 deletions test/suite/initial_guess.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
using CTDirect

# goddard max final altitude (all constraint types formulation)
println("Test: initial guess options")
ocp = Model(variable=true)
Cd = 310
Tmax = 3.5
β = 500
b = 2
r0 = 1
v0 = 0
vmax = 0.1
m0 = 1
mf = 0.6
x0 = [ r0, v0, m0 ]
state!(ocp, 3)
control!(ocp, 1)
variable!(ocp, 1)
time!(ocp, 0, Index(1))
# use all possible types of constraints
# initial condition
constraint!(ocp, :initial, x0, :initial_constraint)
# final condition
constraint!(ocp, :final, Index(3), mf, :final_constraint)
# state constraint
constraint!(ocp, :state, (x,v)->x[2], -Inf, vmax, :state_con_v_ub)
# control constraint
constraint!(ocp, :control, (u,v)->u, -Inf, 1, :control_con_u_ub)
# mixed constraint
constraint!(ocp, :mixed, (x,u,v)->x[3], mf, Inf, :mixed_con_m_lb)
# variable constraint
constraint!(ocp, :variable, v->v, -Inf, 10, :variable_con_tf_ubx)
# state box
constraint!(ocp, :state, 1:2, [r0,v0], [r0+0.2, Inf], :state_box_rv)
# control box
constraint!(ocp, :control, Index(1), 0, Inf, :control_box_lb)
# variable box
constraint!(ocp, :variable, Index(1), 0.01, Inf, :variable_box_tfmin)
objective!(ocp, :mayer, (x0, xf, v) -> xf[1], :max)
function F0(x)
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1))
return [ v, -D/m - 1/r^2, 0 ]
end
function F1(x)
r, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
end
dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
sol0 = solve(ocp, print_level=0)

# default init
@testset verbose = true showtiming = true ":default_init" begin
sol = solve(ocp, print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end

# constant initial guess
x_const = [1.05, 0.2, 0.8]
u_const = 0.5
v_init = 0.15

@testset verbose = true showtiming = true ":constant_init_x" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_u" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_v" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_init))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_xu" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_xv" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_init))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_uv" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_init))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":constant_init_xuv" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
@test sol.objective 1.0125 rtol=1e-2
end

# functional initial guess
x_func = t->[1+t^2, sqrt(t), 1-t]
u_func = t->(cos(t)+1)*0.5

@testset verbose = true showtiming = true ":functional_init_x" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_func))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":functional_init_u" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_func))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":functional_init_xu" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_func, u_init=u_func))
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":mixed_init" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_func, u_init=u_const))
@test sol.objective 1.0125 rtol=1e-2
end

# warm start
@testset verbose = true showtiming = true ":warm_start" begin
sol = solve(ocp, print_level=0, init=OptimalControlInit(sol0))
@test sol.objective 1.0125 rtol=1e-2
end

# set initial guess in DOCP
docp = directTranscription(ocp)
@testset verbose = true showtiming = true ":DOCPInit_constant" begin
setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
sol = solve(docp, print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":DOCPInit_mixed" begin
setDOCPInit(docp, OptimalControlInit(x_init=x_func, u_init=u_const))
sol = solve(docp, print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":DOCPInit_warm_start" begin
setDOCPInit(docp, OptimalControlInit(sol0))
sol = solve(docp, print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end

# pass initial guess to solve
setDOCPInit(docp, OptimalControlInit()) # reset init in docp
@testset verbose = true showtiming = true ":solve_constant_init" begin
sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":solve_mixed_init" begin
sol = solve(docp, init=OptimalControlInit(x_init=x_func, u_init=u_const), print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end
@testset verbose = true showtiming = true ":solve_warm_start" begin
sol = solve(docp, init=OptimalControlInit(sol0), print_level=0)
@test sol.objective 1.0125 rtol=1e-2
end
Loading

0 comments on commit 4e9fff4

Please sign in to comment.