From 7a7de9168e8f66a49f5087d86028fa81b9300c97 Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Thu, 9 Feb 2023 16:49:03 +0100 Subject: [PATCH 1/7] started constant init --- src/direct/init.jl | 28 ++++++++++++++++++++++++++++ src/direct/problem.jl | 30 +++++++++++++++++++++++++----- src/direct/utils.jl | 4 ++-- 3 files changed, 55 insertions(+), 7 deletions(-) create mode 100644 src/direct/init.jl diff --git a/src/direct/init.jl b/src/direct/init.jl new file mode 100644 index 000000000..722f2e491 --- /dev/null +++ b/src/direct/init.jl @@ -0,0 +1,28 @@ +function set_state_at_time_step!(x, i, dim_x, N, xu) + if i > N + error("trying to set x(t_i) for i > N") + else + xu[1+i*dim_x:(i+1)*dim_x] = x[1:dim_x] + end +end + +function set_control_at_time_step!(u, i, dim_x, N, m, xu) + if i > N + error("trying to set (t_i) for i > N") + else + xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] = u[1:m] + end +end + +function constant_init(dim_x, m, N, dim_xu, x, u) + xu = zeros(dim_xu) + if length(x) != dim_x + error("vector x for initialization should be of size ",dim_x) + if length(u) != m + error("vector u for initialization should be of size ",m) + for i in 1:N+1 + set_state_at_time_step(x, i, dim_x, N, xu) + set_control_at_time_step(u, i, dim_x, N, m, xu) + end + return xu +end \ No newline at end of file diff --git a/src/direct/problem.jl b/src/direct/problem.jl index e2dcb8126..e61302683 100644 --- a/src/direct/problem.jl +++ b/src/direct/problem.jl @@ -1,5 +1,9 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer) + # +++todo: add optional arguments x_init and u_init + # if present init_type is set to "constant_variables" + init_type = "default" + # direct_infos t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, @@ -94,7 +98,7 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer) end # bounds for the constraints - function ipopt_l_u_b() + function constraints_bounds() lb = zeros(nc) ub = zeros(nc) index = 1 # counter for the constraints @@ -134,11 +138,27 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer) return lb, ub end - # todo: init a changer - xu0 = 1.1*ones(dim_xu) + # todo: retrieve from ocp parsed constraints + function variables_bounds() + l_var = -Inf*ones(dim_xu) + u_var = Inf*ones(dim_xu) + return lvar, uvar + end + + # todo: pass optional constant init to ADNLPProblem() + function initial_guess() + if init_type == "default" + # default initialization + xu0 = 1.1*ones(dim_xu) + elseif init_type == "constant_variables" + xu0 = constant_init(dim_x, m, N, dim_xu, x_init, u_init) + return xu0 + end + + # variables bounds + lvar, uvar = variables_bounds() + - l_var = -Inf*ones(dim_xu) - u_var = Inf*ones(dim_xu) if has_free_final_time xu0[end] = 1.0 l_var[end] = 1.e-3 diff --git a/src/direct/utils.jl b/src/direct/utils.jl index 934764e94..9a9b9bbcb 100644 --- a/src/direct/utils.jl +++ b/src/direct/utils.jl @@ -4,7 +4,7 @@ function get_state_at_time_step(xu, i, dim_x, N) x(t_i) """ if i > N - error("trying to access at x(t_i) for i > N") + error("trying to get x(t_i) for i > N") end return xu[1+i*dim_x:(i+1)*dim_x] end @@ -15,7 +15,7 @@ function get_control_at_time_step(xu, i, dim_x, N, m) u(t_i) """ if i > N - error("trying to access at (t_i) for i > N") + error("trying to get (t_i) for i > N") end return xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] end From 7a04748d4b7cc734a6e17a9b165c86d5f8355257 Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Thu, 9 Feb 2023 17:07:15 +0100 Subject: [PATCH 2/7] todo: pass constant init vectors to ADNLProble, --- src/direct/problem.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/direct/problem.jl b/src/direct/problem.jl index e61302683..aef5b5750 100644 --- a/src/direct/problem.jl +++ b/src/direct/problem.jl @@ -138,33 +138,38 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer) return lb, ub end - # todo: retrieve from ocp parsed constraints + # todo: retrieve optional bounds from ocp parsed constraints function variables_bounds() + # unbounded case l_var = -Inf*ones(dim_xu) u_var = Inf*ones(dim_xu) - return lvar, uvar + return l_var, u_var end - # todo: pass optional constant init to ADNLPProblem() + # todo: pass optional constant init to ADNLProblem() function initial_guess() if init_type == "default" # default initialization xu0 = 1.1*ones(dim_xu) elseif init_type == "constant_variables" xu0 = constant_init(dim_x, m, N, dim_xu, x_init, u_init) + end return xu0 end # variables bounds - lvar, uvar = variables_bounds() + l_var, u_var = variables_bounds() + # initial guess + xu0 = initial_guess() + # free final time case if has_free_final_time xu0[end] = 1.0 l_var[end] = 1.e-3 end - lb, ub = ipopt_l_u_b() + lb, ub = constraints_bounds() nlp = ADNLPModel(ipopt_objective, xu0, l_var, u_var, ipopt_constraint, lb, ub) From 66abe41a66728e0ebbfa13d10d2922d84d6b45d7 Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Wed, 15 Feb 2023 14:37:32 +0100 Subject: [PATCH 3/7] test ok sur basic et goddard --- src/direct/init.jl | 28 -------------------- src/direct/problem.jl | 51 ++++++++++++++++++++++++++++++------- src/direct/solve.jl | 4 ++- test/test_basic_manual.jl | 7 ++++- test/test_goddard_manual.jl | 6 ++++- 5 files changed, 56 insertions(+), 40 deletions(-) delete mode 100644 src/direct/init.jl diff --git a/src/direct/init.jl b/src/direct/init.jl deleted file mode 100644 index 722f2e491..000000000 --- a/src/direct/init.jl +++ /dev/null @@ -1,28 +0,0 @@ -function set_state_at_time_step!(x, i, dim_x, N, xu) - if i > N - error("trying to set x(t_i) for i > N") - else - xu[1+i*dim_x:(i+1)*dim_x] = x[1:dim_x] - end -end - -function set_control_at_time_step!(u, i, dim_x, N, m, xu) - if i > N - error("trying to set (t_i) for i > N") - else - xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] = u[1:m] - end -end - -function constant_init(dim_x, m, N, dim_xu, x, u) - xu = zeros(dim_xu) - if length(x) != dim_x - error("vector x for initialization should be of size ",dim_x) - if length(u) != m - error("vector u for initialization should be of size ",m) - for i in 1:N+1 - set_state_at_time_step(x, i, dim_x, N, xu) - set_control_at_time_step(u, i, dim_x, N, m, xu) - end - return xu -end \ No newline at end of file diff --git a/src/direct/problem.jl b/src/direct/problem.jl index aef5b5750..77bfab397 100644 --- a/src/direct/problem.jl +++ b/src/direct/problem.jl @@ -1,8 +1,4 @@ -function ADNLProblem(ocp::OptimalControlModel, N::Integer) - - # +++todo: add optional arguments x_init and u_init - # if present init_type is set to "constant_variables" - init_type = "default" +function ADNLProblem(ocp::OptimalControlModel, N::Integer, init=nothing) # direct_infos t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, @@ -146,13 +142,50 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer) return l_var, u_var end - # todo: pass optional constant init to ADNLProblem() + # generate initial guess + function set_state_at_time_step!(x, i, dim_x, N, xu) + if i > N + error("trying to set x(t_i) for i > N") + else + xu[1+i*dim_x:(i+1)*dim_x] = x[1:dim_x] + end + end + + function set_control_at_time_step!(u, i, dim_x, N, m, xu) + if i > N + error("trying to set (t_i) for i > N") + else + xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] = u[1:m] + end + end + function initial_guess() - if init_type == "default" + #println("Initialization: ", init) + + if init === nothing # default initialization xu0 = 1.1*ones(dim_xu) - elseif init_type == "constant_variables" - xu0 = constant_init(dim_x, m, N, dim_xu, x_init, u_init) + else + if length(init) != (n_x + m) + error("vector for initialization should be of size n+m",n_x+m) + end + # split state / control + x_init = zeros(dim_x) + x_init[1:n_x] = init[1:n_x] + u_init = zeros(m) + u_init[1:m] = init[n_x+1:n_x+m] + + # mayer -> lagrange additional state + if hasLagrangeCost + x_init[dim_x] = 0.1 + end + + # constant initialization + xu0 = zeros(dim_xu) + for i in 0:N + set_state_at_time_step!(x_init, i, dim_x, N, xu0) + set_control_at_time_step!(u_init, i, dim_x, N, m, xu0) + end end return xu0 end diff --git a/src/direct/solve.jl b/src/direct/solve.jl index 0d0be4407..eefa88fc2 100644 --- a/src/direct/solve.jl +++ b/src/direct/solve.jl @@ -3,6 +3,7 @@ function solve(ocp::OptimalControlModel, algo::Direct, method::Description; print_level::Integer=__print_level_ipopt(), mu_strategy::String=__mu_strategy_ipopt(), display::Bool=__display(), + init=nothing, #NB. for now, can be nothing or (n+m) vector kwargs...) """ Solve the optimal control problem @@ -21,7 +22,8 @@ function solve(ocp::OptimalControlModel, algo::Direct, method::Description; print_level = display ? 0 : print_level # from OCP to NLP - nlp = ADNLProblem(ocp, grid_size) + nlp = ADNLProblem(ocp, grid_size, init) + #println("nlp x0:", nlp.meta.x0) # solve by IPOPT ipopt_solution = ipopt(nlp, print_level=print_level, mu_strategy=mu_strategy; kwargs...) diff --git a/test/test_basic_manual.jl b/test/test_basic_manual.jl index 97eaca0f1..aa68f45f8 100644 --- a/test/test_basic_manual.jl +++ b/test/test_basic_manual.jl @@ -19,8 +19,13 @@ B = [ 0.0 constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1]) objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2) # default is to minimise +# initial guess (constant state and control functions) +init = [1., 0.5, 0.3] + # solve -sol = solve(ocp,30) +#sol = solve(ocp, grid_size=10, print_level=5, display=false) +sol = solve(ocp, grid_size=10, print_level=5, display=false, init=init) + # plot plot(sol) diff --git a/test/test_goddard_manual.jl b/test/test_goddard_manual.jl index 985bee7c3..cebd6da35 100644 --- a/test/test_goddard_manual.jl +++ b/test/test_goddard_manual.jl @@ -47,6 +47,10 @@ end constraint!(ocp, :dynamics, f) -sol = solve(ocp, 20) +# initial guess (constant state and control functions) +init = [1.01, 0.25, 0.5, 0.4] + +#sol = solve(ocp, grid_size=20, print_level=5, display=false) +sol = solve(ocp, grid_size=20, print_level=5, display=false, init=init) plot(sol) \ No newline at end of file From 27a4a3420538fda153adf9193873ca0bbbdb2a24 Mon Sep 17 00:00:00 2001 From: Olivier Cots <66357348+ocots@users.noreply.github.com> Date: Wed, 15 Feb 2023 15:36:39 +0100 Subject: [PATCH 4/7] Update runtests.jl --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 3bea62551..1b4a68bf4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,11 +21,11 @@ const __display = OptimalControl.CTBase.__display # @testset verbose = true showtiming = true "Optimal control tests" begin for name in ( - #"utils", - #"direct_shooting_CTOptimization", # unconstrained direct simple shooting + "utils", + "direct_shooting_CTOptimization", # unconstrained direct simple shooting "basic", - #"goddard_direct", - #"goddard_indirect", + "goddard_direct", + "goddard_indirect", ) @testset "$name" begin include("test_$name.jl") From cbde146ca779968c8680747906c03fc6b0feb06f Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Wed, 15 Feb 2023 18:27:29 +0100 Subject: [PATCH 5/7] reverted default init to 1.1; added constant init for goddard test; renamed the 2 manual tests --- src/CTDirect/problem.jl | 2 +- test/{test_basic_manual.jl => manual_test_basic.jl} | 0 test/{test_goddard_manual.jl => manual_test_goddard.jl} | 0 test/test_goddard_direct.jl | 3 ++- 4 files changed, 3 insertions(+), 2 deletions(-) rename test/{test_basic_manual.jl => manual_test_basic.jl} (100%) rename test/{test_goddard_manual.jl => manual_test_goddard.jl} (100%) diff --git a/src/CTDirect/problem.jl b/src/CTDirect/problem.jl index 93c49bd4d..bb4604ef4 100644 --- a/src/CTDirect/problem.jl +++ b/src/CTDirect/problem.jl @@ -164,7 +164,7 @@ function ADNLProblem(ocp::OptimalControlModel, N::Integer, init=nothing) if init === nothing # default initialization - xu0 = 0.1*ones(dim_xu) + xu0 = 1.1*ones(dim_xu) else if length(init) != (n_x + m) error("vector for initialization should be of size n+m",n_x+m) diff --git a/test/test_basic_manual.jl b/test/manual_test_basic.jl similarity index 100% rename from test/test_basic_manual.jl rename to test/manual_test_basic.jl diff --git a/test/test_goddard_manual.jl b/test/manual_test_goddard.jl similarity index 100% rename from test/test_goddard_manual.jl rename to test/manual_test_goddard.jl diff --git a/test/test_goddard_direct.jl b/test/test_goddard_direct.jl index 8e7e7b57a..494439ab6 100644 --- a/test/test_goddard_direct.jl +++ b/test/test_goddard_direct.jl @@ -47,7 +47,8 @@ end constraint!(ocp, :dynamics, f) -sol = solve(ocp, grid_size=10, print_level=0) +init = [1.01, 0.25, 0.5, 0.4] +sol = solve(ocp, grid_size=10, print_level=0, init=init) @test objective(sol) ≈ -1.0 atol=1e-1 @test constraints_violation(sol) < 1e-6 From c836f47b380a117065132838cfd03dee457b323d Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Fri, 17 Feb 2023 10:12:42 +0100 Subject: [PATCH 6/7] clean pkgs --- Project.toml | 3 + src/CTBase/CTBase.jl | 71 ----- src/CTBase/algorithms.jl | 4 - src/CTBase/default.jl | 24 -- src/CTBase/model.jl | 276 -------------------- src/CTBase/print.jl | 87 ------ src/CTBase/solutions.jl | 105 -------- src/CTBase/utils.jl | 69 ----- src/CTDirect/CTDirect.jl | 40 --- src/CTDirect/euler/euler_plot.jl | 43 --- src/CTDirect/euler/euler_problem.jl | 143 ---------- src/CTDirect/euler/euler_solution.jl | 76 ------ src/CTDirect/euler/euler_utils.jl | 87 ------ src/CTDirect/plot.jl | 54 ---- src/CTDirect/problem.jl | 211 --------------- src/CTDirect/solution.jl | 67 ----- src/CTDirect/solve.jl | 41 --- src/CTDirect/utils.jl | 89 ------- src/CTDirectShooting/CTDirectShooting.jl | 60 ----- src/CTDirectShooting/init.jl | 142 ---------- src/CTDirectShooting/plot.jl | 119 --------- src/CTDirectShooting/problem.jl | 61 ----- src/CTDirectShooting/solution.jl | 48 ---- src/CTDirectShooting/solve.jl | 84 ------ src/CTDirectShooting/utils.jl | 49 ---- src/OptimalControl.jl | 14 +- src/{resources => }/flows.jl | 0 src/{resources => }/solve.jl | 4 +- test/Project.toml | 1 + test/manual_test_basic.jl | 31 --- test/manual_test_goddard.jl | 56 ---- test/runtests.jl | 4 +- test/test_basic.jl | 42 --- test/test_direct_shooting_CTOptimization.jl | 164 ------------ test/test_goddard_direct.jl | 51 +--- test/test_goddard_indirect.jl | 3 - test/test_utils.jl | 5 - 37 files changed, 15 insertions(+), 2413 deletions(-) delete mode 100644 src/CTBase/CTBase.jl delete mode 100644 src/CTBase/algorithms.jl delete mode 100644 src/CTBase/default.jl delete mode 100644 src/CTBase/model.jl delete mode 100644 src/CTBase/print.jl delete mode 100644 src/CTBase/solutions.jl delete mode 100644 src/CTBase/utils.jl delete mode 100644 src/CTDirect/CTDirect.jl delete mode 100644 src/CTDirect/euler/euler_plot.jl delete mode 100644 src/CTDirect/euler/euler_problem.jl delete mode 100644 src/CTDirect/euler/euler_solution.jl delete mode 100644 src/CTDirect/euler/euler_utils.jl delete mode 100644 src/CTDirect/plot.jl delete mode 100644 src/CTDirect/problem.jl delete mode 100644 src/CTDirect/solution.jl delete mode 100644 src/CTDirect/solve.jl delete mode 100644 src/CTDirect/utils.jl delete mode 100644 src/CTDirectShooting/CTDirectShooting.jl delete mode 100644 src/CTDirectShooting/init.jl delete mode 100644 src/CTDirectShooting/plot.jl delete mode 100644 src/CTDirectShooting/problem.jl delete mode 100644 src/CTDirectShooting/solution.jl delete mode 100644 src/CTDirectShooting/solve.jl delete mode 100644 src/CTDirectShooting/utils.jl rename src/{resources => }/flows.jl (100%) rename src/{resources => }/solve.jl (87%) delete mode 100644 test/manual_test_basic.jl delete mode 100644 test/manual_test_goddard.jl delete mode 100644 test/test_basic.jl delete mode 100644 test/test_direct_shooting_CTOptimization.jl delete mode 100644 test/test_utils.jl diff --git a/Project.toml b/Project.toml index 04083ca54..9ddfb730c 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,9 @@ version = "0.2.0" [deps] ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" +CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd" +CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5" +CTDirectShooting = "e0e6c04b-5022-4cd2-bea2-4a09fff39444" CTOptimization = "22f08de8-270f-4470-8fba-397dbc90d8e0" ControlToolboxTools = "3a0bcf43-9180-47f3-913c-e71e0d69f39f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" diff --git a/src/CTBase/CTBase.jl b/src/CTBase/CTBase.jl deleted file mode 100644 index 81bac4990..000000000 --- a/src/CTBase/CTBase.jl +++ /dev/null @@ -1,71 +0,0 @@ -module CTBase - -# this should be in ControlToolboxTools, which should be renamed CTBase - -# using -using ForwardDiff: jacobian, gradient, ForwardDiff # automatic differentiation -using Parameters # @with_kw: permit to have default values in struct -using Interpolations: linear_interpolation, Line, Interpolations # for default interpolation -#import Base: show # to print an OptimalControlModel -using Printf # to print a OptimalControlModel -using ControlToolboxTools # tools: callbacks, exceptions, functions and more - -# -------------------------------------------------------------------------------------------------- -# Aliases for types -# const AbstractVector{T} = AbstractArray{T,1}. -const MyNumber = Real -const MyVector = AbstractVector{<:MyNumber} - -const Time = MyNumber -const Times = MyVector -const TimesDisc = Union{MyVector,StepRangeLen} - -const States = Vector{<:MyVector} -const Adjoints = Vector{<:MyVector} -const Controls = Vector{<:MyVector} - -const State = MyVector -const Adjoint = MyVector # todo: ajouter type adjoint pour faire par exemple p*f(x, u) au lieu de p'*f(x,u) -const Dimension = Integer - -# -types() = MyNumber, MyVector, Time, Times, TimesDisc, States, Adjoints, Controls, State, Adjoint, Dimension - -# -include("utils.jl") -#include("algorithms.jl") -include("model.jl") -include("print.jl") -include("solutions.jl") -include("default.jl") - -#function solve(ocp::OptimalControlModel, algo::AbstractControlAlgorithm, method::Description; kwargs...) -# error("solve not implemented") -#end - -# -# export only for users - -# utils -export Ad, Poisson - -# model -export Model, time!, constraint!, objective!, state!, control!, remove_constraint!, constraint -export ismin, dynamics, lagrange, criterion, initial_time, final_time -export control_dimension, state_dimension, constraints, initial_condition, final_constraint - -# solution -export time_steps_length, state_dimension, control_dimension -export time_steps, state, control, adjoint, objective -export iterations, success, message, stopping -export constraints_violation - -# export structs -export AbstractOptimalControlModel, OptimalControlModel -export AbstractOptimalControlSolution, DirectSolution, DirectShootingSolution -#export AbstractControlAlgorithm, DirectAlgorithm, DirectShootingAlgorithm - -# solve -#export solve - -end \ No newline at end of file diff --git a/src/CTBase/algorithms.jl b/src/CTBase/algorithms.jl deleted file mode 100644 index 8d38837e0..000000000 --- a/src/CTBase/algorithms.jl +++ /dev/null @@ -1,4 +0,0 @@ -abstract type AbstractControlAlgorithm end - -struct DirectAlgorithm <: AbstractControlAlgorithm end -struct DirectShootingAlgorithm <: AbstractControlAlgorithm end \ No newline at end of file diff --git a/src/CTBase/default.jl b/src/CTBase/default.jl deleted file mode 100644 index 44a223690..000000000 --- a/src/CTBase/default.jl +++ /dev/null @@ -1,24 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# commun -__matrix_dimension_stock() = 1 - -# -------------------------------------------------------------------------------------------------- -# Direct shooting method - default values -__grid_size_direct_shooting() = 201 -__penalty_constraint() = 1e4 # the penalty term in front of final constraints -__iterations() = 100 # number of maximal iterations -__absoluteTolerance() = 10 * eps() # absolute tolerance for the stopping criterion -__optimalityTolerance() = 1e-8 # optimality relative tolerance for the CN1 -__stagnationTolerance() = 1e-8 # step stagnation relative tolerance -__display() = true # print output during resolution -__callbacks() = () -function __init_interpolation() # default for interpolation of the initialization - return (T, U) -> Interpolations.linear_interpolation(T, U, extrapolation_bc = Interpolations.Line()) -end - -# ------------------------------------------------------------------------------------ -# Direct method - default values -# -__grid_size_direct() = 100 -__print_level_ipopt() = 5 -__mu_strategy_ipopt() = "adaptive" \ No newline at end of file diff --git a/src/CTBase/model.jl b/src/CTBase/model.jl deleted file mode 100644 index 9e22c8396..000000000 --- a/src/CTBase/model.jl +++ /dev/null @@ -1,276 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# Abstract Optimal control model -abstract type AbstractOptimalControlModel end - -@with_kw mutable struct OptimalControlModel{time_dependence} <: AbstractOptimalControlModel - initial_time::Union{Time,Nothing}=nothing - final_time::Union{Time,Nothing}=nothing - lagrange::Union{LagrangeFunction{time_dependence},Nothing}=nothing - mayer::Union{Function,Nothing}=nothing - criterion::Union{Symbol,Nothing}=nothing - dynamics::Union{DynamicsFunction{time_dependence},Nothing}=nothing - dynamics!::Union{Function,Nothing}=nothing - state_dimension::Union{Dimension,Nothing}=nothing - control_dimension::Union{Dimension,Nothing}=nothing - constraints::Dict{Symbol, Tuple{Vararg{Any}}}=Dict{Symbol, Tuple{Vararg{Any}}}() -end - -# -isnonautonomous(time_dependence::Symbol) = :nonautonomous == time_dependence -isautonomous(time_dependence::Symbol) = !isnonautonomous(time_dependence) -isnonautonomous(ocp::OptimalControlModel{time_dependence}) where {time_dependence} = isnonautonomous(time_dependence) -isautonomous(ocp::OptimalControlModel) = !isnonautonomous(ocp) - -# Constructors -abstract type Model{td} end # c'est un peu du bricolage ça -function Model{time_dependence}() where {time_dependence} - return OptimalControlModel{time_dependence}() -end -Model() = Model{:autonomous}() # default value - -# ------------------------------------------------------------------------------------------- -# getters -dynamics(ocp::OptimalControlModel) = ocp.dynamics -lagrange(ocp::OptimalControlModel) = ocp.lagrange -mayer(ocp::OptimalControlModel) = ocp.mayer -criterion(ocp::OptimalControlModel) = ocp.criterion -ismin(ocp::OptimalControlModel) = criterion(ocp) == :min -initial_time(ocp::OptimalControlModel) = ocp.initial_time -final_time(ocp::OptimalControlModel) = ocp.final_time -control_dimension(ocp::OptimalControlModel) = ocp.control_dimension -state_dimension(ocp::OptimalControlModel) = ocp.state_dimension -constraints(ocp::OptimalControlModel) = ocp.constraints -function initial_condition(ocp::OptimalControlModel) - cs = constraints(ocp) - x0 = nothing - for (_, c) ∈ cs - type, _, _, val = c - if type == :initial - x0 = val - end - end - return x0 -end -function final_constraint(ocp::OptimalControlModel) - cs = constraints(ocp) - cf = nothing - for (_, c) ∈ cs - type, _, f, val = c - if type == :final - cf = x -> f(x) - val - end - end - return cf -end - -# ------------------------------------------------------------------------------------------- -# -function state!(ocp::OptimalControlModel, n::Dimension) - ocp.state_dimension = n -end - -function control!(ocp::OptimalControlModel, m::Dimension) - ocp.control_dimension = m -end - -# ------------------------------------------------------------------------------------------- -# -function time!(ocp::OptimalControlModel, type::Symbol, time::Time) - type_ = Symbol(type, :_time) - if type_ ∈ [ :initial_time, :final_time ] - setproperty!(ocp, type_, time) - else - error("this time choice is not valid") - end -end - -function time!(ocp::OptimalControlModel, times::Times) - if length(times) != 2 - error("times must be of dimension 2") - end - ocp.initial_time=times[1] - ocp.final_time=times[2] -end - -# ------------------------------------------------------------------------------------------- -# -function constraint!(ocp::OptimalControlModel, type::Symbol, val::State, label::Symbol=gensym(:anonymous)) - if type ∈ [ :initial, :final ] - ocp.constraints[label] = (type, :eq, x -> x, val) - else - error("this constraint is not valid") - end -end - -function constraint!(ocp::OptimalControlModel, type::Symbol, lb::Real, ub::Real, label::Symbol=gensym(:anonymous)) - if type ∈ [ :initial, :final ] - ocp.constraints[label] = (type, :ineq, x -> x, ub, lb) - else - error("this constraint is not valid") - end -end - -function constraint!(ocp::OptimalControlModel{time_dependence}, type::Symbol, f::Function) where {time_dependence} - if type ∈ [ :dynamics, :dynamics! ] - setproperty!(ocp, type, DynamicsFunction{time_dependence}(f)) - else - error("this constraint is not valid") - end -end - -function constraint!(ocp::OptimalControlModel, type::Symbol, f::Function, val::Real, label::Symbol=gensym(:anonymous)) - if type ∈ [ :control, :mixed, :boundary ] - ocp.constraints[label] = (type, :eq, f, val) - else - error("this constraint is not valid") - end -end - -function constraint!(ocp::OptimalControlModel, type::Symbol, f::Function, lb::Real, ub::Real, label::Symbol=gensym(:anonymous)) - if type ∈ [ :control, :mixed, :boundary ] - ocp.constraints[label] = (type, :ineq, f, lb, ub) - else - error("this constraint is not valid") - end -end - -# -function remove_constraint!(ocp::OptimalControlModel, label::Symbol) - delete!(ocp.constraints, label) -end - -# -function constraint(ocp::OptimalControlModel{time_dependence}, label::Symbol) where {time_dependence} - con = ocp.constraints[label] - if length(con) != 4 - error("this constraint is not valid") - end - type, _, f, val = con - if type ∈ [ :initial, :final ] - return x -> f(x) - val - elseif type == :boundary - return (t0, x0, tf, xf) -> f(t0, x0, tf, xf) - val - elseif type == :control - return isautonomous(time_dependence) ? u -> f(u) - val : (t, u) -> f(t, u) - val - elseif type == :mixed - return isautonomous(time_dependence) ? (x, u) -> f(x, u) - val : (t, x, u) -> f(t, x, u) - val - else - error("this constraint is not valid") - end - return nothing -end - -# -function constraint(ocp::OptimalControlModel{time_dependence}, label::Symbol, bound::Symbol) where {time_dependence} - # constraints are all >= 0 - con = ocp.constraints[label] - if length(con) != 5 - error("this constraint is not valid") - end - type, _, f, lb, ub = con - if !( bound ∈ [ :lower, :upper ] ) - error("this constraint is not valid") - end - if (bound == :lower && lb == -Inf) || (bound == :upper && ub == Inf) - error("this constraint is not valid") - end - if type ∈ [ :initial, :final ] - return bound == :lower ? x -> f(x) - lb : x -> ub - f(x) - elseif type == :boundary - return bound == :lower ? (t0, x0, tf, xf) -> f(t0, x0, tf, xf) - lb : (t0, x0, tf, xf) -> ub - f(t0, x0, tf, xf) - elseif type == :control - if isautonomous(time_dependence) - return bound == :lower ? u -> f(u) - lb : u -> ub - f(u) - else - return bound == :lower ? (t, u) -> f(t, u) - lb : (t, u) -> ub - f(t, u) - end - elseif type == :mixed - if isautonomous(time_dependence) - return bound == :lower ? (x, u) -> f(x, u) - lb : (x, u) -> ub - f(x, u) - else - return bound == :lower ? (t, x, u) -> f(t, x, u) - lb : (t, x, u) -> ub - f(t, x, u) - end - else - error("this constraint is not valid") - end - return nothing -end - -# -function nlp_constraints(ocp::OptimalControlModel{time_dependence}) where {time_dependence} - # - constraints = ocp.constraints - n = ocp.state_dimension - - ξf = Vector{ControlFunction}(); ξl = Vector{MyNumber}(); ξu = Vector{MyNumber}() - ψf = Vector{StateConstraintFunction}(); ψl = Vector{MyNumber}(); ψu = Vector{MyNumber}() - ϕf = Vector{Function}(); ϕl = Vector{MyNumber}(); ϕu = Vector{MyNumber}() - - for (_, c) ∈ constraints - if c[1] == :control - push!(ξf, ControlFunction{time_dependence}(c[3])) - append!(ξl, c[4]) - append!(ξu, c[2] == :eq ? c[4] : c[5]) - elseif c[1] == :mixed - push!(ψf, StateConstraintFunction{time_dependence}(c[3])) - append!(ψl, c[4]) - append!(ψu, c[2] == :eq ? c[4] : c[5]) - elseif c[1] == :initial - push!(ϕf, (t0, x0, tf, xf) -> c[3](x0)) - append!(ϕl, c[4]) - append!(ϕu, c[2] == :eq ? c[4] : c[5]) - elseif c[1] == :final - push!(ϕf, (t0, x0, tf, xf) -> c[3](xf)) - append!(ϕl, c[4]) - append!(ϕu, c[2] == :eq ? c[4] : c[5]) - elseif c[1] == :boundary - push!(ϕf, (t0, x0, tf, xf) -> c[3](t0, x0, tf, xf)) - append!(ϕl, c[4]) - append!(ϕu, c[2] == :eq ? c[4] : c[5]) - end - end - -# ξ!(val, u) = [ val[i] = ξf[i](u) for i ∈ 1:length(ξf) ] -# ψ!(val, x, u) = [ val[i] = ψf[i](x, u) for i ∈ 1:length(ψf) ] -# ϕ!(val, t0, x0, tf, xf) = [ val[i] = ϕf[i](t0, x0, tf, xf) for i ∈ 1:length(ϕf) ] - - function ξ(t, u) - val = Vector{MyNumber}() - for i ∈ 1:length(ξf) append!(val, ξf[i](t, u)) end - return val - end - - function ψ(t, x, u) - val = Vector{MyNumber}() - for i ∈ 1:length(ψf) append!(val, ψf[i](t, x, u)) end - return val - end - - function ϕ(t0, x0, tf, xf) - val = Vector{MyNumber}() - for i ∈ 1:length(ϕf) append!(val, ϕf[i](t0, x0, tf, xf)) end - return val - end - - return (ξl, ξ, ξu), (ψl, ψ, ψu), (ϕl, ϕ, ϕu) - -end - -# ------------------------------------------------------------------------------------------- -# -function objective!(ocp::OptimalControlModel{time_dependence}, type::Symbol, f::Function, criterion::Symbol=:min) where {time_dependence} - setproperty!(ocp, :mayer, nothing) - setproperty!(ocp, :lagrange, nothing) - if criterion ∈ [ :min, :max ] - ocp.criterion = criterion - else - error("this criterion is not valid") - end - if type == :mayer - setproperty!(ocp, :mayer, f) - elseif type == :lagrange - setproperty!(ocp, :lagrange, LagrangeFunction{time_dependence}(f)) - else - error("this objective is not valid") - end -end diff --git a/src/CTBase/print.jl b/src/CTBase/print.jl deleted file mode 100644 index 5a37fe551..000000000 --- a/src/CTBase/print.jl +++ /dev/null @@ -1,87 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# Display: text/html ? -# Base.show, Base.print -# pretty print : https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing -function Base.show(io::IO, ::MIME"text/plain", ocp::OptimalControlModel{time_dependence}) where {time_dependence} - - dimx = state_dimension(ocp) === nothing ? "n" : state_dimension(ocp) - dimu = control_dimension(ocp) === nothing ? "m" : control_dimension(ocp) - - printstyled(io, "\nOptimal control problem of the form:\n", bold=true) - println(io, "") - printstyled(io, " minimize ", color=:blue); print("J(t0, tf, x, u) = ") - if ocp.mayer !== nothing - print(io, " g(t0, x(t0), tf, x(tf))") - end - if ocp.mayer !== nothing && lagrange(ocp) !== nothing - print(io, " +") - end - if lagrange(ocp) !== nothing - isnonautonomous(ocp) ? - println(io, '\u222B', " f⁰(t, x(t), u(t)) dt, over [t0, tf]") : - println(io, '\u222B', " f⁰(x(t), u(t)) dt, over [t0, tf]") - end - println(io, "") - printstyled(io, " subject to\n", color=:blue) - println(io, "") - isnonautonomous(ocp) ? println(io, " x", '\u0307', "(t) = f(t, x(t), u(t)), t in [t0, tf] a.e.,") : println(io, " x", '\u0307', "(t) = f(x(t), u(t)), t in [t0, tf] a.e.,") - println(io, "") - # constraints - ξ, ψ, ϕ = nlp_constraints(ocp) - dim_ξ = length(ξ[1]) # dimension of the boundary constraints - dim_ψ = length(ψ[1]) - dim_ϕ = length(ϕ[1]) - has_ξ = !isempty(ξ[1]) - has_ψ = !isempty(ψ[1]) - has_ϕ = !isempty(ϕ[1]) - if has_ξ - isnonautonomous(ocp) ? - println(io, " ξl ≤ ξ(t, u(t)) ≤ ξu, ") : - println(io, " ξl ≤ ξ(u(t)) ≤ ξu, ") - end - if has_ψ - isnonautonomous(ocp) ? - println(io, " ψl ≤ ψ(t, x(t), u(t)) ≤ ψu, ") : - println(io, " ψl ≤ ψ(x(t), u(t)) ≤ ψu, ") - end - if has_ϕ - println(io, " ϕl ≤ ϕ(t0, x(t0), tf, x(tf)) ≤ ϕu, ") - end - println(io, "") - print(io, " where x(t) ", '\u2208', " R", dimx == 1 ? "" : Base.string("^", dimx)) - print(io, " and u(t) ", '\u2208', " R", dimu == 1 ? "" : Base.string("^", dimu), ".") - println(io, "") - println(io, "") - nb_fixed = 0 - s = "" - if initial_time(ocp) !== nothing - s = s * "t0" - nb_fixed += 1 - end - if final_time(ocp) !== nothing - if s == "" - s = s * "tf" - else - s = s * ", tf" - end - nb_fixed += 1 - end - if initial_condition(ocp) !== nothing - if s == "" - s = s * "x0" - else - s = s * " and x0" - end - nb_fixed += 1 - end - if nb_fixed > 1 - s = s * " are fixed." - elseif nb_fixed == 1 - s = s * " is fixed." - end - if nb_fixed > 0 - println(io, " Besides, ", s) - end - #println(io, "") - -end diff --git a/src/CTBase/solutions.jl b/src/CTBase/solutions.jl deleted file mode 100644 index 9828cf378..000000000 --- a/src/CTBase/solutions.jl +++ /dev/null @@ -1,105 +0,0 @@ -# ------------------------------------------------------------------------------------ -# Abstract solution -# -abstract type AbstractOptimalControlSolution end - -# getters -message_aocs = "method not implemented for solutions of type " -error_aocs(sol::AbstractOptimalControlSolution) = error(message_aocs*String(typeof(sol))) -# -state_dimension(sol::AbstractOptimalControlSolution) = error_aocs(sol) -control_dimension(sol::AbstractOptimalControlSolution) = error_aocs(sol) -time_steps(sol::AbstractOptimalControlSolution) = error_aocs(sol) -time_steps_length(sol::AbstractOptimalControlSolution) = error_aocs(sol) -state(sol::AbstractOptimalControlSolution) = error_aocs(sol) -control(sol::AbstractOptimalControlSolution) = error_aocs(sol) -adjoint(sol::AbstractOptimalControlSolution) = error_aocs(sol) -objective(sol::AbstractOptimalControlSolution) = error_aocs(sol) -iterations(sol::AbstractOptimalControlSolution) = error_aocs(sol) -success(sol::AbstractOptimalControlSolution) = error_aocs(sol) -message(sol::AbstractOptimalControlSolution) = error_aocs(sol) -stopping(sol::AbstractOptimalControlSolution) = error_aocs(sol) -constraints_violation(sol::AbstractOptimalControlSolution) = error_aocs(sol) - -# ------------------------------------------------------------------------------------ -# Direct solution -# -mutable struct DirectSolution - T::Vector{<:MyNumber} - X::Matrix{<:MyNumber} - U::Matrix{<:MyNumber} - P::Matrix{<:MyNumber} - P_ξ::Matrix{<:MyNumber} - P_ψ::Matrix{<:MyNumber} - n::Integer - m::Integer - N::Integer - objective::MyNumber - constraints_violation::MyNumber - iterations::Integer - stats # remove later - #type is https://juliasmoothoptimizers.github.io/SolverCore.jl/stable/reference/#SolverCore.GenericExecutionStats -end - -# getters -# todo: return all variables on a common time grid -# trapeze scheme case: state and control on all time steps [t0,...,tN], as well as path constraints -# only exception is the costate, associated to the dynamics equality constraints: N values instead of N+1 -# we could use the basic extension for the final costate P_N := P_(N-1) (or linear extrapolation) -# NB. things will get more complicated with other discretization schemes (time stages on a different grid ...) -# NEED TO CHOOSE A COMMON OUTPUT GRID FOR ALL SCHEMES (note that building the output can be scheme-dependent) -# PM: I propose to use the time steps [t0, ... , t_N] -# - states are ok, and we can evaluate boundary conditions directly -# - control are technically on a different grid (stages) that can coincide with the steps for some schemes. -# Alternately, the averaged control (in the sense of the butcher coefficients) can be computed on each step. cf bocop -# - adjoint are for each equation linking x(t_i+1) and x(t_i), so always N values regardless of the scheme - -state_dimension(sol::DirectSolution) = sol.n -control_dimension(sol::DirectSolution) = sol.m -time_steps_length(sol::DirectSolution) = sol.N -time_steps(sol::DirectSolution) = sol.T -state(sol::DirectSolution) = sol.X -control(sol::DirectSolution) = sol.U -function adjoint(sol::DirectSolution) - N = sol.N - n = sol.n - P = zeros(N+1, n) - P[1:N,1:n] = sol.P[1:N,1:n] - # trivial constant extrapolation for p(t_f) - P[N+1,1:n] = P[N,1:n] - return P -end -objective(sol::DirectSolution) = sol.objective -constraints_violation(sol::DirectSolution) = sol.constraints_violation -iterations(sol::DirectSolution) = sol.iterations - -# ------------------------------------------------------------------------------------ -# Direct shooting solution -# -struct DirectShootingSolution <: AbstractOptimalControlSolution - T::TimesDisc # the times - X::States # the states at the times T - U::Controls # the controls at T - P::Adjoints # the adjoint at T - objective::MyNumber - state_dimension::Dimension # the dimension of the state - control_dimension::Dimension # the dimension of the control - stopping::Symbol # the stopping criterion - message::String # the message corresponding to the stopping criterion - success::Bool # whether or not the method has finished successfully: CN1, stagnation vs iterations max - iterations::Integer # the number of iterations -end - -# getters -state_dimension(sol::DirectShootingSolution) = sol.state_dimension -control_dimension(sol::DirectShootingSolution) = sol.control_dimension -time_steps(sol::DirectShootingSolution) = sol.T -time_steps_length(sol::DirectShootingSolution) = length(time(sol)) -state(sol::DirectShootingSolution) = sol.X -control(sol::DirectShootingSolution) = sol.U -adjoint(sol::DirectShootingSolution) = sol.P -objective(sol::DirectShootingSolution) = sol.objective -iterations(sol::DirectShootingSolution) = sol.iterations -success(sol::DirectShootingSolution) = sol.success -message(sol::DirectShootingSolution) = sol.message -stopping(sol::DirectShootingSolution) = sol.stopping diff --git a/src/CTBase/utils.jl b/src/CTBase/utils.jl deleted file mode 100644 index 18e3cdd9c..000000000 --- a/src/CTBase/utils.jl +++ /dev/null @@ -1,69 +0,0 @@ -# -# method to compute gradient and Jacobian -ctgradient(f::Function, x) = ForwardDiff.gradient(f, x) -ctjacobian(f::Function, x) = ForwardDiff.jacobian(f, x) - -# transform a Vector{<:Vector{<:MyNumber}} to a Vector{<:MyNumber} -function vec2vec(x::Vector{<:Vector{<:MyNumber}}) - y = x[1] - for i in range(2, length(x)) - y = vcat(y, x[i]) - end - return y -end - -# transform a Vector{<:MyNumber} to a Vector{<:Vector{<:MyNumber}} -function vec2vec(x::Vector{<:MyNumber}, n::Integer) - y = [x[1:n]] - for i in n+1:n:length(x)-n+1 - y = vcat(y, [x[i:i+n-1]]) - end - return y -end - -function expand(x::Vector{<:Vector{<:MyNumber}}) - return vec2vec(x) -end - -function expand(x::Vector{<:MyNumber}) - return x -end - -function expand(x::Matrix{<:MyNumber}) - return expand(matrix2vec(x, 1)) -end - -# transform a Matrix{<:MyNumber} to a Vector{<:Vector{<:MyNumber}} -function matrix2vec(x::Matrix{<:MyNumber}, dim::Integer=__matrix_dimension_stock()) - m, n = size(x) - y = nothing - if dim==1 - y = [x[1,:]] - for i in 2:m - y = vcat(y, [x[i,:]]) - end - else - y = [x[:,1]] - for j in 2:n - y = vcat(y, [x[:,j]]) - end - end - return y -end - -# -function Ad(X, f) - return x -> ctgradient(f, x)'*X(x) -end - -function Poisson(f, g) - function fg(x, p) - n = size(x, 1) - ff = z -> f(z[1:n], z[n+1:2n]) - gg = z -> g(z[1:n], z[n+1:2n]) - df = ctgradient(ff, [ x ; p ]) - dg = ctgradient(gg, [ x ; p ]) - return df[n+1:2n]'*dg[1:n] - df[1:n]'*dg[n+1:2n] - end - return fg -end \ No newline at end of file diff --git a/src/CTDirect/CTDirect.jl b/src/CTDirect/CTDirect.jl deleted file mode 100644 index 30b6f6337..000000000 --- a/src/CTDirect/CTDirect.jl +++ /dev/null @@ -1,40 +0,0 @@ -module CTDirect - -# using -# -#include("../CTBase/CTBase.jl"); -using ..CTBase -import ..CTBase: DirectSolution - -# todo: use RecipesBase instead of plot -using Plots -import Plots: plot, plot! # import instead of using to overload the plot and plot! functions - -# nlp modeling and resolution -using NLPModelsIpopt, ADNLPModels - -# tools: Descriptions, callbacks, exceptions, functions and more -using ControlToolboxTools - -# Types -const MyNumber, MyVector, Time, Times, TimesDisc, States, Adjoints, Controls, State, Adjoint, Dimension = CTBase.types() - -# Other declarations -const nlp_constraints = CTBase.nlp_constraints -const __grid_size_direct = CTBase.__grid_size_direct -const __print_level_ipopt = CTBase.__print_level_ipopt -const __mu_strategy_ipopt = CTBase.__mu_strategy_ipopt -const __display = CTBase.__display - -# includes -include("utils.jl") -include("problem.jl") -include("solve.jl") -include("solution.jl") -include("plot.jl") - -# export functions only for user -export direct_solve -export plot - -end \ No newline at end of file diff --git a/src/CTDirect/euler/euler_plot.jl b/src/CTDirect/euler/euler_plot.jl deleted file mode 100644 index 7d8dd2fbf..000000000 --- a/src/CTDirect/euler/euler_plot.jl +++ /dev/null @@ -1,43 +0,0 @@ -function Plots.plot(sol::DirectSolution) - """ - Plot the solution - - input - sol : direct_sol - """ - T = sol.T - X = sol.X - U = sol.U - P = sol.P - n = sol.n - m = sol.m - N = sol.N - - # state - px = Plots.plot(time, X, layout=(n,1)) - Plots.plot!(px[1], title="state") - Plots.plot!(px[n], xlabel="t") - for i ∈ 1:n - Plots.plot!(px[i], ylabel=string("x_",i)) - end - - # costate - pp = Plots.plot(time[1:N], P,layout = (n,1)) - Plots.plot!(pp[1], title="costate") - Plots.plot!(pp[n], xlabel="t") - for i ∈ 1:n - Plots.plot!(pp[i], ylabel=string("p_",i)) - end - - # control - pu = Plots.plot(time[1:N], U, lc=:red, layout=(m,1)) - for i ∈ 1:m - Plots.plot!(pu[i], ylabel=string("u_",i)) - end - Plots.plot!(pu[1], title="control") - Plots.plot!(pu[m], xlabel="t") - - # main plot - Plots.plot(px, pp, pu, layout=(1,3), legend=false) - -end \ No newline at end of file diff --git a/src/CTDirect/euler/euler_problem.jl b/src/CTDirect/euler/euler_problem.jl deleted file mode 100644 index 42f7e3741..000000000 --- a/src/CTDirect/euler/euler_problem.jl +++ /dev/null @@ -1,143 +0,0 @@ -function ADNLProblem(ocp::OptimalControlModel, N::Integer) - - # direct_infos - t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, - dim_x, nc, dim_xu, f_Mayer, has_free_final_time, criterion = direct_infos(ocp, N) - - # IPOPT objective - function ipopt_objective(xu) - obj = 0 - if hasMayerCost - x0 = get_state_at_time_step(xu, 0, dim_x, N) - xf = get_state_at_time_step(xu, N, dim_x, N) - obj = obj + g(t0, x0[1:n_x], tf, xf[1:n_x]) - end - if hasLagrangeCost - obj = obj + xu[(N+1)*dim_x] - end - return criterion==:min ? obj : -obj - end - - # IPOPT constraints - function ipopt_constraint(xu) - """ - compute the constraints for the NLP : - - discretization of the dynamics via the Euler method - - boundary conditions - inputs - ocp :: ocp model - xu :: - layout of the nlp unknown xu for Euler discretization - additional state variable x_{n+1}(t) for the objective (Lagrange to Mayer formulation) - [x_1(t_0), ... , x_{n+1}(t_0), - ... , - x_{1}(t_N), ... , x_{n+1}(t_N), - u_1(t_0), ... , u_m(t_0), - ... , - u_m(t_{N-1}), ..., u_m(t_{N-1})] - return - c :: - """ - if has_free_final_time - tf = xu[end] - else - tf = tf_ - end - h = (tf-t0)/N - c = zeros(eltype(xu),nc) - # - - # state equation - index = 1 # counter for the constraints - for i in 0:N-1 - # state and control at the current state - xi = get_state_at_time_step(xu, i, dim_x, N) - xip1 = get_state_at_time_step(xu, i+1, dim_x, N) - ui = get_control_at_time_step(xu, i, dim_x, N, m) - # state equation - c[index:index+dim_x-1] = xip1 - (xi + h*f_Mayer(xi, ui)) - index = index + dim_x - if has_ξ - c[index:index+dim_ξ-1] = ξ[2](ui) # ui vector - index = index + dim_ξ - end - if has_ψ - c[index:index+dim_ψ-1] = ψ[2](xi[1:n_x],ui) # ui vector - index = index + dim_ψ - end - end - if has_ψ - xf = get_state_at_time_step(xu, N, dim_x, N) - uf = get_control_at_time_step(xu, N-1, dim_x, N, m) - c[index:index+dim_ψ-1] = ψ[2](xf,uf) # ui is false because Euler - index = index + dim_ψ - end - - # boundary conditions - # ------------------- - x0 = get_state_at_time_step(xu, 0, dim_x, N) - xf = get_state_at_time_step(xu, N, dim_x, N) - c[index:index+dim_ϕ-1] = ϕ[2](t0,x0[1:n_x],tf,xf[1:n_x]) # because Lagrange cost possible - index = index + dim_ϕ - if hasLagrangeCost - c[index] = xu[dim_x] - index = index + 1 - end - - return c - end - - # bounds for the constraints - function ipopt_l_u_b() - lb = zeros(nc) - ub = zeros(nc) - index = 1 # counter for the constraints - for i in 0:N-1 - index = index + dim_x - if has_ξ - lb[index:index+dim_ξ-1] = ξ[1] - ub[index:index+dim_ξ-1] = ξ[3] - index = index + dim_ξ - end - if has_ψ - lb[index:index+dim_ψ-1] = ψ[1] - ub[index:index+dim_ψ-1] = ψ[3] - index = index + dim_ψ - end - end - if has_ψ - lb[index:index+dim_ψ-1] = ψ[1] - ub[index:index+dim_ψ-1] = ψ[3] - index = index + dim_ψ - end - # boundary conditions - lb[index:index+dim_ϕ-1] = ϕ[1] - ub[index:index+dim_ϕ-1] = ϕ[3] - index = index + dim_ϕ - if hasLagrangeCost - lb[index] = 0. - ub[index] = 0. - index = index + 1 - end - - return lb, ub - end - - # todo: init a changer - xu0 = 1.1*ones(dim_xu) - - l_var = -Inf*ones(dim_xu) - u_var = Inf*ones(dim_xu) - if has_free_final_time - xu0[end] = 1.0 - l_var[end] = 1.e-3 - end - - lb, ub = ipopt_l_u_b() - - nlp = ADNLPModel(ipopt_objective, xu0, l_var, u_var, ipopt_constraint, lb, ub) - - return nlp - -end diff --git a/src/CTDirect/euler/euler_solution.jl b/src/CTDirect/euler/euler_solution.jl deleted file mode 100644 index 0197b09df..000000000 --- a/src/CTDirect/euler/euler_solution.jl +++ /dev/null @@ -1,76 +0,0 @@ -mutable struct DirectSolution - T::Vector{<:MyNumber} - X::Matrix{<:MyNumber} - U::Matrix{<:MyNumber} - P::Matrix{<:MyNumber} - P_ξ::Matrix{<:MyNumber} - P_ψ::Matrix{<:MyNumber} - n::Integer - m::Integer - N::Integer - stats #stats::SolverCore.GenericExecutionStats -end - -function DirectSolution(ocp::OptimalControlModel, N::Integer, ipopt_solution) - - # direct_infos - t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, - dim_x, nc, dim_xu, f_Mayer, has_free_final_time, criterion = direct_infos(ocp, N) - - function parse_ipopt_sol(stats) - """ - return - X : matrix(N+1,n+1) - U : matrix(N,m) - P : matrix(N,n+1) - """ - # states and controls - xu = stats.solution - X = zeros(N+1,dim_x) - U = zeros(N,m) - for i in 1:N - X[i,:] = get_state_at_time_step(xu, i-1, dim_x, N) - U[i,:] = get_control_at_time_step(xu, i-1, dim_x, N, m) - end - X[N+1,:] = get_state_at_time_step(xu, N, dim_x, N) - - # adjoints - P = zeros(N, dim_x) - lambda = stats.multipliers - P_ξ = zeros(N,dim_ξ) - P_ψ = zeros(N+1,dim_ψ) - index = 1 # counter for the constraints - for i ∈ 1:N - # state equation - P[i,:] = lambda[index:index+dim_x-1] - index = index + dim_x - if has_ξ - P_ξ[i,:] = lambda[index:index+dim_ξ-1] - index = index + dim_ξ - end - if has_ψ - P_ψ[i,:] = lambda[index:index+dim_ψ-1] - index = index + dim_ψ - end - P_ψ[N+1,:] = lambda[index:index+dim_ψ-1] - end - return X, U, P, P_ξ, P_ψ - end - - # state, control, adjoint - X, U, P, P_ξ, P_ψ = parse_ipopt_sol(ipopt_solution) - - # times - if has_free_final_time - tf = stats.solution[end] - else - tf = tf_ - end - T = collect(t0:(tf-t0)/N:tf) - - # DirectSolution - sol = DirectSolution(T, X, U, P, P_ξ, P_ψ, n_x, m, N, ipopt_solution) - - return sol -end diff --git a/src/CTDirect/euler/euler_utils.jl b/src/CTDirect/euler/euler_utils.jl deleted file mode 100644 index 93b035b37..000000000 --- a/src/CTDirect/euler/euler_utils.jl +++ /dev/null @@ -1,87 +0,0 @@ -function get_state_at_time_step(xu, i, dim_x, N) - """ - return - x(t_i) - """ - if i > N - error("trying to access at x(t_i) for i > N") - end - return xu[1+i*dim_x:(i+1)*dim_x] -end - -function get_control_at_time_step(xu, i, dim_x, N, m) - """ - return - u(t_i) - """ - if i > N-1 - error("trying to access at (t_i) for i > N-1") - end - return xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] -end - -function direct_infos(ocp::OptimalControlModel, N::Integer) - - # Parameters of the Optimal Control Problem - # times - t0 = ocp.initial_time - tf = ocp.final_time - has_free_final_time = isnothing(tf) - # multiplier la dynamique par (tf-t0) - # travailler avec le nouveau temps s dans (0., 1.) - # une fonction t(s) - # dimensions - n_x = ocp.state_dimension - m = ocp.control_dimension - # dynamics - f = ocp.dynamics - # constraints - ξ, ψ, ϕ = OptimalControl.nlp_constraints(ocp) - dim_ξ = length(ξ[1]) # dimension of the boundary constraints - dim_ψ = length(ψ[1]) - dim_ϕ = length(ϕ[1]) - has_ξ = !isempty(ξ[1]) - has_ψ = !isempty(ψ[1]) - has_ϕ = !isempty(ϕ[1]) - - #println("has_ξ = ", has_ξ) - #println("has_ψ = ", has_ψ) - #println("has_ϕ = ", has_ϕ) - - hasLagrangeCost = !isnothing(ocp.lagrange) - hasLagrangeCost ? L = ocp.lagrange : nothing - - hasMayerCost = !isnothing(ocp.mayer) - hasMayerCost ? g = ocp.mayer : nothing - #println("hasLagrangien : ", hasLagrangeCost) - #println("Mayer = ", hasMayerCost) - - # Mayer formulation - # use an additional state for the Lagrange cost - # - # remark : we pass u[1] because in our case ocp.dynamics and ocp.lagrange are defined with a scalar u - # and we consider vectors for x and u in the discretized problem. Note that the same would apply for a scalar x. - # question : how determine if u and x are scalar or vector ? - # second member of the ode for the Mayer formulation - - if hasLagrangeCost - dim_x = n_x + 1 - nc = N*(dim_x+dim_ξ+dim_ψ) + dim_ϕ + 1 + dim_ψ # dimension of the constraints - else - dim_x = n_x - nc = N*(dim_x+dim_ξ+dim_ψ) + dim_ϕ + dim_ψ # dimension of the constraints - end - - dim_xu = (N+1)*dim_x+N*m # dimension the the unknown xu - has_free_final_time ? dim_xu = dim_xu + 1 : nothing - - # todo: cas vectoriel sur u a ajouter - f_Mayer(x, u) = hasLagrangeCost ? [f(x[1:n_x], u[1]); L(x[1:n_x], u[1])] : f(x,u[1]) - - criterion = ocp.criterion - - return t0, tf, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, dim_x, nc, dim_xu, - f_Mayer, has_free_final_time, criterion - -end diff --git a/src/CTDirect/plot.jl b/src/CTDirect/plot.jl deleted file mode 100644 index 2f01698cf..000000000 --- a/src/CTDirect/plot.jl +++ /dev/null @@ -1,54 +0,0 @@ -function plot(sol::DirectSolution) - """ - Plot the solution - - input - sol : direct_sol - """ - - # retrieve info from direct solution - n = CTBase.state_dimension(sol) - m = CTBase.control_dimension(sol) - T = CTBase.time_steps(sol) - X = CTBase.state(sol) - U = CTBase.control(sol) - P = CTBase.adjoint(sol) - obj = CTBase.objective(sol) - cons = CTBase.constraints_violation(sol) - iter = CTBase.iterations(sol) - - println("Objective: ",obj," Constraints: ",cons," Iterations: ",iter) - - # state (plot actual state from ocp ie mask additional state for lagrange cost if present) - # use option for layout 'columns' or 'rows' ? - #px = Plots.plot(T, X[:,1:n], layout=(n,1)) - px = Plots.plot(T, X[:,1:n], layout=(1,n)) - Plots.plot!(px[1], title="state") - Plots.plot!(px[n], xlabel="t") - for i ∈ 1:n - Plots.plot!(px[i], ylabel=string("x_",i)) - end - - # costate - #pp = Plots.plot(T[1:N], P[:,1:n],layout = (n,1)) - pp = Plots.plot(T, P[:,1:n],layout = (1,n)) - Plots.plot!(pp[1], title="costate") - Plots.plot!(pp[n], xlabel="t") - for i ∈ 1:n - Plots.plot!(pp[i], ylabel=string("p_",i)) - end - - # control - #pu = Plots.plot(T, U, lc=:red, layout=(m,1)) - pu = Plots.plot(T, U, lc=:red, layout=(1,m)) - for i ∈ 1:m - Plots.plot!(pu[i], ylabel=string("u_",i)) - end - Plots.plot!(pu[1], title="control") - Plots.plot!(pu[m], xlabel="t") - - # main plot - #Plots.plot(px, pp, pu, layout=(1,3), legend=false) #column layout - Plots.plot(px, pp, pu, layout=(3,1), legend=false) #row layout - -end \ No newline at end of file diff --git a/src/CTDirect/problem.jl b/src/CTDirect/problem.jl deleted file mode 100644 index bb4604ef4..000000000 --- a/src/CTDirect/problem.jl +++ /dev/null @@ -1,211 +0,0 @@ -function ADNLProblem(ocp::OptimalControlModel, N::Integer, init=nothing) - - # direct_infos - t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, - dim_x, nc, dim_xu, g, f_Mayer, has_free_final_time, criterion = direct_infos(ocp, N) - - # IPOPT objective - function ipopt_objective(xu) - tf = get_final_time(xu, tf_, has_free_final_time) - obj = 0 - if hasMayerCost - x0 = get_state_at_time_step(xu, 0, dim_x, N) - xf = get_state_at_time_step(xu, N, dim_x, N) - obj = obj + g(t0, x0[1:n_x], tf, xf[1:n_x]) - end - if hasLagrangeCost - obj = obj + xu[(N+1)*dim_x] - end - return ismin(ocp) ? obj : -obj - end - - # IPOPT constraints - function ipopt_constraint(xu) - """ - compute the constraints for the NLP : - - discretization of the dynamics via the trapeze method - - boundary conditions - inputs - ocp :: ocp model - xu :: - layout of the nlp unknown xu for trapeze discretization - additional state variable x_{n+1}(t) for the objective (Lagrange to Mayer formulation) - [x_1(t_0), ... , x_{n+1}(t_0), - ... , - x_{1}(t_N), ... , x_{n+1}(t_N), - u_1(t_0), ... , u_m(t_0), - ... , - u_m(t_N), ..., u_m(t_N)] - return - c :: - """ - tf = get_final_time(xu, tf_, has_free_final_time) - h = (tf-t0)/N - c = zeros(eltype(xu),nc) - # - - # state equation - index = 1 # counter for the constraints - for i in 0:N-1 - ti = t0 + i*h - tip1 = t0 + i*h + h - # state and control at the current state - xi = get_state_at_time_step(xu, i, dim_x, N) - xip1 = get_state_at_time_step(xu, i+1, dim_x, N) - ui = get_control_at_time_step(xu, i, dim_x, N, m) - uip1 = get_control_at_time_step(xu, i+1, dim_x, N, m) - # state equation - c[index:index+dim_x-1] = xip1 - (xi + 0.5*h*(f_Mayer(ti, xi, ui)+f_Mayer(tip1, xip1, uip1))) - index = index + dim_x - if has_ξ - c[index:index+dim_ξ-1] = ξ[2](ti, ui) # ui vector - index = index + dim_ξ - end - if has_ψ - c[index:index+dim_ψ-1] = ψ[2](ti, xi[1:n_x], ui) # ui vector - index = index + dim_ψ - end - end - if has_ξ - uf = get_control_at_time_step(xu, N, dim_x, N, m) - c[index:index+dim_ξ-1] = ξ[2](tf, uf) - index = index + dim_ξ - end - if has_ψ - xf = get_state_at_time_step(xu, N, dim_x, N) - uf = get_control_at_time_step(xu, N-1, dim_x, N, m) - c[index:index+dim_ψ-1] = ψ[2](tf, xf, uf) # ui is false because Euler - index = index + dim_ψ - end - - # boundary conditions - # ------------------- - x0 = get_state_at_time_step(xu, 0, dim_x, N) - xf = get_state_at_time_step(xu, N, dim_x, N) - c[index:index+dim_ϕ-1] = ϕ[2](t0, x0[1:n_x], tf, xf[1:n_x]) # because Lagrange cost possible - index = index + dim_ϕ - if hasLagrangeCost - c[index] = xu[dim_x] - index = index + 1 - end - - return c - end - - # bounds for the constraints - function constraints_bounds() - lb = zeros(nc) - ub = zeros(nc) - index = 1 # counter for the constraints - for i in 0:N-1 - index = index + dim_x # leave 0 for the state equation - if has_ξ - lb[index:index+dim_ξ-1] = ξ[1] - ub[index:index+dim_ξ-1] = ξ[3] - index = index + dim_ξ - end - if has_ψ - lb[index:index+dim_ψ-1] = ψ[1] - ub[index:index+dim_ψ-1] = ψ[3] - index = index + dim_ψ - end - end - if has_ξ - lb[index:index+dim_ξ-1] = ξ[1] - ub[index:index+dim_ξ-1] = ξ[3] - index = index + dim_ξ - end - if has_ψ - lb[index:index+dim_ψ-1] = ψ[1] - ub[index:index+dim_ψ-1] = ψ[3] - index = index + dim_ψ - end - # boundary conditions - lb[index:index+dim_ϕ-1] = ϕ[1] - ub[index:index+dim_ϕ-1] = ϕ[3] - index = index + dim_ϕ - if hasLagrangeCost - lb[index] = 0. - ub[index] = 0. - index = index + 1 - end - - return lb, ub - end - - # todo: retrieve optional bounds from ocp parsed constraints - function variables_bounds() - # unbounded case - l_var = -Inf*ones(dim_xu) - u_var = Inf*ones(dim_xu) - return l_var, u_var - end - - # generate initial guess - function set_state_at_time_step!(x, i, dim_x, N, xu) - if i > N - error("trying to set x(t_i) for i > N") - else - xu[1+i*dim_x:(i+1)*dim_x] = x[1:dim_x] - end - end - - function set_control_at_time_step!(u, i, dim_x, N, m, xu) - if i > N - error("trying to set (t_i) for i > N") - else - xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] = u[1:m] - end - end - - function initial_guess() - #println("Initialization: ", init) - - if init === nothing - # default initialization - xu0 = 1.1*ones(dim_xu) - else - if length(init) != (n_x + m) - error("vector for initialization should be of size n+m",n_x+m) - end - # split state / control - x_init = zeros(dim_x) - x_init[1:n_x] = init[1:n_x] - u_init = zeros(m) - u_init[1:m] = init[n_x+1:n_x+m] - - # mayer -> lagrange additional state - if hasLagrangeCost - x_init[dim_x] = 0.1 - end - - # constant initialization - xu0 = zeros(dim_xu) - for i in 0:N - set_state_at_time_step!(x_init, i, dim_x, N, xu0) - set_control_at_time_step!(u_init, i, dim_x, N, m, xu0) - end - end - return xu0 - end - - # variables bounds - l_var, u_var = variables_bounds() - - # initial guess - xu0 = initial_guess() - - # free final time case - if has_free_final_time - xu0[end] = 1.0 - l_var[end] = 1.e-3 - end - - lb, ub = constraints_bounds() - - nlp = ADNLPModel(ipopt_objective, xu0, l_var, u_var, ipopt_constraint, lb, ub) - - return nlp - -end diff --git a/src/CTDirect/solution.jl b/src/CTDirect/solution.jl deleted file mode 100644 index 189d66ac5..000000000 --- a/src/CTDirect/solution.jl +++ /dev/null @@ -1,67 +0,0 @@ -function DirectSolution(ocp::OptimalControlModel, N::Integer, ipopt_solution) - - # direct_infos - t0, tf_, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, - dim_x, nc, dim_xu, g, f_Mayer, has_free_final_time, criterion = direct_infos(ocp, N) - - function parse_ipopt_sol(stats) - - # states and controls - xu = stats.solution - X = zeros(N+1,dim_x) - U = zeros(N+1,m) - for i in 1:N+1 - X[i,:] = get_state_at_time_step(xu, i-1, dim_x, N) - U[i,:] = get_control_at_time_step(xu, i-1, dim_x, N, m) - end - - # adjoints - P = zeros(N, dim_x) - lambda = stats.multipliers - P_ξ = zeros(N+1,dim_ξ) - P_ψ = zeros(N+1,dim_ψ) - index = 1 # counter for the constraints - for i ∈ 1:N - # state equation - P[i,:] = lambda[index:index+dim_x-1] # use getter - index = index + dim_x - if has_ξ - P_ξ[i,:] = lambda[index:index+dim_ξ-1] # use getter - index = index + dim_ξ - end - if has_ψ - P_ψ[i,:] = lambda[index:index+dim_ψ-1] # use getter - index = index + dim_ψ - end - end - if has_ξ - P_ξ[N+1,:] = lambda[index:index+dim_ξ-1] # use getter - index = index + dim_ξ - end - if has_ψ - P_ψ[N+1,:] = lambda[index:index+dim_ψ-1] # use getter - index = index + dim_ψ - end - return X, U, P, P_ξ, P_ψ - end - - # state, control, adjoint - X, U, P, P_ξ, P_ψ = parse_ipopt_sol(ipopt_solution) - - # times - tf = get_final_time(ipopt_solution.solution, tf_, has_free_final_time) - T = collect(LinRange(t0, tf, N+1)) - - # misc info - objective = ipopt_solution.objective - constraints_violation = ipopt_solution.primal_feas - iterations = ipopt_solution.iter - #status = ipopt_solution.status this is a 'Symbol' not an int... - - # DirectSolution - #sol = DirectSolution(T, X, U, P, P_ξ, P_ψ, n_x, m, N, ipopt_solution) - sol = CTBase.DirectSolution(T, X, U, P, P_ξ, P_ψ, n_x, m, N, objective, constraints_violation, iterations, ipopt_solution) - - return sol -end diff --git a/src/CTDirect/solve.jl b/src/CTDirect/solve.jl deleted file mode 100644 index d9dd978a8..000000000 --- a/src/CTDirect/solve.jl +++ /dev/null @@ -1,41 +0,0 @@ -function direct_solve(ocp::OptimalControlModel, - #algo::DirectAlgorithm, - method::Description; - grid_size::Integer=__grid_size_direct(), - print_level::Integer=__print_level_ipopt(), - mu_strategy::String=__mu_strategy_ipopt(), - display::Bool=__display(), - init=nothing, #NB. for now, can be nothing or (n+m) vector - kwargs...) - """ - Solve the optimal control problem - - Input : - ocp : functional description of the optimal control problem (cf. ocp.jl) - N : number of time steps for the discretization - Int - - Output - sol : solution of the discretized problem - (time, X, U, n, m, N) - """ - - # no display - print_level = display ? print_level : 0 - - # from OCP to NLP - nlp = ADNLProblem(ocp, grid_size, init) - #println("nlp x0:", nlp.meta.x0) - - # solve by IPOPT: more info at - # https://github.com/JuliaSmoothOptimizers/NLPModelsIpopt.jl/blob/main/src/NLPModelsIpopt.jl#L119 - # options of ipopt: https://coin-or.github.io/Ipopt/OPTIONS.html - # callback: https://github.com/jump-dev/Ipopt.jl#solver-specific-callback - ipopt_solution = ipopt(nlp, print_level=print_level, mu_strategy=mu_strategy; kwargs...) - - # from IPOPT solution to DirectSolution - sol = DirectSolution(ocp, grid_size, ipopt_solution) - - return sol - -end diff --git a/src/CTDirect/utils.jl b/src/CTDirect/utils.jl deleted file mode 100644 index 5fa0a1ee2..000000000 --- a/src/CTDirect/utils.jl +++ /dev/null @@ -1,89 +0,0 @@ -function get_state_at_time_step(xu, i, dim_x, N) - """ - return - x(t_i) - """ - if i > N - error("trying to get x(t_i) for i > N") - end - return xu[1+i*dim_x:(i+1)*dim_x] -end - -function get_control_at_time_step(xu, i, dim_x, N, m) - """ - return - u(t_i) - """ - if i > N - error("trying to get (t_i) for i > N") - end - return xu[1+(N+1)*dim_x+i*m:m+(N+1)*dim_x+i*m] -end - -get_final_time(xu, tf_, has_free_final_time) = has_free_final_time ? xu[end] : tf_ - -function direct_infos(ocp::OptimalControlModel, N::Integer) - - # Parameters of the Optimal Control Problem - # times - t0 = ocp.initial_time - tf = ocp.final_time - has_free_final_time = isnothing(tf) - # multiplier la dynamique par (tf-t0) - # travailler avec le nouveau temps s dans (0., 1.) - # une fonction t(s) - # dimensions - n_x = ocp.state_dimension - m = ocp.control_dimension - # dynamics - f = ocp.dynamics - # constraints - ξ, ψ, ϕ = nlp_constraints(ocp) - dim_ξ = length(ξ[1]) # dimension of the boundary constraints - dim_ψ = length(ψ[1]) - dim_ϕ = length(ϕ[1]) - has_ξ = !isempty(ξ[1]) - has_ψ = !isempty(ψ[1]) - has_ϕ = !isempty(ϕ[1]) - - #println("has_ξ = ", has_ξ) - #println("has_ψ = ", has_ψ) - #println("has_ϕ = ", has_ϕ) - - hasLagrangeCost = !isnothing(ocp.lagrange) - L = ocp.lagrange - - hasMayerCost = !isnothing(ocp.mayer) - g = ocp.mayer - #println("hasLagrange : ", hasLagrangeCost) - #println("Mayer = ", hasMayerCost) - - # Mayer formulation - # use an additional state for the Lagrange cost - # - # remark : we pass u[1] because in our case ocp.dynamics and ocp.lagrange are defined with a scalar u - # and we consider vectors for x and u in the discretized problem. Note that the same would apply for a scalar x. - # question : how determine if u and x are scalar or vector ? - # second member of the ode for the Mayer formulation - - if hasLagrangeCost - dim_x = n_x + 1 - nc = N*(dim_x+dim_ξ+dim_ψ) + (dim_ξ + dim_ψ) + dim_ϕ + 1 # dimension of the constraints - else - dim_x = n_x - nc = N*(dim_x+dim_ξ+dim_ψ) + (dim_ξ + dim_ψ) + dim_ϕ # dimension of the constraints - end - - dim_xu = (N+1)*(dim_x+m) # dimension the the unknown xu - has_free_final_time ? dim_xu = dim_xu + 1 : nothing - - # todo: cas vectoriel sur u a ajouter - f_Mayer(t, x, u) = hasLagrangeCost ? [f(t, x[1:n_x], u); L(t, x[1:n_x], u)] : f(t, x, u) - - criterion = ocp.criterion - - return t0, tf, n_x, m, f, ξ, ψ, ϕ, dim_ξ, dim_ψ, dim_ϕ, - has_ξ, has_ψ, has_ϕ, hasLagrangeCost, hasMayerCost, dim_x, nc, dim_xu, - g, f_Mayer, has_free_final_time, criterion - -end diff --git a/src/CTDirectShooting/CTDirectShooting.jl b/src/CTDirectShooting/CTDirectShooting.jl deleted file mode 100644 index f51f08461..000000000 --- a/src/CTDirectShooting/CTDirectShooting.jl +++ /dev/null @@ -1,60 +0,0 @@ -module CTDirectShooting - -# using -# -# ici CTBase est vu comme un sous-module déjà inclus dans un module parent (OptimalControl) -# du coup, on fait du ..CTBase -using ..CTBase -import ..CTBase: DirectShootingSolution - -# -using LinearAlgebra # for the norm for instance -using Printf # to print iterations results for instance - -# todo: use RecipesBase instead of plot -using Plots -import Plots: plot, plot! # import instead of using to overload the plot and plot! functions - -# tools: Descriptions, callbacks, exceptions, functions and more -using ControlToolboxTools -const ControlToolboxCallbacks = Tuple{Vararg{ControlToolboxCallback}} # todo: handle this better - -# flows -using HamiltonianFlows - -# nlp solvers -using CTOptimization -import CTOptimization: solve, CTOptimizationProblem - -# Types -const MyNumber, MyVector, Time, Times, TimesDisc, States, Adjoints, Controls, State, Adjoint, Dimension = CTBase.types() - -# Other declarations -const nlp_constraints = CTBase.nlp_constraints -const __grid_size_direct_shooting = CTBase.__grid_size_direct_shooting -const __display = CTBase.__display -const __penalty_constraint = CTBase.__penalty_constraint -const __callbacks = CTBase.__callbacks -const __init_interpolation = CTBase.__init_interpolation -const __iterations = CTBase.__iterations -const __absoluteTolerance = CTBase.__absoluteTolerance -const __optimalityTolerance = CTBase.__optimalityTolerance -const __stagnationTolerance = CTBase.__stagnationTolerance -const ctgradient = CTBase.ctgradient -const ctjacobian = CTBase.ctjacobian -const expand = CTBase.expand -const vec2vec = CTBase.vec2vec - -# includes -include("init.jl") -include("utils.jl") -include("problem.jl") -include("solve.jl") -include("solution.jl") -include("plot.jl") - -# export functions only for user -export direct_shooting_solve -export plot, plot! - -end \ No newline at end of file diff --git a/src/CTDirectShooting/init.jl b/src/CTDirectShooting/init.jl deleted file mode 100644 index 99afa2965..000000000 --- a/src/CTDirectShooting/init.jl +++ /dev/null @@ -1,142 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# make a CTOptimizationInit (Unconstrained) -# direct shooting - -# -------------------------------------------------------------------------------------------------- -# check if the given grid (to the interface to the solver) is valid -function check_grid_validity(t0::Time, tf::Time, T::TimesDisc) - # T: t0 ≤ t1 ≤ ... ≤ tf - valid = true - valid = (t0==T[1]) & valid - valid = (tf==T[end]) & valid - valid = (T==sort(T)) & valid - return valid -end - -function check_grid_validity(U::Vector{<:MyVector}, T::TimesDisc) - return length(U) == (length(T) - 1) || length(U) == length(T) -end - -# -------------------------------------------------------------------------------------------------- -# -function my_interpolation(interp::Function, T::TimesDisc, U::Controls, T_::TimesDisc) - u_lin = interp(T, U) - return u_lin.(T_) -end - -# -------------------------------------------------------------------------------------------------- -# default values -function __grid(t0::Time, tf::Time, N::Integer=__grid_size_direct_shooting()) - return range(t0, tf, N) -end -function __init(m::Dimension, N::Integer=__grid_size_direct_shooting()) - return expand([zeros(m) for i in 1:N-1]) -end - -# -------------------------------------------------------------------------------------------------- -# -function reduce_one_if_necessary(U::Controls, T::TimesDisc) - if length(U) == length(T) - return U[1:end-1] - else - return U - end -end - -# convert to get a Vector of numbers for the CTOptimization solver -function convert_init(U::Controls) - return expand(U) -end - -# -function finalize_init(U::Controls, T::TimesDisc) - return convert_init(reduce_one_if_necessary(U, T)), T -end - -# -------------------------------------------------------------------------------------------------- -# matrix -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Tuple{TimesDisc, Matrix{<:MyNumber}}, grid, interp::Function) - T = init[1] - U = init[2] - return CTOptimizationInit(t0, tf, m, (T, matrix2vec(U)), grid, interp) -end - -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, U::Matrix{<:MyNumber}, grid, interp::Function) - return CTOptimizationInit(t0, tf, m, matrix2vec(U), grid, interp) -end - -# init=nothing, grid=nothing => init=default, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Nothing, grid::Nothing, args...) - return __init(m), __grid(t0, tf) -end - -# init=nothing, grid=T => init=zeros(m, N-1), grid=T, with N=length(T) (check validity) -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Nothing, grid::TimesDisc, args...) - if !check_grid_validity(t0, tf, grid) - throw(InconsistentArgument("grid argument is inconsistent with ocp argument")) - end - return __init(m, length(grid)), grid -end - -# init=U, grid=nothing => init=U, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, U::Controls, grid::Nothing, interp::Function) - T = __grid(t0, tf, length(U)+1) - T_ = __grid(t0, tf) - U_ = my_interpolation(interp, T[1:end-1], U, T_) - return finalize_init(U_, T_) -end - -# init=U, grid=T => init=U, grid=T (check validity with ocp and with init) -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Controls, grid::TimesDisc, args...) - if !check_grid_validity(t0, tf, grid) - throw(InconsistentArgument("grid argument is inconsistent with ocp argument")) - end - if !check_grid_validity(init, grid) - throw(InconsistentArgument("grid argument is inconsistent with init argument")) - end - return finalize_init(init, grid) -end - -# init=(T,U), grid=nothing => init=U, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() (check validity with ocp and with U) -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Tuple{TimesDisc,Controls}, grid::Nothing, interp::Function) - T_ = __grid(t0, tf) # default grid - return CTOptimizationInit(t0, tf, m, init, T_, interp) -end - -# init=(T1,U), grid=T2 => init=U, grid=T2 (check validity with ocp (T1, T2) and with U (T1)) -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, init::Tuple{TimesDisc,Controls}, grid::TimesDisc, interp::Function) - T1 = init[1] - U = init[2] - T2 = grid - if !check_grid_validity(t0, tf, T2) - throw(InconsistentArgument("grid argument is inconsistent with ocp argument")) - end - if !check_grid_validity(t0, tf, T1) - throw(InconsistentArgument("init[1] argument is inconsistent with ocp argument")) - end - if !check_grid_validity(U, T1) - throw(InconsistentArgument("init[1] argument is inconsistent with init[2] argument")) - end - U_ = my_interpolation(interp, T1[1:length(U)], U, T2) - return finalize_init(U_, T2) -end - -# -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, S::DirectShootingSolution, grid, interp::Function) - return CTOptimizationInit(t0, tf, m, (time_steps(S), control(S)), grid, interp) -end - -# -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, S::DirectSolution, grid, interp::Function) - return CTOptimizationInit(t0, tf, m, (time_steps(S), control(S)), grid, interp) -end - -# -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, u::Function, T::TimesDisc, interp::Function) - return CTOptimizationInit(t0, tf, m, u.(T), T, interp) -end - -# -function CTOptimizationInit(t0::Time, tf::Time, m::Dimension, u::Function, grid::Nothing, interp::Function) - return CTOptimizationInit(t0, tf, m, u, __grid(t0, tf), interp) -end \ No newline at end of file diff --git a/src/CTDirectShooting/plot.jl b/src/CTDirectShooting/plot.jl deleted file mode 100644 index e9d364c39..000000000 --- a/src/CTDirectShooting/plot.jl +++ /dev/null @@ -1,119 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# Plot solution -# print("x", '\u2080'+9) : x₉ -# - -# General plot -""" - Plots.plot(sol::UncFreeXfSolution, args...; - state_style=(), - control_style=(), - adjoint_style=(), kwargs...) - -TBW -""" -function plot(sol::DirectShootingSolution, args...; - state_style=(), control_style=(), adjoint_style=(), kwargs...) - - # todo : gérer le cas dans les labels où m, n > 9 - - n = sol.state_dimension - m = sol.control_dimension - - px = Plots.plot(; xlabel="time", title="state", state_style...) - if n == 1 - Plots.plot!(px, sol, :time, (:state, i); label="x", state_style...) - else - for i in range(1, n) - Plots.plot!(px, sol, :time, (:state, i); label="x" * ('\u2080' + i), state_style...) - end - end - - pu = Plots.plot(; xlabel="time", title="control", control_style...) - if m == 1 - Plots.plot!(pu, sol, :time, (:control, 1); label="u", control_style...) - else - for i in range(1, m) - Plots.plot!(pu, sol, :time, (:control, i); label="u" * ('\u2080' + i), control_style...) - end - end - - pp = Plots.plot(; xlabel="time", title="adjoint", adjoint_style...) - if n == 1 - Plots.plot!(pp, sol, :time, (:adjoint, i); label="p", adjoint_style...) - else - for i in range(1, n) - Plots.plot!(pp, sol, :time, (:adjoint, i); label="p" * ('\u2080' + i), adjoint_style...) - end - end - - ps = Plots.plot(px, pu, pp, args..., layout=(1, 3); kwargs...) - - return ps - -end - -# specific plot -function plot(sol::DirectShootingSolution, - xx::Union{Symbol,Tuple{Symbol,Integer}}, yy::Union{Symbol,Tuple{Symbol,Integer}}, args...; kwargs...) - - x = get(sol, xx) - y = get(sol, yy) - - return Plots.plot(x, y, args...; kwargs...) - -end - -function plot!(p::Plots.Plot{<:Plots.AbstractBackend}, sol::DirectShootingSolution, - xx::Union{Symbol,Tuple{Symbol,Integer}}, yy::Union{Symbol,Tuple{Symbol,Integer}}, args...; kwargs...) - - x = get(sol, xx) - y = get(sol, yy) - - Plots.plot!(p, x, y, args...; kwargs...) - -end -#plot!(p, x, y, args...; kwargs...) = Plots.plot!(p, x, y, args...; kwargs...) - -""" - get(sol::UncFreeXfSolution, xx::Union{Symbol, Tuple{Symbol, Integer}}) - -TBW -""" -function get(sol::DirectShootingSolution, xx::Union{Symbol,Tuple{Symbol,Integer}}) - - T = sol.T - X = sol.X - U = sol.U - P = sol.P - - m = length(T) - - if typeof(xx) == Symbol - vv = xx - if vv == :time - x = T - elseif vv == :state - x = [X[i][1] for i in 1:m] - elseif vv == :adjoint || vv == :costate - x = [P[i][1] for i in 1:m] - else - x = vcat([U[i][1] for i in 1:m-1], U[m-1][1]) - end - else - vv = xx[1] - ii = xx[2] - if vv == :time - x = T - elseif vv == :state - x = [X[i][ii] for i in 1:m] - elseif vv == :adjoint || vv == :costate - x = [P[i][ii] for i in 1:m] - else - x = vcat([U[i][ii] for i in 1:m-1], U[m-1][ii]) - end - end - - return x - -end diff --git a/src/CTDirectShooting/problem.jl b/src/CTDirectShooting/problem.jl deleted file mode 100644 index f2d364b66..000000000 --- a/src/CTDirectShooting/problem.jl +++ /dev/null @@ -1,61 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# make an CTOptimizationProblem (Unconstrained) from OptimalControlModel -# direct simple shooting -function CTOptimizationProblem(ocp::OptimalControlModel, grid::TimesDisc, penalty_constraint::Real) - - # - VFN = VectorField{:nonautonomous} - - # parsing ocp - dy, co, cf, x0, n, m = parse_ocp_direct_shooting(ocp) - - # Jacobian of the constraints - Jcf(x) = ctjacobian(cf, x) - - # penalty term for the final constraints - αₚ = penalty_constraint - - # state flow - f = Flow(VFN(dy)) - - # augmented state flow - fa = Flow(VFN((t, x, u) -> [dy(t, x[1:end-1], u)[:]; co(t, x[1:end-1], u)])) - - # state-costate flow - p⁰ = -1.0 - H(t, x, p, u) = p⁰ * co(t, x, u) + p' * dy(t, x, u) - fh = Flow(Hamiltonian{:nonautonomous}(H)) - - # to compute the gradient of the function by the adjoint method, - # we need the partial derivative of the Hamiltonian wrt to the control - Hu(t, x, p, u) = ctgradient(u -> H(t, x, p, u), u) - - # discretization grid - T = grid - - # gradient of the function J - function ∇J(U::Controls) - xₙ, _ = model_primal_forward(x0, T, U, f) - pₙ = p⁰ * αₚ * transpose(Jcf(xₙ)) * cf(xₙ) - _, _, X, P = model_adjoint_backward(xₙ, pₙ, T, U, fh) - g = [-Hu(T[i], X[i], P[i], U[i]) .* (T[i+1] - T[i]) for i in 1:length(T)-1] - return g - end - # vec2vec permet de passer d'un vecteur de vecteur à simplement un vecteur - ∇J(x::Vector{<:Real}) = vec2vec(∇J(vec2vec(x, m))) # for desent solver - - # function J, that we minimize - function J(U::Controls) - # via augmented system - xₙ, X = model_primal_forward([x0[:]; 0.0], T, U, fa) - cost = xₙ[end] + 0.5 * αₚ * norm(cf(xₙ[1:end-1]))^2 - return cost - end - J(x::Vector{<:Real}) = J(vec2vec(x, m)) # for descent solver - - # CTOptimization (Unconstrained) problem - prob = CTOptimization.CTOptimizationProblem(J, gradient=∇J, dimension=length(grid)*m) - - return prob - -end \ No newline at end of file diff --git a/src/CTDirectShooting/solution.jl b/src/CTDirectShooting/solution.jl deleted file mode 100644 index b49c2135c..000000000 --- a/src/CTDirectShooting/solution.jl +++ /dev/null @@ -1,48 +0,0 @@ -function DirectShootingSolution(sol::CTOptimization.UnconstrainedSolution, - ocp::OptimalControlModel, grid::TimesDisc, penalty_constraint::Real) - - # - VFN = VectorField{:nonautonomous} - - # parsing ocp - dy, co, cf, x0, n, m = parse_ocp_direct_shooting(ocp) - - # control solution - U⁺ = vec2vec(sol.x, m) - - # Jacobian of the constraints - Jcf(x) = ctjacobian(cf, x) - - # penalty term for final constraints - αₚ = penalty_constraint - - # state flow - f = Flow(VFN(dy)) - - # augmented state flow - fa = Flow(VFN((t, x, u) -> [dy(t, x[1:end-1], u)[:]; co(t, x[1:end-1], u)])) - - # flow for state-adjoint - p⁰ = -1.0 - H(t, x, p, u) = p⁰ * co(t, x, u) + p' * dy(t, x, u) - fh = Flow(Hamiltonian{:nonautonomous}(H)) - - # get state and adjoint - T = grid - xₙ, _ = model_primal_forward(x0, T, U⁺, f) - pₙ = p⁰ * αₚ * transpose(Jcf(xₙ)) * cf(xₙ) - _, _, X⁺, P⁺ = model_adjoint_backward(xₙ, pₙ, T, U⁺, fh) - - # function J, that we minimize - function J(U::Controls) - # via augmented system - xₙ, X = model_primal_forward([x0[:]; 0.0], T, U, fa) - cost = xₙ[end] + 0.5 * αₚ * norm(cf(xₙ[1:end-1]))^2 - return cost - end - objective = J(U⁺) - - return CTBase.DirectShootingSolution(T, X⁺, U⁺, P⁺, objective, n, m, - sol.stopping, sol.message, sol.success, sol.iterations) - -end \ No newline at end of file diff --git a/src/CTDirectShooting/solve.jl b/src/CTDirectShooting/solve.jl deleted file mode 100644 index 4c6d2fbc1..000000000 --- a/src/CTDirectShooting/solve.jl +++ /dev/null @@ -1,84 +0,0 @@ -#-------------------------------------------------------------------------------------------------- -# Solver of an prob by unconstrained direct simple shooting -function direct_shooting_solve( - ocp::OptimalControlModel, - #algo::DirectShootingAlgorithm, - method::Description; - init::Union{Nothing,Controls,Tuple{TimesDisc,Controls},Function,DirectShootingSolution,DirectSolution}=nothing, - grid::Union{Nothing,TimesDisc}=nothing, - penalty_constraint::Real=__penalty_constraint(), - display::Bool=__display(), - callbacks::ControlToolboxCallbacks=__callbacks(), - init_interpolation::Function=__init_interpolation(), - kwargs... -) - - # check if we can solve this kind of problems, that is there must be no constraints - # - # parse prob - # make CTOptimizationInit: Vector{<:Real} - # make CTOptimizationProblem: CTOptimizationProblem(f; gradient, dimension) - # resolution with CTOptimization - # make DirectShootingSolution - # - # struct UnconstrainedSolution <: CTOptimizationSolution - # x::Primal - # stopping::Symbol - # message::String - # success::Bool - # iterations::Integer - #end - - # check validity - ξ, ψ, ϕ = nlp_constraints(ocp) - dim_ξ = length(ξ[1]) # dimension of the boundary constraints - dim_ψ = length(ψ[1]) - if dim_ξ != 0 && dim_ψ != 0 - error("direct shooting is implemented for problems without constraints") - end - - # Init - need some parsing - t0 = initial_time(ocp) - tf = final_time(ocp) - m = control_dimension(ocp) - opti_init, grid = CTOptimizationInit(t0, tf, m, init, grid, init_interpolation) - - # Problem - nlp = CTOptimizationProblem(ocp, grid, penalty_constraint) - - # Resolution - # callbacks - cbs_print = get_priority_print_callbacks((PrintCallback(printOCPDescent, priority=0), callbacks...)) - cbs_stop = get_priority_stop_callbacks(callbacks) - # - nlp_sol = CTOptimization.solve( - nlp, - method, - init=opti_init, - iterations=__iterations(), - absoluteTolerance=__absoluteTolerance(), - optimalityTolerance=__optimalityTolerance(), - stagnationTolerance=__stagnationTolerance(), - display=display, - callbacks=(cbs_print..., cbs_stop...); - kwargs... - ) - - # transcription of the solution, from descent to prob - sol = DirectShootingSolution(nlp_sol, ocp, grid, penalty_constraint) - - return sol - -end - -#-------------------------------------------------------------------------------------------------- -# print callback for prob resolution by descent method -function printOCPDescent(i, sᵢ, dᵢ, Uᵢ, gᵢ, fᵢ) - if i == 0 - println("\n Calls ‖∇F(U)‖ ‖U‖ Stagnation \n") - end - @printf("%10d", i) # Iterations - @printf("%16.8e", norm(gᵢ)) # ‖∇F(U)‖ - @printf("%16.8e", norm(Uᵢ)) # ‖U‖ - @printf("%16.8e", norm(Uᵢ) > 1e-14 ? norm(sᵢ * dᵢ) / norm(Uᵢ) : norm(sᵢ * dᵢ)) # Stagnation -end \ No newline at end of file diff --git a/src/CTDirectShooting/utils.jl b/src/CTDirectShooting/utils.jl deleted file mode 100644 index fc6e5ec01..000000000 --- a/src/CTDirectShooting/utils.jl +++ /dev/null @@ -1,49 +0,0 @@ -# -------------------------------------------------------------------------------------------------- -# Utils for the transcription from ocp to descent problem - -function parse_ocp_direct_shooting(ocp::OptimalControlModel) - - # parsing ocp - dy = dynamics(ocp) - co = lagrange(ocp) - cf = final_constraint(ocp) - x0 = initial_condition(ocp) - n = state_dimension(ocp) - m = control_dimension(ocp) - - return dy, co, cf, x0, n, m - -end - -# forward integration of the state -""" - model(x0, T, U, f) - -TBW -""" -function model_primal_forward(x0, T, U, f) - xₙ = x0 - X = [xₙ] - for n in range(1, length(T) - 1) - xₙ = f(T[n], xₙ, T[n+1], U[n]) - X = vcat(X, [xₙ]) # vcat gives a vector of vector - end - return xₙ, X -end - -# backward integration of state and costate -""" - adjoint(xₙ, pₙ, T, U, f) - -TBW -""" -function model_adjoint_backward(xₙ, pₙ, T, U, f) - X = [xₙ] - P = [pₙ] - for n in range(length(T), 2, step=-1) - xₙ, pₙ = f(T[n], xₙ, pₙ, T[n-1], U[n-1]) - X = vcat([xₙ], X) - P = vcat([pₙ], P) - end - return xₙ, pₙ, X, P -end \ No newline at end of file diff --git a/src/OptimalControl.jl b/src/OptimalControl.jl index 6c92e8613..cfd04b8b2 100644 --- a/src/OptimalControl.jl +++ b/src/OptimalControl.jl @@ -4,13 +4,9 @@ module OptimalControl using Reexport # include modules -include("./CTBase/CTBase.jl"); -include("CTDirect/CTDirect.jl"); -include("CTDirectShooting/CTDirectShooting.jl"); -# -@reexport using .CTBase -using .CTDirect -using .CTDirectShooting +@reexport using CTBase +using CTDirect +using CTDirectShooting # tools: callbacks, exceptions, functions and more @reexport using ControlToolboxTools @@ -26,8 +22,8 @@ const MyNumber, MyVector, Time, Times, TimesDisc, States, Adjoints, Controls, St const __display = CTBase.__display # resources -include("resources/flows.jl") -include("resources/solve.jl") +include("flows.jl") +include("solve.jl") # export functions only for user export solve diff --git a/src/resources/flows.jl b/src/flows.jl similarity index 100% rename from src/resources/flows.jl rename to src/flows.jl diff --git a/src/resources/solve.jl b/src/solve.jl similarity index 87% rename from src/resources/solve.jl rename to src/solve.jl index d1bd57766..6765e3fe7 100644 --- a/src/resources/solve.jl +++ b/src/solve.jl @@ -26,9 +26,9 @@ function solve(prob::OptimalControlModel, description...; # if no error before, then the method is correct: no need of else if :direct ∈ method if :shooting ∈ method - return direct_shooting_solve(prob, clean(method); display=display, kwargs...) + return CTDirectShooting.solve(prob, clean(method); display=display, kwargs...) else - return direct_solve(prob, clean(method); display=display, kwargs...) + return CTDirect.solve(prob, clean(method); display=display, kwargs...) end end end diff --git a/test/Project.toml b/test/Project.toml index be3692421..91a5d9a08 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +CTProblemLibrary = "0649932a-8c77-4f67-b1e4-c19ddd080280" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/test/manual_test_basic.jl b/test/manual_test_basic.jl deleted file mode 100644 index c878fba8e..000000000 --- a/test/manual_test_basic.jl +++ /dev/null @@ -1,31 +0,0 @@ -t0 = 0. -tf = 1. - -ocp = Model() - -# -state!(ocp, 2) # dimension of the state -control!(ocp, 1) # dimension of the control -time!(ocp, [t0, tf]) -constraint!(ocp, :initial, [ -1.0, 0.0 ]) -constraint!(ocp, :final, [ 0.0, 0.0 ]) - -# -A = [ 0.0 1.0 - 0.0 0.0 ] -B = [ 0.0 - 1.0 ] - -constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1]) -objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2) # default is to minimise - -# initial guess (constant state and control functions) -init = [1., 0.5, 0.3] - -# solve -#sol = solve(ocp, grid_size=10, print_level=5) -sol = solve(ocp, grid_size=10, print_level=5, init=init) - - -# plot -plot(sol) diff --git a/test/manual_test_goddard.jl b/test/manual_test_goddard.jl deleted file mode 100644 index 912e28f29..000000000 --- a/test/manual_test_goddard.jl +++ /dev/null @@ -1,56 +0,0 @@ -# Parameters -Cd = 310. -Tmax = 3.5 -β = 500. -b = 2. -N = 100 -t0 = 0. -r0 = 1. -v0 = 0. -vmax = 0.1 -m0 = 1. -mf = 0.6 -x0 = [ r0, v0, m0 ] - -epsilon = 1.e-2 - -ocp = Model() - -time!(ocp, :initial, t0) # if not provided, final time is free -state!(ocp, 3) # state dim -control!(ocp, 1) # control dim - -constraint!(ocp, :initial, x0) # vectorial equality -constraint!(ocp, :control, u -> u[1], 0., 1.) # constraints can be labeled or not -constraint!(ocp, :mixed, (x, u) -> x[1], r0, Inf, :state_con1) # epsilon car sinon adjoint state at 0 false because the state is r0 touch the constraint at 0 -constraint!(ocp, :mixed, (x, u) -> x[2], 0., vmax, :state_con2) -constraint!(ocp, :mixed, (x, u) -> x[3], mf, m0, :state_con3) - -objective!(ocp, :mayer, (t0, x0, tf, xf) -> xf[1], :max) - -function F0(x) - r, v, m = x - D = Cd * v^2 * exp(-β*(r-1.)) - F = [ v, -D/m-1.0/r^2, 0. ] - return F -end - -function F1(x) - r, v, m = x - F = [ 0., Tmax/m, -b*Tmax ] - return F -end - -function f(x, u) - return F0(x) + u[1]*F1(x) -end - -constraint!(ocp, :dynamics, f) - -# initial guess (constant state and control functions) -init = [1.01, 0.25, 0.5, 0.4] - -#sol = solve(ocp, grid_size=20, print_level=5) -sol = solve(ocp, grid_size=20, print_level=5, init=init) - -plot(sol) diff --git a/test/runtests.jl b/test/runtests.jl index 1b4a68bf4..7a8a0ad15 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using Plots using LinearAlgebra using MINPACK +using CTProblemLibrary # functions and types that are not exported @@ -21,9 +22,6 @@ const __display = OptimalControl.CTBase.__display # @testset verbose = true showtiming = true "Optimal control tests" begin for name in ( - "utils", - "direct_shooting_CTOptimization", # unconstrained direct simple shooting - "basic", "goddard_direct", "goddard_indirect", ) diff --git a/test/test_basic.jl b/test/test_basic.jl deleted file mode 100644 index 60d3a3ae0..000000000 --- a/test/test_basic.jl +++ /dev/null @@ -1,42 +0,0 @@ -t0 = 0. -tf = 1. - -ocp = Model() - -# -state!(ocp, 2) # dimension of the state -control!(ocp, 1) # dimension of the control -time!(ocp, [t0, tf]) -constraint!(ocp, :initial, [ -1.0, 0.0 ]) -constraint!(ocp, :final, [ 0.0, 0.0 ]) - -# -A = [ 0.0 1.0 - 0.0 0.0 ] -B = [ 0.0 - 1.0 ] - -constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1]) -objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2) # default is to minimise - -# -@test ocp.state_dimension == 2 -@test ocp.control_dimension == 1 -@test typeof(ocp) == OptimalControlModel{:autonomous} -@test ocp.initial_time == t0 -@test ocp.final_time == tf -@test ocp.dynamics(0.0, [0.; 1.], 1.0) ≈ [1.; 1.] atol=1e-8 -@test ocp.lagrange(0.0, [0.; 0.], 1.0) ≈ 0.5 atol=1e-8 - -# solve -sol = solve(ocp, print_level=0) - -# solution -u_sol(t) = 6.0-12.0*t -U = sol.U -T = sol.T -dT = T[2:end]-T[1:end-1] - -@test sum(dT .* abs.(U[1:end-1] - u_sol.(T[1:end-1]))) ≈ 0 atol=1e-1 -@test objective(sol) ≈ 6.0 atol=1e-1 -@test constraints_violation(sol) < 1e-6 diff --git a/test/test_direct_shooting_CTOptimization.jl b/test/test_direct_shooting_CTOptimization.jl deleted file mode 100644 index b8db8cf00..000000000 --- a/test/test_direct_shooting_CTOptimization.jl +++ /dev/null @@ -1,164 +0,0 @@ -# -m = 1 -n = 2 -# description of the optimal control problem -t0 = 0 -tf = 1 -x0 = [ -1, 0 ] -xf = [ 0, 0 ] -ocp = Model() -state!(ocp, n) # dimension of the state -control!(ocp, m) # dimension of the control -time!(ocp, [t0, tf]) -constraint!(ocp, :initial, x0) -constraint!(ocp, :final, xf) -A = [ 0 1 - 0 0 ] -B = [ 0 - 1 ] -constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1]) -objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2) - -# solution -u_sol(t) = 6.0-12.0*t - -# -------------------------------------------------------------------------------------------------- -# - -# init=nothing, grid=nothing => init=zeros(m, N-1), grid=range(t0, tf, N), with N=__grid_size_direct_shooting() -init_, grid_ = CTOptimizationInit(t0, tf, m, nothing, nothing) -@test init_ == __init(m) -@test grid_ == __grid(t0, tf) - -# init=nothing, grid=T => init=zeros(m, N-1), grid=T, with N=length(T) (check validity) -N = floor(Int64, __grid_size_direct_shooting()/2) # pour avoir un N différent de __grid_size_direct_shooting() -grid = range(t0, tf, N) -init_, grid_ = CTOptimizationInit(t0, tf, m, nothing, grid) -N = length(grid) -@test init_ == __init(m, N) -@test grid_ == grid -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, nothing, range(t0-1, tf, N)) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, nothing, range(t0, tf+1, N)) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, nothing, range(tf, t0, N)) - -# init=U, grid=nothing => init=U, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() -N = floor(Int64, __grid_size_direct_shooting()/2) # pour avoir un N différent de __grid_size_direct_shooting() -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] -init_, grid_ = CTOptimizationInit(t0, tf, m, U, nothing, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:__grid_size_direct_shooting()-1]) atol=1e-4 -@test grid_ == __grid(t0, tf) - -# init=U, grid=T => init=U, grid=T (check validity with ocp and with init) -N = floor(Int64, __grid_size_direct_shooting()/2) # pour avoir un N différent de __grid_size_direct_shooting() -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] -init_, grid_ = CTOptimizationInit(t0, tf, m, U, T) -@test init_ == convert_init(U) -@test grid_ == T -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, U, range(t0-1, tf, N)) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, U, range(t0, tf+1, N)) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, U, range(tf, t0, N)) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, U, range(t0, tf, 2*N)) - -# init=(T,U), grid=nothing => init=U, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() (check validity with ocp and with U) -N = floor(Int64, __grid_size_direct_shooting()/2) # pour avoir un N différent de __grid_size_direct_shooting() -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] -init_, grid_ = CTOptimizationInit(t0, tf, m, (T,U), nothing, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:__grid_size_direct_shooting()-1]) atol=1e-4 -@test grid_ == __grid(t0, tf) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0-1, tf, N), U), nothing, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0, tf+1, N), U), nothing, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(tf, t0, N), U), nothing, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0, tf, 2*N), U), nothing, __init_interpolation()) - -# init=(T1,U), grid=T2 => init=U, grid=T2 (check validity with ocp (T1, T2) and with U (T1)) -N1 = floor(Int64, __grid_size_direct_shooting()/2) # pour avoir un N différent de __grid_size_direct_shooting() -T1 = range(t0, tf, N1) -U = [[u_sol(T[i])-1.0] for i = 1:N1-1] -N2 = floor(Int64, __grid_size_direct_shooting()/4) # pour avoir un N différent de __grid_size_direct_shooting() -T2 = range(t0, tf, N2) -init_, grid_ = CTOptimizationInit(t0, tf, m, (T1,U), T2, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:N2-1]) atol=1e-4 -@test grid_ == T2 -# T1 -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0-1, tf, N), U), T2, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0, tf+1, N), U), T2, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(tf, t0, N), U), T2, __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (range(t0, tf, 2*N), U), T2, __init_interpolation()) -# T2 -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (T1, U), range(t0-1, tf, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (T1, U), range(t0, tf+1, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, (T1, U), range(tf, t0, N), __init_interpolation()) - -# init=S, grid=nothing => init=S.U, grid=range(t0, tf, N), with N=__grid_size_direct_shooting() -N = floor(Int64, __grid_size_direct_shooting()/2) -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] -sol = solve(ocp, :descent, init=U, grid=T, iterations=0, display=false) -init_, grid_ = CTOptimizationInit(t0, tf, m, sol, nothing, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:__grid_size_direct_shooting()-1]) atol=1e-4 -@test grid_ == __grid(t0, tf) - -# init=S, grid=T => init=S.U, grid=T (check validity with ocp) -N = floor(Int64, __grid_size_direct_shooting()/2) -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] -sol = solve(ocp, :descent, init=U, grid=T, iterations=0, display=false) -N2 = floor(Int64, __grid_size_direct_shooting()/4) # pour avoir un N différent de __grid_size_direct_shooting() -T2 = range(t0, tf, N2) -init_, grid_ = CTOptimizationInit(t0, tf, m, sol, T2, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:N2-1]) atol=1e-4 -@test grid_ == T2 - -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, sol, range(t0-1, tf, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, sol, range(t0, tf+1, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, sol, range(tf, t0, N), __init_interpolation()) - -# init=u, grid=nothing => init=u(T), grid=T=range(t0, tf, N), with N=__grid_size_direct_shooting() -u_init(t) = [u_sol(t)-1.0] -init_, grid_ = CTOptimizationInit(t0, tf, m, u_init, nothing, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(grid_[i])-1.0] for i = 1:__grid_size_direct_shooting()-1]) atol=1e-4 -@test grid_ == __grid(t0, tf) - -# init=u, grid=T => init=u(T), grid=T (check validity with ocp) -u_init(t) = [u_sol(t)-1.0] -N = floor(Int64, __grid_size_direct_shooting()/2) -T = range(t0, tf, N) -init_, grid_ = CTOptimizationInit(t0, tf, m, u_init, T, __init_interpolation()) -@test init_ ≈ convert_init([[u_sol(T[i])-1.0] for i = 1:N-1]) atol=1e-4 -@test grid_ == T - -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, u_init, range(t0-1, tf, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, u_init, range(t0, tf+1, N), __init_interpolation()) -@test_throws InconsistentArgument CTOptimizationInit(t0, tf, m, u_init, range(tf, t0, N), __init_interpolation()) - -# -------------------------------------------------------------------------------------------------- -# resolution -# -# init, grid -N = floor(Int64, __grid_size_direct_shooting()/2) -T = range(t0, tf, N) -U = [[u_sol(T[i])-1.0] for i = 1:N-1] - -N_ = floor(Int64, __grid_size_direct_shooting()/4) -T_ = range(t0, tf, N_) - -u_init(t) = [u_sol(t)-1.0] - -# resolution with different init -common_args = (iterations=5, display=false) -sol = solve(ocp, :descent, init=nothing, grid=nothing; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=nothing, grid=T; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=U, grid=nothing; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=U, grid=T; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=(T,U), grid=nothing; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=(T,U), grid=T_; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=sol, grid=nothing; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=u_init, grid=nothing; common_args...); @test typeof(sol) == DirectShootingSolution; -sol = solve(ocp, :descent, init=u_init, grid=T; common_args...); @test typeof(sol) == DirectShootingSolution; - -# plots -@test typeof(plot(sol)) == Plots.Plot{Plots.GRBackend} -@test typeof(plot(sol, :time, (:control, 1))) == Plots.Plot{Plots.GRBackend} diff --git a/test/test_goddard_direct.jl b/test/test_goddard_direct.jl index 494439ab6..736e318b1 100644 --- a/test/test_goddard_direct.jl +++ b/test/test_goddard_direct.jl @@ -1,51 +1,6 @@ -# Parameters -Cd = 310. -Tmax = 3.5 -β = 500. -b = 2. -N = 100 -t0 = 0. -r0 = 1. -v0 = 0. -vmax = 0.1 -m0 = 1. -mf = 0.6 -x0 = [ r0, v0, m0 ] - -epsilon = 1.e-2 - -ocp = Model() - -time!(ocp, :initial, t0) # if not provided, final time is free -state!(ocp, 3) # state dim -control!(ocp, 1) # control dim - -constraint!(ocp, :initial, x0) # vectorial equality -constraint!(ocp, :control, u -> u[1], 0., 1.) # constraints can be labeled or not -constraint!(ocp, :mixed, (x, u) -> x[1], r0, Inf, :state_con1) # epsilon car sinon adjoint state at 0 false because the state is r0 touch the constraint at 0 -constraint!(ocp, :mixed, (x, u) -> x[2], 0., vmax, :state_con2) -constraint!(ocp, :mixed, (x, u) -> x[3], mf, m0, :state_con3) - -objective!(ocp, :mayer, (t0, x0, tf, xf) -> xf[1], :max) - -function F0(x) - r, v, m = x - D = Cd * v^2 * exp(-β*(r-1.)) - F = [ v, -D/m-1.0/r^2, 0. ] - return F -end - -function F1(x) - r, v, m = x - F = [ 0., Tmax/m, -b*Tmax ] - return F -end - -function f(x, u) - return F0(x) + u[1]*F1(x) -end - -constraint!(ocp, :dynamics, f) +# goddard with state constraint - maximize altitude +prob = Problem(:goddard, :state_constraint) +ocp = prob.model init = [1.01, 0.25, 0.5, 0.4] sol = solve(ocp, grid_size=10, print_level=0, init=init) diff --git a/test/test_goddard_indirect.jl b/test/test_goddard_indirect.jl index e954f843c..091bc3fdb 100644 --- a/test/test_goddard_indirect.jl +++ b/test/test_goddard_indirect.jl @@ -39,9 +39,6 @@ constraint!(ocp, :dynamics, f) @test typeof(ocp) == OptimalControlModel{:autonomous} @test ocp.initial_time == t0 -# -------------------------------------------------------- -# Direct - # -------------------------------------------------------- # Indirect diff --git a/test/test_utils.jl b/test/test_utils.jl deleted file mode 100644 index 9ef5dc6df..000000000 --- a/test/test_utils.jl +++ /dev/null @@ -1,5 +0,0 @@ -v = [1.0; 2.0; 3.0; 4.0; 5.0; 6.0] -n = 2 -u = vec2vec(v, n) -w = vec2vec(u) -@test v == w \ No newline at end of file From ea52bc4ddb5e6100784ccf74c7473d9867d0d0b2 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Sun, 19 Feb 2023 15:44:07 +0100 Subject: [PATCH 7/7] rm ControlToolboxTools --- Project.toml | 12 +----------- src/OptimalControl.jl | 6 ------ test/Project.toml | 2 -- 3 files changed, 1 insertion(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 9ddfb730c..732ea8aee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,23 +1,13 @@ name = "OptimalControl" uuid = "5f98b655-cc9a-415a-b60e-744165666948" authors = ["Olivier Cots "] -version = "0.2.0" +version = "0.2.1" [deps] -ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd" CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5" CTDirectShooting = "e0e6c04b-5022-4cd2-bea2-4a09fff39444" -CTOptimization = "22f08de8-270f-4470-8fba-397dbc90d8e0" -ControlToolboxTools = "3a0bcf43-9180-47f3-913c-e71e0d69f39f" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HamiltonianFlows = "5fb78580-10a0-4606-82bc-07a60f425ab3" -Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" -Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] diff --git a/src/OptimalControl.jl b/src/OptimalControl.jl index cfd04b8b2..b5f7906f8 100644 --- a/src/OptimalControl.jl +++ b/src/OptimalControl.jl @@ -8,16 +8,10 @@ using Reexport using CTDirect using CTDirectShooting -# tools: callbacks, exceptions, functions and more -@reexport using ControlToolboxTools - # flows @reexport using HamiltonianFlows import HamiltonianFlows: Flow -# Types -const MyNumber, MyVector, Time, Times, TimesDisc, States, Adjoints, Controls, State, Adjoint, Dimension = CTBase.types() - # Other declarations const __display = CTBase.__display diff --git a/test/Project.toml b/test/Project.toml index 91a5d9a08..98ab2e8de 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,5 @@ [deps] CTProblemLibrary = "0649932a-8c77-4f67-b1e4-c19ddd080280" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"