diff --git a/Project.toml b/Project.toml index 6ee27c9..cbde643 100644 --- a/Project.toml +++ b/Project.toml @@ -23,7 +23,7 @@ CTSolveExt = ["NLPModelsIpopt", "HSL"] [compat] ADNLPModels = "0.8" -CTBase = "0.10" +CTBase = "0.11" DocStringExtensions = "0.9" HSL = "0.4" julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index 4a79b27..4b44440 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,6 +7,3 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" - -[compat] -julia = "1.9" diff --git a/docs/make.jl b/docs/make.jl index 743a213..d9f62a9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,6 +6,7 @@ using HSL using JLD2 using JSON3 using Plots +using CommonSolve DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true) DocMeta.setdocmeta!(CTDirect, :DocTestSetup, :(using CTDirect); recursive = true) diff --git a/docs/src/api.md b/docs/src/api.md index cad02f0..cb660de 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -4,17 +4,10 @@ ```@index Pages = ["api.md"] -Modules = [CTDirect] +Modules = [CTDirect, CommonSolve] Order = [:module, :constant, :type, :function, :macro] ``` -## Available methods - -```@example -using CTDirect -available_methods() -``` - ## Documentation ```@autodocs @@ -23,8 +16,6 @@ Order = [:module, :constant, :type, :function, :macro] Private = false ``` - - ```@docs solve save_OCP_solution diff --git a/ext/CTSolveExt.jl b/ext/CTSolveExt.jl index 62a7aa7..3dd77f3 100644 --- a/ext/CTSolveExt.jl +++ b/ext/CTSolveExt.jl @@ -42,7 +42,7 @@ function CommonSolve.solve(docp::DOCP; else # use given initial guess ocp = docp.ocp - x0 = CTDirect.DOCP_initial_guess(docp, _OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) + x0 = CTDirect.DOCP_initial_guess(docp, OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) docp_solution = ipopt(nlp, x0=x0, print_level=print_level, mu_strategy=mu_strategy, sb="yes", linear_solver=linear_solver; kwargs...) end diff --git a/src/CTDirect.jl b/src/CTDirect.jl index 255930e..842a24a 100644 --- a/src/CTDirect.jl +++ b/src/CTDirect.jl @@ -22,7 +22,7 @@ include("solution.jl") include("solve.jl") # re exports -export solve # CommonSolve +export solve # from CommonSolve, extended in CTSolveExt # exports export available_methods @@ -39,6 +39,4 @@ export export_OCP_solution export read_OCP_solution export OCP_Solution_discrete -export _OptimalControlInit #temp, later from CTBase - end \ No newline at end of file diff --git a/src/problem.jl b/src/problem.jl index 89af1cb..dc95dca 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -177,9 +177,9 @@ function variables_bounds(docp) # state box offset = 0 - if dim_state_box(ocp) > 0 + if dim_state_range(ocp) > 0 for i in 0:N - for j in 1:dim_state_box(ocp) + for j in 1:dim_state_range(ocp) indice = docp.state_box[2][j] l_var[offset+indice] = docp.state_box[1][j] u_var[offset+indice] = docp.state_box[3][j] @@ -190,9 +190,9 @@ function variables_bounds(docp) # control box offset = (N+1) * docp.dim_NLP_x - if dim_control_box(ocp) > 0 + if dim_control_range(ocp) > 0 for i in 0:N - for j in 1:dim_control_box(ocp) + for j in 1:dim_control_range(ocp) indice = docp.control_box[2][j] l_var[offset+indice] = docp.control_box[1][j] u_var[offset+indice] = docp.control_box[3][j] @@ -203,8 +203,8 @@ function variables_bounds(docp) # variable box offset = (N+1) * (docp.dim_NLP_x + docp.dim_NLP_u) - if dim_variable_box(ocp) > 0 - for j in 1:dim_variable_box(ocp) + if dim_variable_range(ocp) > 0 + for j in 1:dim_variable_range(ocp) indice = docp.variable_box[2][j] l_var[offset+indice] = docp.variable_box[1][j] u_var[offset+indice] = docp.variable_box[3][j] diff --git a/src/solution.jl b/src/solution.jl index a35edaa..233564f 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -128,20 +128,20 @@ function OCPSolutionFromDOCP_raw(docp, T, X, U, v, P; end # box constraints multipliers - if dim_state_box(ocp) > 0 + if dim_state_range(ocp) > 0 # remove additional state for lagrange cost mbox_x_l = ctinterpolate(T, matrix2vec(mult_state_box_lower[:,1:ocp.state_dimension], 1)) mbox_x_u = ctinterpolate(T, matrix2vec(mult_state_box_upper[:,1:ocp.state_dimension], 1)) sol.infos[:mult_state_box_lower] = t -> mbox_x_l(t) sol.infos[:mult_state_box_upper] = t -> mbox_x_u(t) end - if dim_control_box(ocp) > 0 + if dim_control_range(ocp) > 0 mbox_u_l = ctinterpolate(T, matrix2vec(mult_control_box_lower, 1)) mbox_u_u = ctinterpolate(T, matrix2vec(mult_control_box_upper, 1)) sol.infos[:mult_control_box_lower] = t -> mbox_u_l(t) sol.infos[:mult_control_box_upper] = t -> mbox_u_u(t) end - if dim_variable_box(ocp) > 0 + if dim_variable_range(ocp) > 0 sol.infos[:mult_variable_box_lower] = mult_variable_box_lower sol.infos[:mult_variable_box_upper] = mult_variable_box_upper end @@ -193,7 +193,7 @@ function parse_DOCP_solution(docp, solution, multipliers_constraints, multiplier mult_state_box_upper[i,:] = vget_state_at_time_step(mult_U, docp, i-1) mult_control_box_upper[i,:] = vget_control_at_time_step(mult_U, docp, i-1) end - if dim_variable_box(ocp) > 0 + if dim_variable_range(ocp) > 0 mult_variable_box_lower = get_variable(mult_L, docp) mult_variable_box_upper = get_variable(mult_U, docp) end diff --git a/src/solve.jl b/src/solve.jl index 45ac3c0..b9973f8 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -29,7 +29,7 @@ function directTranscription(ocp::OptimalControlModel, docp = DOCP(ocp, grid_size, time_grid) # set initial guess - x0 = DOCP_initial_guess(docp, _OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) + x0 = DOCP_initial_guess(docp, OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) # set bounds docp.var_l, docp.var_u = variables_bounds(docp) @@ -67,6 +67,6 @@ function setInitialGuess(docp::DOCP, init) nlp = getNLP(docp) ocp = docp.ocp - nlp.meta.x0 .= DOCP_initial_guess(docp, _OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) + nlp.meta.x0 .= DOCP_initial_guess(docp, OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) end diff --git a/src/utils.jl b/src/utils.jl index 2b6c9fe..ae605d6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -188,238 +188,6 @@ function set_variable!(xu, v_init, docp) end end -####################################################### -# local init MOVE TO CTBASE -####################################################### -""" -$(TYPEDSIGNATURES) - -Check if actual dimension is equal to target dimension, error otherwise -""" -function checkDim(actual_dim, target_dim) - if !isnothing(target_dim) && actual_dim != target_dim - error("Init dimension mismatch: got ",actual_dim," instead of ",target_dim ) - end -end - -""" -$(TYPEDSIGNATURES) - -Return true if argument is a vector of vectors -""" -function isaVectVect(data) - return (data isa Vector) && (data[1] isa ctVector) -end - -""" -$(TYPEDSIGNATURES) - -Convert matrix to vector of vectors (could be expanded) -""" -function formatData(data) - if data isa Matrix - return matrix2vec(data,1) - else - return data - end -end - -""" -$(TYPEDSIGNATURES) - -Convert matrix time-grid to vector -""" -function formatTimeGrid(time) - if isnothing(time) - return nothing - elseif time isa ctVector - return time - else - return vec(time) - end -end - -""" -$(TYPEDSIGNATURES) - -Build functional initialization: default case -""" -function buildFunctionalInit(data::Nothing, time, dim) - # fallback to method-dependent default initialization - return t-> nothing -end - -""" -$(TYPEDSIGNATURES) - -Build functional initialization: function case -""" -function buildFunctionalInit(data::Function, time, dim) - # functional initialization - checkDim(length(data(0)),dim) - return t -> data(t) -end - -""" -$(TYPEDSIGNATURES) - -Build functional initialization: constant / 1D interpolation -""" -function buildFunctionalInit(data::ctVector, time, dim) - if !isnothing(time) && (length(data) == length(time)) - # interpolation vs time, dim 1 case - itp = ctinterpolate(time, data) - return t -> itp(t) - else - # constant initialization - checkDim(length(data), dim) - return t -> data - end -end - -""" -$(TYPEDSIGNATURES) - -Build functional initialization: general interpolation case -""" -function buildFunctionalInit(data, time, dim) - if isaVectVect(data) - # interpolation vs time, general case - itp = ctinterpolate(time, data) - checkDim(length(itp(0)), dim) - return t -> itp(t) - else - error("Unrecognized initialization argument: ",typeof(data)) - end -end - -""" -$(TYPEDSIGNATURES) - -Build vector initialization: default / vector case -""" -function buildVectorInit(data, dim) - if isnothing(data) - return data - else - checkDim(length(data),dim) - return data - end -end - -""" -$(TYPEDSIGNATURES) - -Initial guess for OCP, contains -- functions of time for the state and control variables -- vector for optimization variables -Initialization data for each field can be left to default or: -- vector for optimization variables -- constant / vector / function for state and control -- existing solution ('warm start') for all fields - -# Constructors: - -- `_OptimalControlInit()`: default initialization -- `_OptimalControlInit(state, control, variable, time)`: constant vector, function handles and / or matrices / vectors interpolated along given time grid -- `_OptimalControlInit(sol)`: from existing solution - -# Examples - -```julia-repl -julia> init = _OptimalControlInit() -julia> init = _OptimalControlInit(state=[0.1, 0.2], control=0.3) -julia> init = _OptimalControlInit(state=[0.1, 0.2], control=0.3, variable=0.5) -julia> init = _OptimalControlInit(state=[0.1, 0.2], controlt=t->sin(t), variable=0.5) -julia> init = _OptimalControlInit(state=[[0, 0], [1, 2], [5, -1]], time=[0, .3, 1.], controlt=t->sin(t)) -julia> init = _OptimalControlInit(sol) -``` - -""" -mutable struct _OptimalControlInit - - state_init::Function - control_init::Function - variable_init::Union{Nothing, ctVector} - #costate_init::Function - #multipliers_init::Union{Nothing, ctVector} - - """ - $(TYPEDSIGNATURES) - - _OptimalControlInit base constructor with separate explicit arguments - """ - function _OptimalControlInit(; state=nothing, control=nothing, variable=nothing, time=nothing, state_dim=nothing, control_dim=nothing, variable_dim=nothing) - - init = new() - - # some matrix / vector conversions - time = formatTimeGrid(time) - state = formatData(state) - control = formatData(control) - - # set initialization for x, u, v - init.state_init = buildFunctionalInit(state, time, state_dim) - init.control_init = buildFunctionalInit(control, time, control_dim) - init.variable_init = buildVectorInit(variable, variable_dim) - - return init - - end - - """ - $(TYPEDSIGNATURES) - - _OptimalControlInit constructor with arguments grouped as named tuple or dict - """ - function _OptimalControlInit(init_data; state_dim=nothing, control_dim=nothing, variable_dim=nothing) - - # trivial case: default init - x_init = nothing - u_init = nothing - v_init = nothing - t_init = nothing - x_dim = nothing - u_dim = nothing - v_dim = nothing - - # parse arguments - if !isnothing(init_data) - for key in keys(init_data) - if key == :state - x_init = init_data[:state] - elseif key == :control - u_init = init_data[:control] - elseif key == :variable - v_init = init_data[:variable] - elseif key == :time - t_init = init_data[:time] - else - error("Unknown key in initialization data (allowed: state, control, variable, time, state_dim, control_dim, variable_dim): ", key) - end - end - end - - # call base constructor - return _OptimalControlInit(state=x_init, control=u_init, variable=v_init, time=t_init, state_dim=state_dim, control_dim=control_dim, variable_dim=variable_dim) - - end - - """ - $(TYPEDSIGNATURES) - - _OptimalControlInit constructor with solution as argument (warm start) - """ - function _OptimalControlInit(sol::OptimalControlSolution; unused_kwargs...) - return _OptimalControlInit(state=sol.state, control=sol.control, variable=sol.variable, state_dim=sol.state_dimension, control_dim=sol.control_dimension, variable_dim=sol.variable_dimension) - end - -end - -####################################################### -# end of part for CTBase -####################################################### - """ $(TYPEDSIGNATURES) @@ -427,7 +195,7 @@ $(TYPEDSIGNATURES) Build initial guess for discretized problem """ function DOCP_initial_guess(docp, - init::_OptimalControlInit=_OptimalControlInit()) + init::OptimalControlInit=OptimalControlInit()) # default initialization # note: internal variables (lagrange cost, k_i for RK schemes) will keep these default values diff --git a/test/suite/initial_guess.jl b/test/suite/initial_guess.jl index b9239c3..3c39975 100644 --- a/test/suite/initial_guess.jl +++ b/test/suite/initial_guess.jl @@ -38,7 +38,7 @@ v_const = 0.15 x_func = t->[t^2, sqrt(t)] u_func = t->(cos(10*t)+1)*0.5 -# interpolated initial gues +# interpolated initial guess x_vec = [[0, 0], [1, 2], [5, -1]] x_matrix = [0 0; 1 2; 5 -1] u_vec = [0, 0.3, .1] @@ -129,13 +129,12 @@ end # 2 Setting the initial guess at the DOCP level docp = directTranscription(ocp) # mixed init -setInitialGuess(docp, (time=t_vec, state=x_vec, control=u_func, variable=v_const)) -dsol = solve(docp, print_level=0, max_iter=maxiter) @testset verbose = true showtiming = true ":docp_mixed_init" begin + setInitialGuess(docp, (time=t_vec, state=x_vec, control=u_func, variable=v_const)) + dsol = solve(docp, print_level=0, max_iter=maxiter) sol = OCPSolutionFromDOCP(docp, dsol) @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.times[end])) && check_v(sol, v_const)) end - # warm start @testset verbose = true showtiming = true ":docp_warm_start" begin setInitialGuess(docp, sol0)