From af7dbe3a084d19408de3730a5ce251bcfb1e7de8 Mon Sep 17 00:00:00 2001 From: Olivier Cots <66357348+ocots@users.noreply.github.com> Date: Mon, 26 Aug 2024 09:16:54 +0200 Subject: [PATCH] format code (#218) * Update README.md * format code * fix merge --------- Co-authored-by: Olivier Cots --- .github/workflows/Formatter.yml | 33 ++ README.md | 9 - benchmark/benchmark.jl | 83 +++-- benchmark/prof.jl | 72 +++-- docs/make.jl | 13 +- ext/CTDirectExt.jl | 19 +- ext/CTSolveExtIpopt.jl | 51 +++- ext/CTSolveExtMadNLP.jl | 40 ++- src/CTDirect.jl | 2 +- src/default.jl | 2 +- src/problem.jl | 179 +++++++---- src/solution.jl | 322 ++++++++++++++------ src/solve.jl | 71 +++-- src/utils.jl | 38 +-- test/manual_test.jl | 2 +- test/old/bocop_goddard.jl | 43 ++- test/old/local_test_exponential.jl | 22 +- test/old/manual_simple_integrator_energy.jl | 4 +- test/old/manual_test_all_constraints.jl | 151 +++++++-- test/old/manual_test_exponential.jl | 14 +- test/old/manual_test_integrator.jl | 10 +- test/old/test_goddard.jl | 36 +-- test/old/test_integrator.jl | 42 +-- test/old/test_integrator_constraints.jl | 195 ++++++------ test/old/test_scalar_vector_CTBase.jl | 97 +++--- test/old/test_simple_integrator_energy.jl | 10 +- test/old/test_utils.jl | 10 +- test/problems/beam.jl | 12 +- test/problems/bioreactor.jl | 49 +-- test/problems/bolza.jl | 3 +- test/problems/double_integrator.jl | 62 ++-- test/problems/fuller.jl | 12 +- test/problems/goddard.jl | 95 +++--- test/problems/insurance.jl | 26 +- test/problems/jackson.jl | 18 +- test/problems/parametric.jl | 18 +- test/problems/robbins.jl | 14 +- test/problems/simple_integrator.jl | 14 +- test/problems/swimmer.jl | 129 ++++++-- test/problems/vanderpol.jl | 14 +- test/refine_grid.jl | 130 +++++--- test/runtests.jl | 2 +- test/scalar_vector.jl | 14 +- test/suite/test_abstract_ocp.jl | 8 +- test/suite/test_constraints.jl | 112 ++++--- test/suite/test_continuation.jl | 31 +- test/suite/test_grid.jl | 36 ++- test/suite/test_initial_guess.jl | 136 +++++++-- test/suite/test_misc.jl | 8 +- test/suite/test_nlp.jl | 50 +-- test/suite/test_objective.jl | 20 +- 51 files changed, 1642 insertions(+), 941 deletions(-) create mode 100644 .github/workflows/Formatter.yml diff --git a/.github/workflows/Formatter.yml b/.github/workflows/Formatter.yml new file mode 100644 index 00000000..6f666abe --- /dev/null +++ b/.github/workflows/Formatter.yml @@ -0,0 +1,33 @@ +name: Formatter + +# Modified from https://github.com/julia-actions/julia-format +on: + schedule: + - cron: '0 0 * * *' + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/cache@v2 + - name: Install JuliaFormatter and format + run: | + julia -e 'import Pkg; Pkg.add("JuliaFormatter")' + julia -e 'using JuliaFormatter; format(".")' + # https://github.com/marketplace/actions/create-pull-request + # https://github.com/peter-evans/create-pull-request#reference-example + - name: Create Pull Request + id: cpr + uses: peter-evans/create-pull-request@v3 + with: + token: ${{ secrets.GITHUB_TOKEN }} + commit-message: ":robot: Format .jl files" + title: '[AUTO] JuliaFormatter.jl run' + branch: auto-juliaformatter-pr + delete-branch: true + labels: formatting, automated pr, no changelog + - name: Check outputs + run: | + echo "Pull Request Number - ${{ steps.cpr.outputs.pull-request-number }}" + echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}" \ No newline at end of file diff --git a/README.md b/README.md index d21a5bd3..d805a0fb 100644 --- a/README.md +++ b/README.md @@ -50,15 +50,6 @@ pkg> add CTDirect ## Contributing -If you think you found a bug or if you have a feature request or suggestion, feel free to open an [issue](https://github.com/control-toolbox/CTDirect.jl/issues). -Before opening a pull request, start an issue or a discussion on the topic, please. - -Any contributions are welcomed, check out [how to contribute to a Github project](https://docs.github.com/en/get-started/exploring-projects-on-github/contributing-to-a-project). First good issues, if any, may be found [here](https://github.com/control-toolbox/CTDirect.jl/contribute). - -If you want to ask a question, feel free to start a discussion [here](https://github.com/orgs/control-toolbox/discussions). This forum is for general discussion about this repository and the [control-toolbox organization](https://github.com/control-toolbox), so questions about any of our packages are welcome. - -## Contributing - [issue-url]: https://github.com/control-toolbox/CTDirect.jl/issues [first-good-issue-url]: https://github.com/control-toolbox/CTDirect.jl/contribute diff --git a/benchmark/benchmark.jl b/benchmark/benchmark.jl index 986e528d..5ed06f1d 100644 --- a/benchmark/benchmark.jl +++ b/benchmark/benchmark.jl @@ -8,14 +8,21 @@ using MadNLPMumps ####################################################### # load examples library -problem_path = pwd()*"/test/problems" -for problem_file in filter(contains(r".jl$"), readdir(problem_path; join=true)) +problem_path = pwd() * "/test/problems" +for problem_file in filter(contains(r".jl$"), readdir(problem_path; join = true)) include(problem_file) end -function bench(;nlp_solver = :ipopt, linear_solver = nothing, - tol=1e-8, grid_size=1000, precompile = true, display=false, verbose=true) +function bench(; + nlp_solver = :ipopt, + linear_solver = nothing, + tol = 1e-8, + grid_size = 1000, + precompile = true, + display = false, + verbose = true, +) ####################################################### # set (non) linear solvers and backends @@ -30,45 +37,77 @@ function bench(;nlp_solver = :ipopt, linear_solver = nothing, linear_solver = "UmfpackSolver" end - verbose && @printf("Profile: NLP Solver %s with linear solver %s\n", nlp_solver, linear_solver) + verbose && + @printf("Profile: NLP Solver %s with linear solver %s\n", nlp_solver, linear_solver) # blas backend (cf using MKL above, should be option...) verbose && @printf("Blas config: %s\n", LinearAlgebra.BLAS.lbt_get_config()) # settings - verbose && @printf("Settings: tol=%g grid_size=%d precompile=%s\n\n", tol, grid_size, precompile) + verbose && @printf( + "Settings: tol=%g grid_size=%d precompile=%s\n\n", + tol, + grid_size, + precompile + ) # load problems for benchmark names_list = ["beam", "bioreactor_1day", "fuller", "goddard", "jackson", "vanderpol"] problem_list = [] for problem_name in names_list ocp_data = getfield(Main, Symbol(problem_name))() - push!(problem_list,ocp_data) + push!(problem_list, ocp_data) end ####################################################### # precompile if required if precompile - t_precomp = 0. + t_precomp = 0.0 verbose && print("Precompilation: ") for problem in problem_list - verbose && @printf("%s ",problem[:name]) - t = @elapsed direct_solve(problem[:ocp], nlp_solver, linear_solver=linear_solver, max_iter=0, display=display) + verbose && @printf("%s ", problem[:name]) + t = @elapsed direct_solve( + problem[:ocp], + nlp_solver, + linear_solver = linear_solver, + max_iter = 0, + display = display, + ) t_precomp += t end - verbose && @printf("\nPrecompilation total time %6.2f\n\n",t_precomp) + verbose && @printf("\nPrecompilation total time %6.2f\n\n", t_precomp) end ####################################################### # solve examples with timer and objective check t_list = [] for problem in problem_list - t = @elapsed sol = direct_solve(problem[:ocp], nlp_solver, init=problem[:init], display=display, linear_solver=linear_solver, grid_size=grid_size, tol=tol) - if !isnothing(problem[:obj]) && !isapprox(sol.objective, problem[:obj], rtol=5e-2) - error("Objective mismatch for ", problem[:name], ": ", sol.objective, " instead of ", problem[:obj]) + t = @elapsed sol = direct_solve( + problem[:ocp], + nlp_solver, + init = problem[:init], + display = display, + linear_solver = linear_solver, + grid_size = grid_size, + tol = tol, + ) + if !isnothing(problem[:obj]) && !isapprox(sol.objective, problem[:obj], rtol = 5e-2) + error( + "Objective mismatch for ", + problem[:name], + ": ", + sol.objective, + " instead of ", + problem[:obj], + ) else - verbose && @printf("%-30s completed in %6.2f s after %4d iterations\n",problem[:name],t,sol.iterations) - append!(t_list,t) + verbose && @printf( + "%-30s completed in %6.2f s after %4d iterations\n", + problem[:name], + t, + sol.iterations + ) + append!(t_list, t) end end @@ -83,16 +122,16 @@ function bench(;nlp_solver = :ipopt, linear_solver = nothing, end # +++ put repeat directly in bench() -function bench_average(; repeat=2, verbose=false, kwargs...) +function bench_average(; repeat = 2, verbose = false, kwargs...) # execute series of benchmark runs t_list = [] - for i in 1:repeat - t = bench(;verbose=verbose, kwargs...) + for i = 1:repeat + t = bench(; verbose = verbose, kwargs...) append!(t_list, t) @printf("Run %d / %d: time (s) = %6.2f\n", i, repeat, t) end - + # print / return average total time avg_time = sum(t_list) / length(t_list) @printf("Average time (s): %6.2f\n", avg_time) @@ -101,11 +140,11 @@ function bench_average(; repeat=2, verbose=false, kwargs...) end -function bench_series(;grid_size_list=[250, 500, 1000, 2500, 5000, 10000], kwargs...) +function bench_series(; grid_size_list = [250, 500, 1000, 2500, 5000, 10000], kwargs...) println(grid_size_list) t_list = [] for grid_size in grid_size_list - t = bench_average(;grid_size=grid_size, kwargs...) + t = bench_average(; grid_size = grid_size, kwargs...) append!(t_list, t) @printf("Grid size %d: time (s) = %6.2f\n\n", grid_size, t) end diff --git a/benchmark/prof.jl b/benchmark/prof.jl index fb39c58a..67c5860f 100644 --- a/benchmark/prof.jl +++ b/benchmark/prof.jl @@ -21,48 +21,60 @@ prob = include("../problems/fuller.jl") #prob = include("../problems/goddard.jl") ocp = prob[:ocp] grid_size = 100 -docp, nlp = direct_transcription(ocp, grid_size=grid_size) +docp, nlp = direct_transcription(ocp, grid_size = grid_size) println("Load problem ", prob[:name]) # full solve if test_time - if precompile - println("Precompilation") - solve(ocp, grid_size=grid_size, display=false, max_iter=2) - end - println("Timed solve") - @timev sol = solve(ocp, grid_size=grid_size, print_level=0) + if precompile + println("Precompilation") + solve(ocp, grid_size = grid_size, display = false, max_iter = 2) + end + println("Timed solve") + @timev sol = solve(ocp, grid_size = grid_size, print_level = 0) end if precompile - println("Precompilation") - if test == :objective - CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) - else - CTDirect.DOCP_constraints!(zeros(docp.dim_NLP_constraints), CTDirect.DOCP_initial_guess(docp), docp) - end + println("Precompilation") + if test == :objective + CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) + else + CTDirect.DOCP_constraints!( + zeros(docp.dim_NLP_constraints), + CTDirect.DOCP_initial_guess(docp), + docp, + ) + end end if test_code_warntype - if test == :objective - # NB. Pb with the mayer part: obj is type unstable (Any) because ocp.mayer is Union(Mayer,nothing), even for mayer problems (also, we should not even enter this code part for lagrange problems since has_mayer us defined as const in DOCP oO ...). - @code_warntype CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) - else - # OK ! - @code_warntype CTDirect.DOCP_constraints!(zeros(docp.dim_NLP_constraints), CTDirect.DOCP_initial_guess(docp), docp) - end + if test == :objective + # NB. Pb with the mayer part: obj is type unstable (Any) because ocp.mayer is Union(Mayer,nothing), even for mayer problems (also, we should not even enter this code part for lagrange problems since has_mayer us defined as const in DOCP oO ...). + @code_warntype CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) + else + # OK ! + @code_warntype CTDirect.DOCP_constraints!( + zeros(docp.dim_NLP_constraints), + CTDirect.DOCP_initial_guess(docp), + docp, + ) + end end if test_jet - if test == :objective - # 4 possible errors - # due to the ocp.mayer type problem cf above - @report_opt CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) - else - # 50 possible errors: some getindex (Integer vs Int...) - # all variables x,u,v - @report_opt CTDirect.DOCP_constraints!(zeros(docp.dim_NLP_constraints), CTDirect.DOCP_initial_guess(docp), docp) - end + if test == :objective + # 4 possible errors + # due to the ocp.mayer type problem cf above + @report_opt CTDirect.DOCP_objective(CTDirect.DOCP_initial_guess(docp), docp) + else + # 50 possible errors: some getindex (Integer vs Int...) + # all variables x,u,v + @report_opt CTDirect.DOCP_constraints!( + zeros(docp.dim_NLP_constraints), + CTDirect.DOCP_initial_guess(docp), + docp, + ) + end end #= @@ -73,4 +85,4 @@ if test_ipopt println("\n Discrete solution") println(dsol) end -=# \ No newline at end of file +=# diff --git a/docs/make.jl b/docs/make.jl index 91c6d222..9ceb1e35 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,19 +19,12 @@ makedocs( format = Documenter.HTML( prettyurls = false, size_threshold_ignore = ["api-ctbase.md"], - assets=[ + assets = [ asset("https://control-toolbox.org/assets/css/documentation.css"), asset("https://control-toolbox.org/assets/js/documentation.js"), ], ), - pages = [ - "Introduction" => "index.md", - "API" => "api.md", - "Developers" => "dev-api.md", - ] + pages = ["Introduction" => "index.md", "API" => "api.md", "Developers" => "dev-api.md"], ) -deploydocs( - repo = "github.com/control-toolbox/CTDirect.jl.git", - devbranch = "main" -) +deploydocs(repo = "github.com/control-toolbox/CTDirect.jl.git", devbranch = "main") diff --git a/ext/CTDirectExt.jl b/ext/CTDirectExt.jl index bd8177b3..b1b08d69 100644 --- a/ext/CTDirectExt.jl +++ b/ext/CTDirectExt.jl @@ -13,17 +13,17 @@ $(TYPEDSIGNATURES) Save OCP solution in JLD2 format """ -function JLD2.save(sol::OptimalControlSolution; filename_prefix="solution") +function JLD2.save(sol::OptimalControlSolution; filename_prefix = "solution") save_object(filename_prefix * ".jld2", sol) return nothing end - + """ $(TYPEDSIGNATURES) Load OCP solution in JLD2 format """ -function JLD2.load(filename_prefix="solution") +function JLD2.load(filename_prefix = "solution") return load_object(filename_prefix * ".jld2") end @@ -32,8 +32,11 @@ $(TYPEDSIGNATURES) Export OCP solution in JSON format """ -function CTDirect.export_ocp_solution(sol::OptimalControlSolution; filename_prefix="solution") -# +++ redo this, start with basics +function CTDirect.export_ocp_solution( + sol::OptimalControlSolution; + filename_prefix = "solution", +) + # +++ redo this, start with basics open(filename_prefix * ".json", "w") do io #JSON3.pretty(io, CTDirect.OCPDiscreteSolution(sol)) end @@ -45,11 +48,11 @@ $(TYPEDSIGNATURES) Read OCP solution in JSON format """ -function CTDirect.import_ocp_solution(filename_prefix="solution") -# +++ add constructor from json blob +function CTDirect.import_ocp_solution(filename_prefix = "solution") + # +++ add constructor from json blob json_string = read(filename_prefix * ".json", String) return OptimalControlSolution(JSON3.read(json_string)) end -end \ No newline at end of file +end diff --git a/ext/CTSolveExtIpopt.jl b/ext/CTSolveExtIpopt.jl index e19f61d1..088f6a9c 100644 --- a/ext/CTSolveExtIpopt.jl +++ b/ext/CTSolveExtIpopt.jl @@ -12,41 +12,60 @@ $(TYPEDSIGNATURES) Solve a discretized optimal control problem with Ipopt """ -function CTDirect.solve_docp(tag::CTDirect.IpoptTag, docp::CTDirect.DOCP, nlp; - display::Bool=CTBase.__display(), - max_iter::Integer=CTDirect.__max_iterations(), - tol::Real=CTDirect.__tolerance(), - print_level::Integer=CTDirect.__ipopt_print_level(), - mu_strategy::String=CTDirect.__ipopt_mu_strategy(), - linear_solver::String=CTDirect.__ipopt_linear_solver(), - kwargs...) +function CTDirect.solve_docp( + tag::CTDirect.IpoptTag, + docp::CTDirect.DOCP, + nlp; + display::Bool = CTBase.__display(), + max_iter::Integer = CTDirect.__max_iterations(), + tol::Real = CTDirect.__tolerance(), + print_level::Integer = CTDirect.__ipopt_print_level(), + mu_strategy::String = CTDirect.__ipopt_mu_strategy(), + linear_solver::String = CTDirect.__ipopt_linear_solver(), + kwargs..., +) # check HSL requirements if linear_solver in ["ma27", "ma57", "ma77", "ma86", "ma97"] && !LIBHSL_isfunctional() - println("WARNING: HSL not available, defaulting given linear solver ",linear_solver, " to MUMPS") + println( + "WARNING: HSL not available, defaulting given linear solver ", + linear_solver, + " to MUMPS", + ) linear_solver = "mumps" end # check SPRAL requirements - if linear_solver == "spral" && (!haskey(ENV, "OMP_CANCELLATION") || !haskey(ENV, "OMP_PROC_BIND")) - println("WARNING: missing required environment variables for SPRAL (OMP_CANCELLATION=TRUE and OMP_PROC_BIND=TRUE), defaulting to MUMPS") + if linear_solver == "spral" && + (!haskey(ENV, "OMP_CANCELLATION") || !haskey(ENV, "OMP_PROC_BIND")) + println( + "WARNING: missing required environment variables for SPRAL (OMP_CANCELLATION=TRUE and OMP_PROC_BIND=TRUE), defaulting to MUMPS", + ) linear_solver = "mumps" end # disable output if needed - print_level = display ? print_level : 0 + print_level = display ? print_level : 0 # preallocate solver solver = IpoptSolver(nlp) # solve discretized problem with NLP solver - docp_solution = solve!(solver, nlp, print_level=print_level, mu_strategy=mu_strategy, tol=tol, max_iter=max_iter, sb="yes", linear_solver=linear_solver; kwargs...) + docp_solution = solve!( + solver, + nlp, + print_level = print_level, + mu_strategy = mu_strategy, + tol = tol, + max_iter = max_iter, + sb = "yes", + linear_solver = linear_solver; + kwargs..., + ) # return DOCP solution return docp_solution - -end end - +end diff --git a/ext/CTSolveExtMadNLP.jl b/ext/CTSolveExtMadNLP.jl index ce1568c0..83ae7525 100644 --- a/ext/CTSolveExtMadNLP.jl +++ b/ext/CTSolveExtMadNLP.jl @@ -12,21 +12,25 @@ $(TYPEDSIGNATURES) Solve a discretized optimal control problem DOCP """ -function CTDirect.solve_docp(tag::CTDirect.MadNLPTag, docp::CTDirect.DOCP, nlp; - display::Bool=CTBase.__display(), - max_iter::Integer=CTDirect.__max_iterations(), - tol::Real=CTDirect.__tolerance(), - linear_solver::String=CTDirect.__madnlp_linear_solver(), - kwargs...) +function CTDirect.solve_docp( + tag::CTDirect.MadNLPTag, + docp::CTDirect.DOCP, + nlp; + display::Bool = CTBase.__display(), + max_iter::Integer = CTDirect.__max_iterations(), + tol::Real = CTDirect.__tolerance(), + linear_solver::String = CTDirect.__madnlp_linear_solver(), + kwargs..., +) # disable output if needed - print_level = display ? MadNLP.INFO : MadNLP.ERROR + print_level = display ? MadNLP.INFO : MadNLP.ERROR # preallocate solver (NB. need to pass printlevel here) - solver = MadNLPSolver(nlp, print_level=print_level) + solver = MadNLPSolver(nlp, print_level = print_level) # solve discretized problem with NLP solver - docp_solution = solve!(solver, tol=tol, max_iter=max_iter; kwargs...) + docp_solution = solve!(solver, tol = tol, max_iter = max_iter; kwargs...) # return DOCP solution return docp_solution @@ -38,16 +42,24 @@ function CTBase.OptimalControlSolution(docp, docp_solution::MadNLP.MadNLPExecuti # adjust objective sign for maximization problems if is_min(docp.ocp) objective = docp_solution.objective - else - objective = - docp_solution.objective + else + objective = -docp_solution.objective end # call lower level constructor - return OptimalControlSolution(docp, primal=docp_solution.solution, dual=docp_solution.multipliers, objective=objective, iterations=docp_solution.iter, constraints_violation=docp_solution.primal_feas, message="MadNLP", mult_LB=docp_solution.multipliers_L, mult_UB=docp_solution.multipliers_U) + return OptimalControlSolution( + docp, + primal = docp_solution.solution, + dual = docp_solution.multipliers, + objective = objective, + iterations = docp_solution.iter, + constraints_violation = docp_solution.primal_feas, + message = "MadNLP", + mult_LB = docp_solution.multipliers_L, + mult_UB = docp_solution.multipliers_U, + ) end end - - diff --git a/src/CTDirect.jl b/src/CTDirect.jl index 3dfe02d3..74e17ab4 100644 --- a/src/CTDirect.jl +++ b/src/CTDirect.jl @@ -28,4 +28,4 @@ export export_ocp_solution export import_ocp_solution export direct_solve -end \ No newline at end of file +end diff --git a/src/default.jl b/src/default.jl index 21f39a36..83d68064 100644 --- a/src/default.jl +++ b/src/default.jl @@ -68,4 +68,4 @@ $(TYPEDSIGNATURES) Used to set the default value of the linear solver of MadNLP for the direct method. The default value is `umfpack`. """ -__madnlp_linear_solver() = "umfpack" \ No newline at end of file +__madnlp_linear_solver() = "umfpack" diff --git a/src/problem.jl b/src/problem.jl index c2e3278e..6d9e030e 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -15,14 +15,14 @@ struct DOCP ## OCP ocp::OptimalControlModel - control_constraints - state_constraints - mixed_constraints - boundary_constraints - variable_constraints - control_box - state_box - variable_box + control_constraints::Any + state_constraints::Any + mixed_constraints::Any + boundary_constraints::Any + variable_constraints::Any + control_box::Any + state_box::Any + variable_box::Any # faster than calling ocp functions each time ? has_free_t0::Bool @@ -59,16 +59,20 @@ struct DOCP con_u::Vector{Float64} # constructor - function DOCP(ocp::OptimalControlModel, grid_size::Integer, time_grid) + function DOCP(ocp::OptimalControlModel, grid_size::Integer, time_grid) # time grid if time_grid == nothing - NLP_normalized_time_grid = collect(LinRange(0, 1, grid_size+1)) + NLP_normalized_time_grid = collect(LinRange(0, 1, grid_size + 1)) dim_NLP_steps = grid_size else # check strictly increasing - if !issorted(time_grid,lt=<=) - throw(ArgumentError("given time grid is not strictly increasing. Aborting...")) + if !issorted(time_grid, lt = <=) + throw( + ArgumentError( + "given time grid is not strictly increasing. Aborting...", + ), + ) return nothing end # normalize input grid if needed @@ -76,7 +80,7 @@ struct DOCP #println("INFO: normalizing given time grid...") t0 = time_grid[1] tf = time_grid[end] - time_grid = (time_grid .- t0) ./ (tf - t0) + time_grid = (time_grid .- t0) ./ (tf - t0) end NLP_normalized_time_grid = time_grid dim_NLP_steps = length(time_grid) - 1 @@ -111,7 +115,14 @@ struct DOCP # NLP constraints # parse NLP constraints (and initialize dimensions) - control_constraints, state_constraints, mixed_constraints, boundary_constraints, variable_constraints, control_box, state_box, variable_box = CTBase.nlp_constraints!(ocp) + control_constraints, + state_constraints, + mixed_constraints, + boundary_constraints, + variable_constraints, + control_box, + state_box, + variable_box = CTBase.nlp_constraints!(ocp) dim_x_box = dim_state_range(ocp) dim_u_box = dim_control_range(ocp) @@ -124,14 +135,51 @@ struct DOCP dim_boundary_cons = dim_boundary_constraints(ocp) # lagrange to mayer transformation - dim_NLP_constraints = N * (dim_NLP_x + dim_path_cons) + dim_path_cons + dim_boundary_cons + dim_v_cons + dim_NLP_constraints = + N * (dim_NLP_x + dim_path_cons) + dim_path_cons + dim_boundary_cons + dim_v_cons if has_lagrange - dim_NLP_constraints += 1 + dim_NLP_constraints += 1 end # call constructor with const fields - docp = new(ocp, control_constraints, state_constraints, mixed_constraints, boundary_constraints, variable_constraints, control_box, state_box, variable_box, has_free_t0, has_free_tf, has_lagrange, has_mayer, has_variable, - has_maximization, dim_x_box, dim_u_box, dim_v_box, dim_path_cons, dim_x_cons, dim_u_cons, dim_v_cons,dim_mixed_cons, dim_boundary_cons, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_OCP_x, dim_NLP_steps, NLP_normalized_time_grid, dim_NLP_variables, dim_NLP_constraints,-Inf * ones(dim_NLP_variables), Inf * ones(dim_NLP_variables), zeros(dim_NLP_constraints), zeros(dim_NLP_constraints)) + docp = new( + ocp, + control_constraints, + state_constraints, + mixed_constraints, + boundary_constraints, + variable_constraints, + control_box, + state_box, + variable_box, + has_free_t0, + has_free_tf, + has_lagrange, + has_mayer, + has_variable, + has_maximization, + dim_x_box, + dim_u_box, + dim_v_box, + dim_path_cons, + dim_x_cons, + dim_u_cons, + dim_v_cons, + dim_mixed_cons, + dim_boundary_cons, + dim_NLP_x, + dim_NLP_u, + dim_NLP_v, + dim_OCP_x, + dim_NLP_steps, + NLP_normalized_time_grid, + dim_NLP_variables, + dim_NLP_constraints, + -Inf * ones(dim_NLP_variables), + Inf * ones(dim_NLP_variables), + zeros(dim_NLP_constraints), + zeros(dim_NLP_constraints), + ) return docp @@ -162,18 +210,18 @@ function constraints_bounds!(docp::DOCP) ub = docp.con_u index = 1 # counter for the constraints - for i in 0:docp.dim_NLP_steps-1 + for i = 0:docp.dim_NLP_steps-1 # skip (ie leave 0) for equality dynamics constraint index = index + docp.dim_NLP_x # path constraints index = setPathBounds!(docp, index, lb, ub) end - + # path constraints at final time - index = setPathBounds!(docp, index, lb, ub) + index = setPathBounds!(docp, index, lb, ub) # boundary and variable constraints - index = setPointBounds!(docp, index, lb, ub) + index = setPointBounds!(docp, index, lb, ub) return lb, ub end @@ -196,29 +244,29 @@ function variables_bounds!(docp::DOCP) # build ordered bounds vectors for state and control x_lb = -Inf * ones(docp.dim_OCP_x) x_ub = Inf * ones(docp.dim_OCP_x) - for j in 1:docp.dim_x_box + for j = 1:docp.dim_x_box indice = docp.state_box[2][j] x_lb[indice] = docp.state_box[1][j] x_ub[indice] = docp.state_box[3][j] end u_lb = -Inf * ones(docp.dim_NLP_u) u_ub = Inf * ones(docp.dim_NLP_u) - for j in 1:docp.dim_u_box + for j = 1:docp.dim_u_box indice = docp.control_box[2][j] u_lb[indice] = docp.control_box[1][j] u_ub[indice] = docp.control_box[3][j] end # apply bounds for NLP variables - for i in 0:N + for i = 0:N set_variables_at_time_step!(var_l, x_lb, u_lb, docp, i) set_variables_at_time_step!(var_u, x_ub, u_ub, docp, i) end # variable box - offset = (N+1) * (docp.dim_NLP_x + docp.dim_NLP_u) + offset = (N + 1) * (docp.dim_NLP_x + docp.dim_NLP_u) if docp.dim_v_box > 0 - for j in 1:docp.dim_v_box + for j = 1:docp.dim_v_box indice = docp.variable_box[2][j] var_l[offset+indice] = docp.variable_box[1][j] var_u[offset+indice] = docp.variable_box[3][j] @@ -237,7 +285,7 @@ Compute the objective for the DOCP problem. # DOCP objective function DOCP_objective(xu, docp::DOCP) - obj = 0. + obj = 0.0 N = docp.dim_NLP_steps ocp = docp.ocp @@ -246,23 +294,23 @@ function DOCP_objective(xu, docp::DOCP) v = get_variable(xu, docp) t0 = get_initial_time(xu, docp) tf = get_final_time(xu, docp) - x0,u0,xl0 = get_variables_at_time_step(xu, docp, 0) - xf,uf,xlf = get_variables_at_time_step(xu, docp, N) + x0, u0, xl0 = get_variables_at_time_step(xu, docp, 0) + xf, uf, xlf = get_variables_at_time_step(xu, docp, N) #obj = obj + ocp.mayer(t0, tf, x0, xf, v) obj = obj + ocp.mayer(x0, xf, v) end - + # lagrange cost if docp.has_lagrange - xf,uf,xlf = get_variables_at_time_step(xu, docp, N) + xf, uf, xlf = get_variables_at_time_step(xu, docp, N) obj = obj + xlf end # maximization problem if docp.has_maximization obj = -obj - end - + end + return obj end @@ -272,7 +320,7 @@ $(TYPEDSIGNATURES) Compute the constraints C for the DOCP problem (modeled as LB <= C(X) <= UB). """ -function DOCP_constraints!(c, xu, docp::DOCP) +function DOCP_constraints!(c, xu, docp::DOCP) """ compute the constraints for the NLP : - discretization of the dynamics via the trapeze method @@ -292,7 +340,7 @@ function DOCP_constraints!(c, xu, docp::DOCP) args_ip1 = ArgsAtTimeStep(xu, docp, 1, v) index = 1 # counter for the constraints - for i in 0:docp.dim_NLP_steps-1 + for i = 0:docp.dim_NLP_steps-1 # state equation index = setStateEquation!(docp, c, index, (args_i, args_ip1)) @@ -300,9 +348,9 @@ function DOCP_constraints!(c, xu, docp::DOCP) index = setPathConstraints!(docp, c, index, args_i, v) # smart update for next iteration - if i < docp.dim_NLP_steps-1 + if i < docp.dim_NLP_steps - 1 args_i = args_ip1 - args_ip1 = ArgsAtTimeStep(xu, docp, i+2, v) + args_ip1 = ArgsAtTimeStep(xu, docp, i + 2, v) end end @@ -316,7 +364,7 @@ function DOCP_constraints!(c, xu, docp::DOCP) # needed even for inplace version, AD error otherwise # may be because actual return would be index above ? - return c + return c end @@ -330,12 +378,12 @@ Useful values at a time step: time, state, control, dynamics... """ struct ArgsAtTimeStep - time - state - control - dynamics - lagrange_state - lagrange_cost + time::Any + state::Any + control::Any + dynamics::Any + lagrange_state::Any + lagrange_cost::Any function ArgsAtTimeStep(xu, docp::DOCP, i::Int, v) @@ -368,19 +416,25 @@ $(TYPEDSIGNATURES) Set the constraints corresponding to the state equation """ function setStateEquation!(docp::DOCP, c, index::Int, args_trapeze) - + ocp = docp.ocp args_i = args_trapeze[1] args_ip1 = args_trapeze[2] hi = args_ip1.time - args_i.time - + # trapeze rule - c[index:index+docp.dim_OCP_x-1] .= args_ip1.state .- (args_i.state .+ 0.5*hi*(args_i.dynamics .+ args_ip1.dynamics)) + c[index:index+docp.dim_OCP_x-1] .= + args_ip1.state .- + (args_i.state .+ 0.5 * hi * (args_i.dynamics .+ args_ip1.dynamics)) if docp.has_lagrange - c[index+docp.dim_OCP_x] = args_ip1.lagrange_state - (args_i.lagrange_state + 0.5*hi*(args_i.lagrange_cost + args_ip1.lagrange_cost)) + c[index+docp.dim_OCP_x] = + args_ip1.lagrange_state - ( + args_i.lagrange_state + + 0.5 * hi * (args_i.lagrange_cost + args_ip1.lagrange_cost) + ) end - + index = index + docp.dim_NLP_x return index end @@ -407,7 +461,7 @@ function setPathConstraints!(docp::DOCP, c, index::Int, args::ArgsAtTimeStep, v) # pure state constraints if docp.dim_x_cons > 0 - c[index:index+docp.dim_x_cons-1] = docp.state_constraints[2](ti, xi ,v) + c[index:index+docp.dim_x_cons-1] = docp.state_constraints[2](ti, xi, v) index = index + docp.dim_x_cons end @@ -416,7 +470,7 @@ function setPathConstraints!(docp::DOCP, c, index::Int, args::ArgsAtTimeStep, v) c[index:index+docp.dim_mixed_cons-1] = docp.mixed_constraints[2](ti, xi, ui, v) index = index + docp.dim_mixed_cons end - + return index end @@ -450,7 +504,7 @@ function setPathBounds!(docp::DOCP, index::Int, lb, ub) ub[index:index+docp.dim_mixed_cons-1] = docp.mixed_constraints[3] index = index + docp.dim_mixed_cons end - + return index end @@ -461,7 +515,14 @@ $(TYPEDSIGNATURES) Set the boundary and variable constraints """ -function setPointConstraints!(docp::DOCP, c, index::Int, args_0::ArgsAtTimeStep, args_f::ArgsAtTimeStep, v) +function setPointConstraints!( + docp::DOCP, + c, + index::Int, + args_0::ArgsAtTimeStep, + args_f::ArgsAtTimeStep, + v, +) ocp = docp.ocp @@ -516,8 +577,8 @@ function setPointBounds!(docp::DOCP, index::Int, lb, ub) # null initial condition for lagrangian cost state if docp.has_lagrange - lb[index] = 0. - ub[index] = 0. + lb[index] = 0.0 + ub[index] = 0.0 index = index + 1 end @@ -537,7 +598,7 @@ function DOCP_constraints_check!(cb, constraints, docp) # check constraints vs bounds # by construction only one of the two can be active - for i in 1:docp.dim_NLP_constraints + for i = 1:docp.dim_NLP_constraints if constraints[i] < docp.con_l[i] cb[i] = constraints[i] - docp.con_l[i] end @@ -557,7 +618,7 @@ Check the variables box constraints violation for the DOCP problem. function DOCP_variables_check!(vb, variables, docp) # check variables vs bounds # by construction only one of the two can be active - for i in 1:docp.dim_NLP_variables + for i = 1:docp.dim_NLP_variables if variables[i] < docp.var_l[i] vb[i] = variables[i] - docp.var_l[i] # < 0 end @@ -566,4 +627,4 @@ function DOCP_variables_check!(vb, variables, docp) end end return nothing -end \ No newline at end of file +end diff --git a/src/solution.jl b/src/solution.jl index c8f662c3..f5f6559e 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -9,13 +9,23 @@ function CTBase.OptimalControlSolution(docp, docp_solution) # retrieve data (could pass some status info too (get_status ?)) if docp.has_maximization - objective = - docp_solution.objective - else + objective = -docp_solution.objective + else objective = docp_solution.objective end # call lower level constructor - return OptimalControlSolution(docp, primal=docp_solution.solution, dual=docp_solution.multipliers, objective=objective, iterations=docp_solution.iter, constraints_violation=docp_solution.primal_feas, message=String(docp_solution.solver_specific[:internal_msg]), mult_LB=docp_solution.multipliers_L, mult_UB=docp_solution.multipliers_U) + return OptimalControlSolution( + docp, + primal = docp_solution.solution, + dual = docp_solution.multipliers, + objective = objective, + iterations = docp_solution.iter, + constraints_violation = docp_solution.primal_feas, + message = String(docp_solution.solver_specific[:internal_msg]), + mult_LB = docp_solution.multipliers_L, + mult_UB = docp_solution.multipliers_U, + ) end @@ -25,26 +35,37 @@ $(TYPEDSIGNATURES) Build OCP functional solution from the DOCP discrete solution, given as a vector. Costate will be retrieved from dual variables (multipliers) if available. """ -function CTBase.OptimalControlSolution(docp; primal=Vector(), dual=nothing, objective=nothing, iterations=nothing, constraints_violation=nothing, message=nothing, mult_LB=nothing, mult_UB=nothing) +function CTBase.OptimalControlSolution( + docp; + primal = Vector(), + dual = nothing, + objective = nothing, + iterations = nothing, + constraints_violation = nothing, + message = nothing, + mult_LB = nothing, + mult_UB = nothing, +) # time grid N = docp.dim_NLP_steps - T = zeros(N+1) - for i=1:N+1 + T = zeros(N + 1) + for i = 1:N+1 T[i] = get_unnormalized_time(primal, docp, docp.NLP_normalized_time_grid[i]) end # recover primal variables - X, U, v, box_multipliers = parse_DOCP_solution_primal(docp, primal, mult_LB=mult_LB, mult_UB=mult_UB) + X, U, v, box_multipliers = + parse_DOCP_solution_primal(docp, primal, mult_LB = mult_LB, mult_UB = mult_UB) # recompute / check objective objective_r = DOCP_objective(primal, docp) if docp.has_maximization - objective_r = - objective_r + objective_r = -objective_r end if isnothing(objective) objective = objective_r - elseif abs((objective-objective_r)/objective) > 1e-2 + elseif abs((objective - objective_r) / objective) > 1e-2 println("WARNING: recomputed objective mismatch ", objective, objective_r) end @@ -53,13 +74,26 @@ function CTBase.OptimalControlSolution(docp; primal=Vector(), dual=nothing, obje DOCP_constraints!(constraints, primal, docp) # recover costate and constraints with their multipliers - P, constraints_types, constraints_mult = parse_DOCP_solution_dual(docp, dual, constraints) + P, constraints_types, constraints_mult = + parse_DOCP_solution_dual(docp, dual, constraints) # call lowest level constructor - return OptimalControlSolution(docp, T, X, U, v, P, objective=objective, iterations=iterations, constraints_violation=constraints_violation, message=message, - constraints_types=constraints_types, constraints_mult=constraints_mult, - box_multipliers=box_multipliers) - + return OptimalControlSolution( + docp, + T, + X, + U, + v, + P, + objective = objective, + iterations = iterations, + constraints_violation = constraints_violation, + message = message, + constraints_types = constraints_types, + constraints_mult = constraints_mult, + box_multipliers = box_multipliers, + ) + end @@ -68,12 +102,12 @@ $(TYPEDSIGNATURES) Recover OCP primal variables from DOCP solution """ -function parse_DOCP_solution_primal(docp, solution; mult_LB=nothing, mult_UB=nothing) +function parse_DOCP_solution_primal(docp, solution; mult_LB = nothing, mult_UB = nothing) # state and control variables N = docp.dim_NLP_steps - X = zeros(N+1,docp.dim_NLP_x) - U = zeros(N+1,docp.dim_NLP_u) + X = zeros(N + 1, docp.dim_NLP_x) + U = zeros(N + 1, docp.dim_NLP_u) # multipliers for box constraints if isnothing(mult_LB) || length(mult_LB) == 0 @@ -82,12 +116,12 @@ function parse_DOCP_solution_primal(docp, solution; mult_LB=nothing, mult_UB=not if isnothing(mult_UB) || length(mult_UB) == 0 mult_UB = zeros(docp.dim_NLP_variables) end - mult_state_box_lower = zeros(N+1,docp.dim_NLP_x) - mult_state_box_upper = zeros(N+1,docp.dim_NLP_x) - mult_control_box_lower = zeros(N+1,docp.dim_NLP_u) - mult_control_box_upper = zeros(N+1,docp.dim_NLP_u) - mult_variable_box_lower = zeros(N+1,docp.dim_NLP_v) - mult_variable_box_upper = zeros(N+1,docp.dim_NLP_v) + mult_state_box_lower = zeros(N + 1, docp.dim_NLP_x) + mult_state_box_upper = zeros(N + 1, docp.dim_NLP_x) + mult_control_box_lower = zeros(N + 1, docp.dim_NLP_u) + mult_control_box_upper = zeros(N + 1, docp.dim_NLP_u) + mult_variable_box_lower = zeros(N + 1, docp.dim_NLP_v) + mult_variable_box_upper = zeros(N + 1, docp.dim_NLP_v) # retrieve optimization variables v = get_variable(solution, docp) @@ -95,17 +129,23 @@ function parse_DOCP_solution_primal(docp, solution; mult_LB=nothing, mult_UB=not mult_variable_box_upper = get_variable(mult_UB, docp) # loop over time steps - for i in 1:N+1 + for i = 1:N+1 # state and control variables at current step - X[i,:], U[i,:] = get_NLP_variables_at_time_step(solution, docp, i-1) + X[i, :], U[i, :] = get_NLP_variables_at_time_step(solution, docp, i - 1) # box multipliers - mult_state_box_lower[i,:], mult_control_box_lower[i,:] = get_NLP_variables_at_time_step(mult_LB, docp, i-1) - mult_state_box_upper[i,:], mult_control_box_upper[i,:] = get_NLP_variables_at_time_step(mult_UB, docp, i-1) + mult_state_box_lower[i, :], mult_control_box_lower[i, :] = + get_NLP_variables_at_time_step(mult_LB, docp, i - 1) + mult_state_box_upper[i, :], mult_control_box_upper[i, :] = + get_NLP_variables_at_time_step(mult_UB, docp, i - 1) end - box_multipliers=((mult_state_box_lower,mult_state_box_upper), (mult_control_box_lower,mult_control_box_upper), (mult_variable_box_lower,mult_variable_box_upper)) + box_multipliers = ( + (mult_state_box_lower, mult_state_box_upper), + (mult_control_box_lower, mult_control_box_upper), + (mult_variable_box_lower, mult_variable_box_upper), + ) return X, U, v, box_multipliers end @@ -129,31 +169,31 @@ function parse_DOCP_solution_dual(docp, multipliers, constraints) dbc = dim_boundary_constraints(ocp) dvc = dim_variable_constraints(ocp) - sol_control_constraints = zeros(N+1,dcc) - sol_state_constraints = zeros(N+1,dsc) - sol_mixed_constraints = zeros(N+1,dmc) + sol_control_constraints = zeros(N + 1, dcc) + sol_state_constraints = zeros(N + 1, dsc) + sol_mixed_constraints = zeros(N + 1, dmc) sol_boundary_constraints = zeros(dbc) sol_variable_constraints = zeros(dvc) - mul_control_constraints = zeros(N+1,dcc) - mul_state_constraints = zeros(N+1,dsc) - mul_mixed_constraints = zeros(N+1,dmc) + mul_control_constraints = zeros(N + 1, dcc) + mul_state_constraints = zeros(N + 1, dsc) + mul_mixed_constraints = zeros(N + 1, dmc) mul_boundary_constraints = zeros(dbc) mul_variable_constraints = zeros(dvc) # if called with multipliers = nothing, fill with zeros - if isnothing(multipliers) + if isnothing(multipliers) multipliers = zeros(docp.dim_NLP_constraints) end # loop over time steps i_c = 1 i_m = 1 - for i in 1:N+1 + for i = 1:N+1 # state equation multiplier for costate (except last step) - if i < N+1 - P[i,:] = multipliers[i_m : i_m+docp.dim_NLP_x-1] + if i < N + 1 + P[i, :] = multipliers[i_m:i_m+docp.dim_NLP_x-1] i_c += docp.dim_NLP_x # skip dynamics constraints i_m += docp.dim_NLP_x end @@ -161,22 +201,22 @@ function parse_DOCP_solution_dual(docp, multipliers, constraints) # path constraints and multipliers # pure control constraints if dcc > 0 - sol_control_constraints[i,:] = constraints[i_c : i_c+dcc-1] - mul_control_constraints[i,:] = multipliers[i_m : i_m+dcc-1] + sol_control_constraints[i, :] = constraints[i_c:i_c+dcc-1] + mul_control_constraints[i, :] = multipliers[i_m:i_m+dcc-1] i_c += dcc i_m += dcc end # pure state constraints if dsc > 0 - sol_state_constraints[i,:] = constraints[i_c : i_c+dsc-1] - mul_state_constraints[i,:] = multipliers[i_m : i_m+dsc-1] + sol_state_constraints[i, :] = constraints[i_c:i_c+dsc-1] + mul_state_constraints[i, :] = multipliers[i_m:i_m+dsc-1] i_c += dsc i_m += dsc end # mixed constraints if dmc > 0 - sol_mixed_constraints[i,:] = constraints[i_c : i_c+dmc-1] - mul_mixed_constraints[i,:] = multipliers[i_m : i_m+dmc-1] + sol_mixed_constraints[i, :] = constraints[i_c:i_c+dmc-1] + mul_mixed_constraints[i, :] = multipliers[i_m:i_m+dmc-1] i_c += dmc i_m += dmc end @@ -185,21 +225,33 @@ function parse_DOCP_solution_dual(docp, multipliers, constraints) # pointwise constraints: boundary then variables if dbc > 0 - sol_boundary_constraints[:] = constraints[i_c : i_c+dbc-1] - mul_boundary_constraints[:] = multipliers[i_m : i_m+dbc-1] + sol_boundary_constraints[:] = constraints[i_c:i_c+dbc-1] + mul_boundary_constraints[:] = multipliers[i_m:i_m+dbc-1] i_c += dbc i_m += dbc end if dvc > 0 - sol_variable_constraints[:] = constraints[i_c : i_c+dvc-1] - mul_variable_constraints[:] = multipliers[i_m : i_m+dvc-1] + sol_variable_constraints[:] = constraints[i_c:i_c+dvc-1] + mul_variable_constraints[:] = multipliers[i_m:i_m+dvc-1] i_c += dvc i_m += dvc - end + end # return tuples - constraints_types = (sol_control_constraints, sol_state_constraints, sol_mixed_constraints, sol_boundary_constraints, sol_variable_constraints) - constraints_mult = (mul_control_constraints, mul_state_constraints, mul_mixed_constraints, mul_boundary_constraints, mul_variable_constraints) + constraints_types = ( + sol_control_constraints, + sol_state_constraints, + sol_mixed_constraints, + sol_boundary_constraints, + sol_variable_constraints, + ) + constraints_mult = ( + mul_control_constraints, + mul_state_constraints, + mul_mixed_constraints, + mul_boundary_constraints, + mul_variable_constraints, + ) return P, constraints_types, constraints_mult end @@ -210,11 +262,23 @@ $(TYPEDSIGNATURES) Build OCP functional solution from DOCP vector solution (given as raw variables and multipliers plus some optional infos) """ -function CTBase.OptimalControlSolution(docp, T, X, U, v, P; - objective=0, iterations=0, constraints_violation=0, message="No msg", stopping=nothing, success=nothing, - constraints_types=(nothing,nothing,nothing,nothing,nothing), - constraints_mult=(nothing,nothing,nothing,nothing,nothing), - box_multipliers=((nothing,nothing),(nothing,nothing),(nothing,nothing))) +function CTBase.OptimalControlSolution( + docp, + T, + X, + U, + v, + P; + objective = 0, + iterations = 0, + constraints_violation = 0, + message = "No msg", + stopping = nothing, + success = nothing, + constraints_types = (nothing, nothing, nothing, nothing, nothing), + constraints_mult = (nothing, nothing, nothing, nothing, nothing), + box_multipliers = ((nothing, nothing), (nothing, nothing), (nothing, nothing)), +) ocp = docp.ocp dim_x = state_dimension(ocp) @@ -223,22 +287,24 @@ function CTBase.OptimalControlSolution(docp, T, X, U, v, P; # check that time grid is strictly increasing # if not proceed with list of indexes as time grid - if !issorted(T,lt=<=) - println("WARNING: time grid at solution is not strictly increasing, replacing with list of indices...") + if !issorted(T, lt = <=) + println( + "WARNING: time grid at solution is not strictly increasing, replacing with list of indices...", + ) println(T) - T = LinRange(0,docp.dim_NLP_steps,docp.dim_NLP_steps+1) + T = LinRange(0, docp.dim_NLP_steps, docp.dim_NLP_steps + 1) end # variables: remove additional state for lagrange cost - x = ctinterpolate(T, matrix2vec(X[:,1:dim_x], 1)) - p = ctinterpolate(T[1:end-1], matrix2vec(P[:,1:dim_x], 1)) + x = ctinterpolate(T, matrix2vec(X[:, 1:dim_x], 1)) + p = ctinterpolate(T[1:end-1], matrix2vec(P[:, 1:dim_x], 1)) u = ctinterpolate(T, matrix2vec(U, 1)) - + # force scalar output when dimension is 1 - fx = (dim_x==1) ? deepcopy(t->x(t)[1]) : deepcopy(t->x(t)) - fu = (dim_u==1) ? deepcopy(t->u(t)[1]) : deepcopy(t->u(t)) - fp = (dim_x==1) ? deepcopy(t->p(t)[1]) : deepcopy(t->p(t)) - var = (dim_v==1) ? v[1] : v + fx = (dim_x == 1) ? deepcopy(t -> x(t)[1]) : deepcopy(t -> x(t)) + fu = (dim_u == 1) ? deepcopy(t -> u(t)[1]) : deepcopy(t -> u(t)) + fp = (dim_x == 1) ? deepcopy(t -> p(t)[1]) : deepcopy(t -> p(t)) + var = (dim_v == 1) ? v[1] : v # misc infos infos = Dict{Symbol,Any}() @@ -246,29 +312,86 @@ function CTBase.OptimalControlSolution(docp, T, X, U, v, P; # +++ put interpolations here directly and reuse vectors ? # nonlinear constraints and multipliers - (control_constraints, state_constraints, mixed_constraints, boundary_constraints, variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_boundary_constraints, mult_variable_constraints) = set_constraints_and_multipliers(T, constraints_types, constraints_mult) + ( + control_constraints, + state_constraints, + mixed_constraints, + boundary_constraints, + variable_constraints, + mult_control_constraints, + mult_state_constraints, + mult_mixed_constraints, + mult_boundary_constraints, + mult_variable_constraints, + ) = set_constraints_and_multipliers(T, constraints_types, constraints_mult) # box constraints multipliers - (mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper) = set_box_multipliers(T, box_multipliers, dim_x, dim_u) + ( + mult_state_box_lower, + mult_state_box_upper, + mult_control_box_lower, + mult_control_box_upper, + mult_variable_box_lower, + mult_variable_box_upper, + ) = set_box_multipliers(T, box_multipliers, dim_x, dim_u) # build and return solution if docp.has_variable - return OptimalControlSolution(ocp; - state=fx, control=fu, objective=objective, costate=fp, time_grid=T, variable=var, iterations=iterations, stopping=stopping, message=message, success=success, infos=infos, control_constraints=control_constraints, state_constraints=state_constraints, mixed_constraints=mixed_constraints, boundary_constraints=boundary_constraints, variable_constraints=variable_constraints, - mult_control_constraints=mult_control_constraints, mult_state_constraints=mult_state_constraints, mult_mixed_constraints=mult_mixed_constraints, mult_boundary_constraints=mult_boundary_constraints, mult_variable_constraints=mult_variable_constraints, - mult_state_box_lower=mult_state_box_lower, - mult_state_box_upper=mult_state_box_upper, - mult_control_box_lower=mult_control_box_lower, - mult_control_box_upper=mult_control_box_upper, - mult_variable_box_lower=mult_variable_box_lower, - mult_variable_box_upper=mult_variable_box_upper) + return OptimalControlSolution( + ocp; + state = fx, + control = fu, + objective = objective, + costate = fp, + time_grid = T, + variable = var, + iterations = iterations, + stopping = stopping, + message = message, + success = success, + infos = infos, + control_constraints = control_constraints, + state_constraints = state_constraints, + mixed_constraints = mixed_constraints, + boundary_constraints = boundary_constraints, + variable_constraints = variable_constraints, + mult_control_constraints = mult_control_constraints, + mult_state_constraints = mult_state_constraints, + mult_mixed_constraints = mult_mixed_constraints, + mult_boundary_constraints = mult_boundary_constraints, + mult_variable_constraints = mult_variable_constraints, + mult_state_box_lower = mult_state_box_lower, + mult_state_box_upper = mult_state_box_upper, + mult_control_box_lower = mult_control_box_lower, + mult_control_box_upper = mult_control_box_upper, + mult_variable_box_lower = mult_variable_box_lower, + mult_variable_box_upper = mult_variable_box_upper, + ) else - return OptimalControlSolution(ocp; - state=fx, control=fu, objective=objective, costate=fp, time_grid=T, iterations=iterations, stopping=stopping, message=message, success=success, infos=infos, control_constraints=control_constraints, state_constraints=state_constraints, mixed_constraints=mixed_constraints, boundary_constraints=boundary_constraints, - mult_control_constraints=mult_control_constraints, mult_state_constraints=mult_state_constraints, mult_mixed_constraints=mult_mixed_constraints, mult_boundary_constraints=mult_boundary_constraints, - mult_state_box_lower=mult_state_box_lower, - mult_state_box_upper=mult_state_box_upper, - mult_control_box_lower=mult_control_box_lower, - mult_control_box_upper=mult_control_box_upper) + return OptimalControlSolution( + ocp; + state = fx, + control = fu, + objective = objective, + costate = fp, + time_grid = T, + iterations = iterations, + stopping = stopping, + message = message, + success = success, + infos = infos, + control_constraints = control_constraints, + state_constraints = state_constraints, + mixed_constraints = mixed_constraints, + boundary_constraints = boundary_constraints, + mult_control_constraints = mult_control_constraints, + mult_state_constraints = mult_state_constraints, + mult_mixed_constraints = mult_mixed_constraints, + mult_boundary_constraints = mult_boundary_constraints, + mult_state_box_lower = mult_state_box_lower, + mult_state_box_upper = mult_state_box_upper, + mult_control_box_lower = mult_control_box_lower, + mult_control_box_upper = mult_control_box_upper, + ) end end @@ -297,7 +420,18 @@ function set_constraints_and_multipliers(T, constraints_types, constraints_mult) variable_constraints = constraints_types[5] mult_variable_constraints = constraints_mult[5] - return (control_constraints, state_constraints, mixed_constraints, boundary_constraints, variable_constraints, mult_control_constraints, mult_state_constraints, mult_mixed_constraints, mult_boundary_constraints, mult_variable_constraints) + return ( + control_constraints, + state_constraints, + mixed_constraints, + boundary_constraints, + variable_constraints, + mult_control_constraints, + mult_state_constraints, + mult_mixed_constraints, + mult_boundary_constraints, + mult_variable_constraints, + ) end @@ -311,11 +445,19 @@ function set_box_multipliers(T, box_multipliers, dim_x, dim_u) # state box mult_state_box_lower, mult_state_box_upper = set_box_block(T, box_multipliers[1], dim_x) # control box - mult_control_box_lower, mult_control_box_upper = set_box_block(T, box_multipliers[2], dim_u) + mult_control_box_lower, mult_control_box_upper = + set_box_block(T, box_multipliers[2], dim_u) # variable box mult_variable_box_lower, mult_variable_box_upper = box_multipliers[3] - return (mult_state_box_lower, mult_state_box_upper, mult_control_box_lower, mult_control_box_upper, mult_variable_box_lower, mult_variable_box_upper) + return ( + mult_state_box_lower, + mult_state_box_upper, + mult_control_box_lower, + mult_control_box_upper, + mult_variable_box_lower, + mult_variable_box_upper, + ) end """ @@ -327,10 +469,10 @@ Process data related to a box type for solution building function set_box_block(T, mults, dim) mult_l, mult_u = mults if !isnothing(mult_l) && !isnothing(mult_u) - m_l = ctinterpolate(T, matrix2vec(mult_l[:,1:dim], 1)) - m_u = ctinterpolate(T, matrix2vec(mult_u[:,1:dim], 1)) + m_l = ctinterpolate(T, matrix2vec(mult_l[:, 1:dim], 1)) + m_u = ctinterpolate(T, matrix2vec(mult_u[:, 1:dim], 1)) end - return t -> m_l(t), t -> m_u(t) + return t -> m_l(t), t -> m_u(t) end diff --git a/src/solve.jl b/src/solve.jl index 76a37503..74e01d5a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -19,31 +19,45 @@ $(TYPEDSIGNATURES) Discretize an optimal control problem into a nonlinear optimization problem (ie direct transcription) """ -function direct_transcription(ocp::OptimalControlModel, +function direct_transcription( + ocp::OptimalControlModel, description...; - init=CTBase.__ocp_init(), - grid_size=__grid_size(), - time_grid=__time_grid()) + init = CTBase.__ocp_init(), + grid_size = __grid_size(), + time_grid = __time_grid(), +) # build DOCP docp = DOCP(ocp, grid_size, time_grid) - + # set bounds in DOCP variables_bounds!(docp) constraints_bounds!(docp) # set initial guess - x0 = DOCP_initial_guess(docp, OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) + x0 = DOCP_initial_guess( + docp, + OptimalControlInit( + init, + state_dim = ocp.state_dimension, + control_dim = ocp.control_dimension, + variable_dim = ocp.variable_dimension, + ), + ) # call NLP problem constructor - nlp = ADNLPModel!(x -> DOCP_objective(x, docp), + nlp = ADNLPModel!( + x -> DOCP_objective(x, docp), x0, - docp.var_l, docp.var_u, - (c, x) -> DOCP_constraints!(c, x, docp), - docp.con_l, docp.con_u, - backend = :optimized) + docp.var_l, + docp.var_u, + (c, x) -> DOCP_constraints!(c, x, docp), + docp.con_l, + docp.con_u, + backend = :optimized, + ) -return docp, nlp + return docp, nlp end @@ -56,7 +70,15 @@ Set initial guess in the DOCP function set_initial_guess(docp::DOCP, nlp, init) ocp = docp.ocp - nlp.meta.x0 .= DOCP_initial_guess(docp, OptimalControlInit(init, state_dim=ocp.state_dimension, control_dim=ocp.control_dimension, variable_dim=ocp.variable_dimension)) + nlp.meta.x0 .= DOCP_initial_guess( + docp, + OptimalControlInit( + init, + state_dim = ocp.state_dimension, + control_dim = ocp.control_dimension, + variable_dim = ocp.variable_dimension, + ), + ) end @@ -65,17 +87,26 @@ $(TYPEDSIGNATURES) Solve an OCP with a direct method """ -function direct_solve(ocp::OptimalControlModel, description::Symbol...; - init=CTBase.__ocp_init(), - grid_size::Int=CTDirect.__grid_size(), - time_grid=CTDirect.__time_grid(), - kwargs...) +function direct_solve( + ocp::OptimalControlModel, + description::Symbol...; + init = CTBase.__ocp_init(), + grid_size::Int = CTDirect.__grid_size(), + time_grid = CTDirect.__time_grid(), + kwargs..., +) method = getFullDescription(description, available_methods()) #println(method) # build discretized OCP, including initial guess - docp, nlp = direct_transcription(ocp, description, init=init, grid_size=grid_size, time_grid=time_grid) + docp, nlp = direct_transcription( + ocp, + description, + init = init, + grid_size = grid_size, + time_grid = time_grid, + ) # solve DOCP if :ipopt ∈ method @@ -90,7 +121,7 @@ function direct_solve(ocp::OptimalControlModel, description::Symbol...; # build and return OCP solution return OptimalControlSolution(docp, docp_solution) - end +end # placeholders (see CTSolveExt*** extensions) diff --git a/src/utils.jl b/src/utils.jl index fba016b8..7c870b3e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -9,14 +9,14 @@ function get_variable(xu, docp) nx = docp.dim_NLP_x m = docp.dim_NLP_u N = docp.dim_NLP_steps - offset = (nx+m) * (N+1) + offset = (nx + m) * (N + 1) if docp.dim_NLP_v == 1 - return xu[offset + 1] + return xu[offset+1] #return xu[end] else #return xu[end-docp.dim_NLP_v+1:end] - return xu[offset + 1: offset + docp.dim_NLP_v] + return xu[offset+1:offset+docp.dim_NLP_v] end else return Float64[] @@ -35,8 +35,8 @@ function get_single_variable(xu, docp, i::Int) nx = docp.dim_NLP_x m = docp.dim_NLP_u N = docp.dim_NLP_steps - offset = (nx+m) * (N+1) - return xu[offset + i] + offset = (nx + m) * (N + 1) + return xu[offset+i] else error("Tring to access variable in variable independent problem") end @@ -55,7 +55,7 @@ function get_variables_at_time_step(xu, docp, i) n = docp.dim_OCP_x m = docp.dim_NLP_u offset = (nx + m) * i - + # retrieve scalar/vector OCP state (w/o lagrange state) if n == 1 xi = xu[offset+1] @@ -67,7 +67,7 @@ function get_variables_at_time_step(xu, docp, i) if docp.has_lagrange xli = xu[offset+nx] else - xli = 0. + xli = 0.0 end # retrieve scalar/vector control @@ -86,7 +86,7 @@ function get_NLP_variables_at_time_step(xu, docp, i) nx = docp.dim_NLP_x m = docp.dim_NLP_u offset = (nx + m) * i - + xi = xu[offset+1:offset+nx] ui = xu[offset+nx+1:offset+nx+m] @@ -181,7 +181,7 @@ Get actual (un-normalized) time value """ function get_unnormalized_time(xu, docp, t_normalized) t0 = get_initial_time(xu, docp) - tf = get_final_time(xu, docp) + tf = get_final_time(xu, docp) return t0 + t_normalized * (tf - t0) end @@ -233,10 +233,10 @@ function set_variables_at_time_step!(xu, x_init, u_init, docp, i) # NB. only set first the actual state variables from the OCP (not the possible additional state for lagrange cost) if !isnothing(x_init) - xu[offset + 1 : offset + n] .= x_init + xu[offset+1:offset+n] .= x_init end if !isnothing(u_init) - xu[offset + nx + 1 : offset + nx + m] .= u_init + xu[offset+nx+1:offset+nx+m] .= u_init end end @@ -246,7 +246,7 @@ $(TYPEDSIGNATURES) Set optimization variables in the NLP variables (for initial guess) """ function set_variable!(xu, v_init, docp) - xu[end-docp.dim_NLP_v+1 : end] .= v_init + xu[end-docp.dim_NLP_v+1:end] .= v_init end @@ -255,8 +255,7 @@ $(TYPEDSIGNATURES) Build initial guess for discretized problem """ -function DOCP_initial_guess(docp, - init::OptimalControlInit=OptimalControlInit()) +function DOCP_initial_guess(docp, init::OptimalControlInit = OptimalControlInit()) # default initialization # note: internal variables (lagrange cost, k_i for RK schemes) will keep these default values @@ -269,9 +268,15 @@ function DOCP_initial_guess(docp, end # set state / control variables if provided - for i in 0:docp.dim_NLP_steps + for i = 0:docp.dim_NLP_steps ti = get_time_at_time_step(NLP_X, docp, i) - set_variables_at_time_step!(NLP_X, init.state_init(ti), init.control_init(ti), docp, i) + set_variables_at_time_step!( + NLP_X, + init.state_init(ti), + init.control_init(ti), + docp, + i, + ) end return NLP_X @@ -280,4 +285,3 @@ end # placeholders (see CTDirectExt) function export_ocp_solution end function import_ocp_solution end - diff --git a/test/manual_test.jl b/test/manual_test.jl index 3f66eec2..900bd335 100644 --- a/test/manual_test.jl +++ b/test/manual_test.jl @@ -3,4 +3,4 @@ include("deps.jl") using Plots using Printf -using Test \ No newline at end of file +using Test diff --git a/test/old/bocop_goddard.jl b/test/old/bocop_goddard.jl index 9e1cb825..4ba4d41b 100644 --- a/test/old/bocop_goddard.jl +++ b/test/old/bocop_goddard.jl @@ -16,25 +16,25 @@ v0 = 0 vmax = 0.1 m0 = 1 mf = 0.6 -x0 = [ r0, v0, m0 ] +x0 = [r0, v0, m0] ocp = Model() time!(ocp, :initial, t0) state!(ocp, 3, ["r", "v", "m"]) control!(ocp, 1) constraint!(ocp, :initial, x0, :initial_constraint) -constraint!(ocp, :state, 1:3, [1,0,0.6], [Inf,0.1,1], :box_state) +constraint!(ocp, :state, 1:3, [1, 0, 0.6], [Inf, 0.1, 1], :box_state) constraint!(ocp, :control, Index(1), 0, 1, :box_control) -objective!(ocp, :mayer, (t0, x0, tf, xf) -> xf[1], :max) +objective!(ocp, :mayer, (t0, x0, tf, xf) -> xf[1], :max) function F0(x) r, v, m = x - D = Cd * v^2 * exp(-β*(r - 1)) - return [ v, -D/m - 1/r^2, 0 ] + D = Cd * v^2 * exp(-β * (r - 1)) + return [v, -D / m - 1 / r^2, 0] end function F1(x) r, v, m = x - return [ 0, Tmax/m, -b*Tmax ] + return [0, Tmax / m, -b * Tmax] end -f(x, u) = F0(x) + u*F1(x) +f(x, u) = F0(x) + u * F1(x) constraint!(ocp, :dynamics, f) # settings @@ -42,10 +42,31 @@ init = [1.01, 0.05, 0.8, 0.1] N = 10 # dummy run then 2 runs -@time sol = solve(ocp, grid_size=N, print_level=0, tol=1e-12, mu_strategy="adaptive", init=init) -println("N=",N, " Obj ",sol.objective," Iter ", sol.iterations) -@time sol = solve(ocp, grid_size=N, print_level=0, tol=1e-12, mu_strategy="adaptive", init=init) -@time sol = solve(ocp, grid_size=N, print_level=0, tol=1e-12, mu_strategy="adaptive", init=init) +@time sol = solve( + ocp, + grid_size = N, + print_level = 0, + tol = 1e-12, + mu_strategy = "adaptive", + init = init, +) +println("N=", N, " Obj ", sol.objective, " Iter ", sol.iterations) +@time sol = solve( + ocp, + grid_size = N, + print_level = 0, + tol = 1e-12, + mu_strategy = "adaptive", + init = init, +) +@time sol = solve( + ocp, + grid_size = N, + print_level = 0, + tol = 1e-12, + mu_strategy = "adaptive", + init = init, +) # benchmark. Not very practical -_- #BenchmarkTools.DEFAULT_PARAMETERS.samples = 2 diff --git a/test/old/local_test_exponential.jl b/test/old/local_test_exponential.jl index 77f7285a..f32f22b6 100644 --- a/test/old/local_test_exponential.jl +++ b/test/old/local_test_exponential.jl @@ -79,39 +79,39 @@ sol = solve(ocp2, grid_size=100, print_level=5, tol=1e-12) =# # try min tf lagrange and mayer -ocp3 = Model(variable=true) +ocp3 = Model(variable = true) state!(ocp3, 2) control!(ocp3, 1) variable!(ocp3, 1) time!(ocp3, 0, Index(1)) #variable!(ocp3, 2) #time!(ocp3, Index(1), Index(2)) -constraint!(ocp3, :initial, [0,0], :initial_constraint) -constraint!(ocp3, :final, [1,0], :final_constraint) +constraint!(ocp3, :initial, [0, 0], :initial_constraint) +constraint!(ocp3, :final, [1, 0], :final_constraint) constraint!(ocp3, :control, -1, 1, :control_constraint) constraint!(ocp3, :variable, 0.1, 10, :variable_constraint) #constraint!(ocp3, :variable, [0.1, 0.1], [10, 10], :variable_constraint) -dynamics!(ocp3, (x, u, v) -> [x[2], u]) +dynamics!(ocp3, (x, u, v) -> [x[2], u]) #objective!(ocp3, :lagrange, (x, u, v) -> 1) # min free tf ok #objective!(ocp3, :mayer, (x0, xf, v) -> v[1]) # min free tf ok objective!(ocp3, :mayer, (x0, xf, v) -> v) # scalar variable #objective!(ocp3, :mayer, (x0, xf, v) -> - v[1]) # max t0 free to and tf ok #objective!(ocp3, :mayer, (x0, xf, v) -> v[1], :max) # max t0 free to and tf ok -sol3 = solve(ocp3, grid_size=100, print_level=5, tol=1e-12) +sol3 = solve(ocp3, grid_size = 100, print_level = 5, tol = 1e-12) plot(sol3) # min tf, abstract definition @def ocp4 begin tf ∈ R, variable - t ∈ [ 0, tf ], time + t ∈ [0, tf], time x ∈ R², state u ∈ R, control -1 ≤ u(t) ≤ 1 - x(0) == [ 0, 0 ] - x(tf) == [ 1, 0 ] - 0.1 ≤ tf ≤ Inf - ẋ(t) == [ x₂(t), u(t) ] + x(0) == [0, 0] + x(tf) == [1, 0] + 0.1 ≤ tf ≤ Inf + ẋ(t) == [x₂(t), u(t)] tf → min end -sol4 = solve(ocp4, grid_size=100, print_level=5, tol=1e-12) +sol4 = solve(ocp4, grid_size = 100, print_level = 5, tol = 1e-12) plot(sol4) diff --git a/test/old/manual_simple_integrator_energy.jl b/test/old/manual_simple_integrator_energy.jl index d33b1acc..ed1896a8 100644 --- a/test/old/manual_simple_integrator_energy.jl +++ b/test/old/manual_simple_integrator_energy.jl @@ -15,7 +15,7 @@ ocp = prob.model # solve #sol = solve(ocp, grid_size=10, print_level=5) -sol = solve(ocp, grid_size=100, print_level=5, tol=1e-12, init=nothing) +sol = solve(ocp, grid_size = 100, print_level = 5, tol = 1e-12, init = nothing) p1 = plot(sol) @@ -27,4 +27,4 @@ sol = solve(ocp, grid_size=100, print_level=5, tol=1e-12, mu_strategy="adaptive" # plot p2 = plot(sol) =# -plot(p1) \ No newline at end of file +plot(p1) diff --git a/test/old/manual_test_all_constraints.jl b/test/old/manual_test_all_constraints.jl index 1fde99df..792d8642 100644 --- a/test/old/manual_test_all_constraints.jl +++ b/test/old/manual_test_all_constraints.jl @@ -26,7 +26,14 @@ check_control_box = false # solve problem println("Solving test problem...") init = [1.01, 0.05, 0.8, 0.1] -sol = solve(ocp, grid_size=50, print_level=0, tol=1e-8, mu_strategy="adaptive", init=init) +sol = solve( + ocp, + grid_size = 50, + print_level = 0, + tol = 1e-8, + mu_strategy = "adaptive", + init = init, +) t0 = sol.times[1] tf = last(sol.times) t = sol.times @@ -38,12 +45,26 @@ end if plot_state_constraints ncx = sol.infos[:dim_state_constraints] if ncx > 0 - PCX = Array{Plots.Plot, 1}(undef, ncx); - for i in 1:ncx - PCX[i] = plot(t -> sol.infos[:state_constraints](t)[i], t0, tf, label="state_constraint v <= 0.1", legend=:topleft) - PCX[i] = plot!(twinx(), t -> sol.infos[:mult_state_constraints](t)[i], t0, tf, color=:red, label="multiplier", xticks=:none) #2nd scale + PCX = Array{Plots.Plot,1}(undef, ncx) + for i = 1:ncx + PCX[i] = plot( + t -> sol.infos[:state_constraints](t)[i], + t0, + tf, + label = "state_constraint v <= 0.1", + legend = :topleft, + ) + PCX[i] = plot!( + twinx(), + t -> sol.infos[:mult_state_constraints](t)[i], + t0, + tf, + color = :red, + label = "multiplier", + xticks = :none, + ) #2nd scale end - plot(PCX..., layout = (ncx,1)) + plot(PCX..., layout = (ncx, 1)) end end @@ -51,12 +72,26 @@ end if plot_control_constraints ncu = sol.infos[:dim_control_constraints] if ncu > 0 - PCU = Array{Plots.Plot, 1}(undef, ncu); - for i in 1:ncu - PCU[i] = plot(t -> sol.infos[:control_constraints](t)[i], t0, tf, label="control_constraint u <= 1", legend=:topleft) - PCU[i] = plot!(twinx(), t -> sol.infos[:mult_control_constraints](t)[i], t0, tf, color=:red, label="multiplier", xticks=:none) #2nd scale + PCU = Array{Plots.Plot,1}(undef, ncu) + for i = 1:ncu + PCU[i] = plot( + t -> sol.infos[:control_constraints](t)[i], + t0, + tf, + label = "control_constraint u <= 1", + legend = :topleft, + ) + PCU[i] = plot!( + twinx(), + t -> sol.infos[:mult_control_constraints](t)[i], + t0, + tf, + color = :red, + label = "multiplier", + xticks = :none, + ) #2nd scale end - plot(PCU..., layout = (ncu,1)) + plot(PCU..., layout = (ncu, 1)) end end @@ -64,31 +99,93 @@ end if plot_mixed_constraints ncxu = sol.infos[:dim_mixed_constraints] if ncxu > 0 - PCXU = Array{Plots.Plot, 1}(undef, ncxu); - for i in 1:ncxu - PCXU[i] = plot(t -> sol.infos[:mixed_constraints](t)[i], t0, tf, label="mixed_constraint m >= 0.6", legend=:topleft) - PCXU[i] = plot!(twinx(), t -> sol.infos[:mult_mixed_constraints](t)[i], t0, tf, color=:red, label="multiplier", xticks=:none) #2nd scale + PCXU = Array{Plots.Plot,1}(undef, ncxu) + for i = 1:ncxu + PCXU[i] = plot( + t -> sol.infos[:mixed_constraints](t)[i], + t0, + tf, + label = "mixed_constraint m >= 0.6", + legend = :topleft, + ) + PCXU[i] = plot!( + twinx(), + t -> sol.infos[:mult_mixed_constraints](t)[i], + t0, + tf, + color = :red, + label = "multiplier", + xticks = :none, + ) #2nd scale end - plot(PCXU..., layout = (ncxu,1)) + plot(PCXU..., layout = (ncxu, 1)) end end # state box +++ not generic ! if plot_state_box - PX = Array{Plots.Plot, 1}(undef, 2); # only boxes on r, v - for i in 1:2 - PX[i] = plot(t -> sol.state(t)[i], t0, tf, label="state with box mult (green: LB, red: UB)", legend=:topleft) - PX[i] = plot!(twinx(),t -> sol.infos[:mult_state_box_lower](t)[i], t0, tf, color=:green, xticks=:none, label=:none, linestyle=:dash) - PX[i] = plot!(twinx(),t -> sol.infos[:mult_state_box_upper](t)[i], t0, tf, color=:red, xticks=:none, label=:none, linestyle=:dash) + PX = Array{Plots.Plot,1}(undef, 2) # only boxes on r, v + for i = 1:2 + PX[i] = plot( + t -> sol.state(t)[i], + t0, + tf, + label = "state with box mult (green: LB, red: UB)", + legend = :topleft, + ) + PX[i] = plot!( + twinx(), + t -> sol.infos[:mult_state_box_lower](t)[i], + t0, + tf, + color = :green, + xticks = :none, + label = :none, + linestyle = :dash, + ) + PX[i] = plot!( + twinx(), + t -> sol.infos[:mult_state_box_upper](t)[i], + t0, + tf, + color = :red, + xticks = :none, + label = :none, + linestyle = :dash, + ) end PPX = plot(PX..., layout = (2, 1)) end # control box +++ not generic ! if plot_control_box - PU = plot(t -> sol.control(t)[1], t0, tf, label="control with box mult (green: LB u>=0, red: UB unused)", legend=:topleft) - PU = plot!(twinx(),t -> sol.infos[:mult_control_box_lower](t)[1], t0, tf, color=:green, xticks=:none, label=:none, linestyle=:dash) - PU = plot!(twinx(),t -> sol.infos[:mult_control_box_upper](t)[1], t0, tf, color=:red, xticks=:none, label=:none, linestyle=:dash) + PU = plot( + t -> sol.control(t)[1], + t0, + tf, + label = "control with box mult (green: LB u>=0, red: UB unused)", + legend = :topleft, + ) + PU = plot!( + twinx(), + t -> sol.infos[:mult_control_box_lower](t)[1], + t0, + tf, + color = :green, + xticks = :none, + label = :none, + linestyle = :dash, + ) + PU = plot!( + twinx(), + t -> sol.infos[:mult_control_box_upper](t)[1], + t0, + tf, + color = :red, + xticks = :none, + label = :none, + linestyle = :dash, + ) plot(PU) end @@ -98,11 +195,11 @@ end if check_state_constraints a = sol.infos[:mult_state_constraints](t)[1] println(a) - println(all(>=(0), a) ) + println(all(>=(0), a)) b = sol.infos[:state_constraints](t)[1] println(b) - println(a.*b) - println(norm(a.*b)) + println(a .* b) + println(norm(a .* b)) end # control constraint u <= 1 diff --git a/test/old/manual_test_exponential.jl b/test/old/manual_test_exponential.jl index 4301570a..010f0703 100644 --- a/test/old/manual_test_exponential.jl +++ b/test/old/manual_test_exponential.jl @@ -7,12 +7,12 @@ using CTProblems #remove_constraint!(ocp,:control_constraint) #constraint!(ocp, :control, 0, 1, :control_constraint2) -n=1 -m=1 -t0=0 -tf=1 -x0=-1 -xf=0 +n = 1 +m = 1 +t0 = 0 +tf = 1 +x0 = -1 +xf = 0 ocp = Model() state!(ocp, n) # dimension of the state control!(ocp, m) # dimension of the control @@ -25,5 +25,5 @@ constraint!(ocp, :dynamics, (x, u) -> -x + u) objective!(ocp, :lagrange, (x, u) -> abs(u)) #objective!(ocp, :lagrange, (x, u) -> sqrt(u*u)) -sol = solve(ocp, grid_size=20, print_level=5) +sol = solve(ocp, grid_size = 20, print_level = 5) p1 = plot(sol) diff --git a/test/old/manual_test_integrator.jl b/test/old/manual_test_integrator.jl index 4b6faf9e..392551be 100644 --- a/test/old/manual_test_integrator.jl +++ b/test/old/manual_test_integrator.jl @@ -8,20 +8,20 @@ prob = Problem(:integrator, :dim2, :energy) ocp = prob.model # initial guess (constant state and control functions) -init = [1., 0.5, 0.3] +init = [1.0, 0.5, 0.3] # solve #sol = solve(ocp, grid_size=10, print_level=5) -sol = solve(ocp, grid_size=100, print_level=5, tol=1e-12, init=init) +sol = solve(ocp, grid_size = 100, print_level = 5, tol = 1e-12, init = init) p1 = plot(sol) # control box constraint!(ocp, :control, -4.01, 4.01, :control_con1) -sol = solve(ocp, grid_size=100, print_level=5, tol=1e-12, init=init) +sol = solve(ocp, grid_size = 100, print_level = 5, tol = 1e-12, init = init) # plot -p2 = plot(sol) +p2 = plot(sol) -plot(p1,p2) \ No newline at end of file +plot(p1, p2) diff --git a/test/old/test_goddard.jl b/test/old/test_goddard.jl index d7862438..d13aa3b9 100644 --- a/test/old/test_goddard.jl +++ b/test/old/test_goddard.jl @@ -2,9 +2,9 @@ println("Goddard test") prob = Problem(:goddard, :all_constraint) init = [1.01, 0.05, 0.8, 0.1] -sol = solve(prob.model, grid_size=10, print_level=0, init=init) +sol = solve(prob.model, grid_size = 10, print_level = 0, init = init) @testset verbose = true showtiming = true ":goddard :all_constraints" begin - @test sol.objective ≈ prob.solution.objective rtol=1e-2 + @test sol.objective ≈ prob.solution.objective rtol = 1e-2 end # goddard, abstract def @@ -18,45 +18,45 @@ v0 = 0 vmax = 0.1 m0 = 1 mf = 0.6 -x0 = [ r0, v0, m0 ] +x0 = [r0, v0, m0] @def ocp begin tf, variable - t ∈ [ t0, tf ], time + t ∈ [t0, tf], time x ∈ R³, state u ∈ R, control - + r = x₁ v = x₂ m = x₃ - - x(t0) == [ r0, v0, m0 ] - 0 ≤ u(t) ≤ 1 - r0 ≤ r(t) ≤ Inf, (1) - 0 ≤ v(t) ≤ vmax, (2) - mf ≤ m(t) ≤ m0, (3) + + x(t0) == [r0, v0, m0] + 0 ≤ u(t) ≤ 1 + r0 ≤ r(t) ≤ Inf, (1) + 0 ≤ v(t) ≤ vmax, (2) + mf ≤ m(t) ≤ m0, (3) ẋ(t) == F0(x(t)) + u(t) * F1(x(t)) - + r(tf) → max - + end F0(x) = begin r, v, m = x - D = Cd * v^2 * exp(-β*(r - 1)) - F = [ v, -D/m - 1/r^2, 0 ] + D = Cd * v^2 * exp(-β * (r - 1)) + F = [v, -D / m - 1 / r^2, 0] return F end F1(x) = begin r, v, m = x - F = [ 0, Tmax/m, -b*Tmax ] + F = [0, Tmax / m, -b * Tmax] return F end -sol = solve(ocp, grid_size=10, print_level=5) +sol = solve(ocp, grid_size = 10, print_level = 5) @testset verbose = true showtiming = true ":goddard :all_constraints" begin - @test sol.objective ≈ prob.solution.objective rtol=1e-2 + @test sol.objective ≈ prob.solution.objective rtol = 1e-2 end diff --git a/test/old/test_integrator.jl b/test/old/test_integrator.jl index 2b924e0e..dcb3cf6b 100644 --- a/test/old/test_integrator.jl +++ b/test/old/test_integrator.jl @@ -1,6 +1,6 @@ # double integrator - energy min println("Double integrator test") -prob = Problem(:integrator, :dim2, :energy) +prob = Problem(:integrator, :dim2, :energy) ocp = prob.model @testset verbose = true showtiming = true ":integrator :dim2 :energy" begin @@ -14,30 +14,30 @@ ocp = prob.model @test ctd.dim_control_box == 0 @test ctd.dim_state_box == 0 @test ctd.has_control_constraints == false - @test ctd.has_state_constraints == false - @test ctd.has_mixed_constraints == false + @test ctd.has_state_constraints == false + @test ctd.has_mixed_constraints == false @test ctd.has_boundary_conditions == true - @test ctd.has_control_box == false - @test ctd.has_state_box == false - @test ctd.has_lagrange_cost == true - @test ctd.has_mayer_cost == false - @test ctd.dim_NLP_state == ocp.state_dimension + 1 - @test ctd.mayer === nothing - f_Mayer_test(t,x,u)=[ocp.dynamics(t,x,u);ocp.lagrange(t,x,u)] - @test ctd.dynamics_lagrange_to_mayer(0,[0;2],1) == f_Mayer_test(0,[0;2],1) - @test ctd.has_free_final_time == false - @test ctd.dim_NLP_variables == (N+1)*4 - @test ctd.dim_NLP_constraints == 3*N+5 + @test ctd.has_control_box == false + @test ctd.has_state_box == false + @test ctd.has_lagrange_cost == true + @test ctd.has_mayer_cost == false + @test ctd.dim_NLP_state == ocp.state_dimension + 1 + @test ctd.mayer === nothing + f_Mayer_test(t, x, u) = [ocp.dynamics(t, x, u); ocp.lagrange(t, x, u)] + @test ctd.dynamics_lagrange_to_mayer(0, [0; 2], 1) == f_Mayer_test(0, [0; 2], 1) + @test ctd.has_free_final_time == false + @test ctd.dim_NLP_variables == (N + 1) * 4 + @test ctd.dim_NLP_constraints == 3 * N + 5 end @testset verbose = true showtiming = true "constraints" begin lb, ub = CTDirect.constraints_bounds(ctd) l_var, u_var = CTDirect.variables_bounds(ctd) - true_lb = zeros(3*N) # test without boundary conditions because of the dictionary + true_lb = zeros(3 * N) # test without boundary conditions because of the dictionary # true_lb[end-4:end] = [-1,0,0,0,0] @test lb[1:3*N] == true_lb @test ub[1:3*N] == true_lb - @test l_var == -Inf*ones((N+1)*4) + @test l_var == -Inf * ones((N + 1) * 4) @test u_var == -l_var end @@ -45,18 +45,18 @@ ocp = prob.model # Tests on the numerical solution @testset verbose = true showtiming = true "numerical solution" begin - sol = solve(ocp, grid_size=100, print_level=0) + sol = solve(ocp, grid_size = 100, print_level = 0) u_sol(t) = prob.solution.control(t)[1] u = t -> sol.control(t)[1] x_sol(t) = prob.solution.state(t) x = t -> sol.state(t)[1:2] T = sol.times - dT = T[2:end]-T[1:end-1] + dT = T[2:end] - T[1:end-1] N = length(T) - @test distance_infty(x,x_sol,T) ≈ 0 atol = 0.01 - @test distance_L2(u,u_sol,T) ≈ 0 atol=1e-1 - @test sol.objective ≈ prob.solution.objective atol=1e-2 + @test distance_infty(x, x_sol, T) ≈ 0 atol = 0.01 + @test distance_L2(u, u_sol, T) ≈ 0 atol = 1e-1 + @test sol.objective ≈ prob.solution.objective atol = 1e-2 end diff --git a/test/old/test_integrator_constraints.jl b/test/old/test_integrator_constraints.jl index c996272c..bcef1838 100644 --- a/test/old/test_integrator_constraints.jl +++ b/test/old/test_integrator_constraints.jl @@ -1,115 +1,116 @@ # double integrator - energy min println("Double integrator test") -prob = Problem(:integrator, :dim2, :energy) +prob = Problem(:integrator, :dim2, :energy) ocp = prob.model # solve println("Is solvable ? ", CTDirect.is_solvable(ocp)) @testset verbose = true showtiming = true ":integrator :dim2 :energy with constraints" begin - + @testset verbose = true showtiming = true "control_box_constraints" begin - umax = 5. - N = 2 - constraint!(ocp, :control, -umax, umax, :control_con1) - ctd = CTDirect.CTDirect_data(ocp, N, nothing) - @test ctd.dim_control_constraints == 0 - @test ctd.has_control_box == true - @test ctd.has_control_constraints == false - @test ctd.dim_control_box == 1 - @test ctd.dim_NLP_state == 3 + umax = 5.0 + N = 2 + constraint!(ocp, :control, -umax, umax, :control_con1) + ctd = CTDirect.CTDirect_data(ocp, N, nothing) + @test ctd.dim_control_constraints == 0 + @test ctd.has_control_box == true + @test ctd.has_control_constraints == false + @test ctd.dim_control_box == 1 + @test ctd.dim_NLP_state == 3 - lb, ub = CTDirect.constraints_bounds(ctd) - l_var, u_var = CTDirect.variables_bounds(ctd) - true_lb = zeros(3*N) # test without boundary conditions because of the dictionary - # true_lb[end-4:end] = [-1,0,0,0,0] - true_l_var = -Inf*ones((N+1)*4) - true_l_var[3*(N+1)+1:4*(N+1)] .= -umax - @test lb[1:end-5] == true_lb - @test ub[1:end-5] == -true_lb - @test l_var == true_l_var - @test u_var == -l_var - end + lb, ub = CTDirect.constraints_bounds(ctd) + l_var, u_var = CTDirect.variables_bounds(ctd) + true_lb = zeros(3 * N) # test without boundary conditions because of the dictionary + # true_lb[end-4:end] = [-1,0,0,0,0] + true_l_var = -Inf * ones((N + 1) * 4) + true_l_var[3*(N+1)+1:4*(N+1)] .= -umax + @test lb[1:end-5] == true_lb + @test ub[1:end-5] == -true_lb + @test l_var == true_l_var + @test u_var == -l_var + end - @testset verbose = true showtiming = true "control_constraints" begin - umax = 5 - remove_constraint!(ocp, :control_con1) - constraint!(ocp, :control, u -> u, -umax, umax, :control_con2) - N = 3 - ctd = CTDirect.CTDirect_data(ocp, N, nothing) - @test ctd.dim_control_constraints == 1 - @test ctd.has_control_box == false - @test ctd.has_control_constraints == true - @test ctd.dim_control_box == 0 - @test ctd.dim_control_constraints == 1 + @testset verbose = true showtiming = true "control_constraints" begin + umax = 5 + remove_constraint!(ocp, :control_con1) + constraint!(ocp, :control, u -> u, -umax, umax, :control_con2) + N = 3 + ctd = CTDirect.CTDirect_data(ocp, N, nothing) + @test ctd.dim_control_constraints == 1 + @test ctd.has_control_box == false + @test ctd.has_control_constraints == true + @test ctd.dim_control_box == 0 + @test ctd.dim_control_constraints == 1 - lb, ub = CTDirect.constraints_bounds(ctd) - l_var, u_var = CTDirect.variables_bounds(ctd) - true_lb = zeros(4*N+1) # test without boundary conditions because of the dictionary - true_lb[4:4:end] .= -umax - true_lb[end] = -umax - # true_lb[end-4:end] = [-1,0,0,0,0] - true_l_var = -Inf*ones((N+1)*4) - @test lb[1:end-5] == true_lb - @test ub[1:end-5] == -true_lb - @test l_var == true_l_var - @test u_var == -l_var - end + lb, ub = CTDirect.constraints_bounds(ctd) + l_var, u_var = CTDirect.variables_bounds(ctd) + true_lb = zeros(4 * N + 1) # test without boundary conditions because of the dictionary + true_lb[4:4:end] .= -umax + true_lb[end] = -umax + # true_lb[end-4:end] = [-1,0,0,0,0] + true_l_var = -Inf * ones((N + 1) * 4) + @test lb[1:end-5] == true_lb + @test ub[1:end-5] == -true_lb + @test l_var == true_l_var + @test u_var == -l_var + end - @testset verbose = true showtiming = true "mixed_constraints" begin - umax = 5 - remove_constraint!(ocp, :control_con2) - constraint!(ocp, :mixed, (x,u) -> u, -umax, umax, :control_con3) - N = 3 - ctd = CTDirect.CTDirect_data(ocp, N, nothing) - @test ctd.dim_control_constraints == 0 - @test ctd.has_control_box == false - @test ctd.has_control_constraints == false - @test ctd.has_mixed_constraints == true - @test ctd.dim_control_constraints == 0 - @test ctd.dim_mixed_constraints == 1 + @testset verbose = true showtiming = true "mixed_constraints" begin + umax = 5 + remove_constraint!(ocp, :control_con2) + constraint!(ocp, :mixed, (x, u) -> u, -umax, umax, :control_con3) + N = 3 + ctd = CTDirect.CTDirect_data(ocp, N, nothing) + @test ctd.dim_control_constraints == 0 + @test ctd.has_control_box == false + @test ctd.has_control_constraints == false + @test ctd.has_mixed_constraints == true + @test ctd.dim_control_constraints == 0 + @test ctd.dim_mixed_constraints == 1 - lb, ub = CTDirect.constraints_bounds(ctd) - l_var, u_var = CTDirect.variables_bounds(ctd) - true_lb = zeros(4*N+1) # test without boundary conditions because of the dictionary - true_lb[4:4:end] .= -umax - true_lb[end] = -umax - # true_lb[end-4:end] = [-1,0,0,0,0] - true_l_var = -Inf*ones((N+1)*4) - @test lb[1:end-5] == true_lb - @test ub[1:end-5] == -true_lb - @test l_var == true_l_var - @test u_var == -l_var - end + lb, ub = CTDirect.constraints_bounds(ctd) + l_var, u_var = CTDirect.variables_bounds(ctd) + true_lb = zeros(4 * N + 1) # test without boundary conditions because of the dictionary + true_lb[4:4:end] .= -umax + true_lb[end] = -umax + # true_lb[end-4:end] = [-1,0,0,0,0] + true_l_var = -Inf * ones((N + 1) * 4) + @test lb[1:end-5] == true_lb + @test ub[1:end-5] == -true_lb + @test l_var == true_l_var + @test u_var == -l_var + end - @testset verbose = true showtiming = true "state_box_constraints" begin - x_min1 = 1; x_min2 = -1.; x_max1 = 2; x_max2 = 5. - N = 2 - remove_constraint!(ocp, :control_con3) - constraint!(ocp, :state, [x_min1, x_min2], [x_max1, x_max2], :state_con1) - ctd = CTDirect.CTDirect_data(ocp, N, nothing) - @test ctd.dim_control_constraints == 0 - @test ctd.has_control_box == false - @test ctd.has_control_constraints == false - @test ctd.dim_control_box == 0 - @test ctd.has_state_box == true + @testset verbose = true showtiming = true "state_box_constraints" begin + x_min1 = 1 + x_min2 = -1.0 + x_max1 = 2 + x_max2 = 5.0 + N = 2 + remove_constraint!(ocp, :control_con3) + constraint!(ocp, :state, [x_min1, x_min2], [x_max1, x_max2], :state_con1) + ctd = CTDirect.CTDirect_data(ocp, N, nothing) + @test ctd.dim_control_constraints == 0 + @test ctd.has_control_box == false + @test ctd.has_control_constraints == false + @test ctd.dim_control_box == 0 + @test ctd.has_state_box == true - lb, ub = CTDirect.constraints_bounds(ctd) - l_var, u_var = CTDirect.variables_bounds(ctd) - true_lb = zeros(3*N) # test without boundary conditions because of the dictionary - # true_lb[end-4:end] = [-1,0,0,0,0] - true_l_var = -Inf*ones((N+1)*4) - true_l_var[1:3:end-(N+1)] .= x_min1 - true_l_var[2:3:end-(N+1)] .= x_min2 - true_u_var = Inf*ones((N+1)*4) - true_u_var[1:3:end-(N+1)] .= x_max1 - true_u_var[2:3:end-(N+1)] .= x_max2 - @test lb[1:end-5] == true_lb - @test ub[1:end-5] == -true_lb - @test l_var == true_l_var - @test u_var == true_u_var - end + lb, ub = CTDirect.constraints_bounds(ctd) + l_var, u_var = CTDirect.variables_bounds(ctd) + true_lb = zeros(3 * N) # test without boundary conditions because of the dictionary + # true_lb[end-4:end] = [-1,0,0,0,0] + true_l_var = -Inf * ones((N + 1) * 4) + true_l_var[1:3:end-(N+1)] .= x_min1 + true_l_var[2:3:end-(N+1)] .= x_min2 + true_u_var = Inf * ones((N + 1) * 4) + true_u_var[1:3:end-(N+1)] .= x_max1 + true_u_var[2:3:end-(N+1)] .= x_max2 + @test lb[1:end-5] == true_lb + @test ub[1:end-5] == -true_lb + @test l_var == true_l_var + @test u_var == true_u_var + end end - - diff --git a/test/old/test_scalar_vector_CTBase.jl b/test/old/test_scalar_vector_CTBase.jl index dbfd9b92..7bc48a48 100644 --- a/test/old/test_scalar_vector_CTBase.jl +++ b/test/old/test_scalar_vector_CTBase.jl @@ -3,40 +3,40 @@ using CTBase using Test -n=1 -m=1 -t0=0 -tf=1 -x0=2 -xf=0 +n = 1 +m = 1 +t0 = 0 +tf = 1 +x0 = 2 +xf = 0 ocp = Model() state!(ocp, n) control!(ocp, m) time!(ocp, [t0, tf]) constraint!(ocp, :initial, x0, :initial_constraint) constraint!(ocp, :control, -1, 1, :control_constraint) -dynamics!(ocp, (x, u) -> -x+u) +dynamics!(ocp, (x, u) -> -x + u) objective!(ocp, :lagrange, (x, u) -> x[1]^2 + u^2) f = ocp.dynamics lag1 = ocp.lagrange v = Real[] function test_scalar_vector(f) -u_vector = false -x_vector = false -v = Real[] -if n == 1 && m == 1 - x = zeros(n) - u = zeros(m) - try - f(0,x,u,v) - u_vector = true - x_vector = true - catch - + u_vector = false + x_vector = false + v = Real[] + if n == 1 && m == 1 + x = zeros(n) + u = zeros(m) + try + f(0, x, u, v) + u_vector = true + x_vector = true + catch + + end end -end -return x_vector, u_vector + return x_vector, u_vector end x_vector, u_vector = test_scalar_vector(f) @@ -44,44 +44,41 @@ println("u_vector = ", u_vector) println("x_vector = ", x_vector) @testset verbose = true showtiming = true "x and u scalar " begin - @test f(t0,2,1,v) == -1 # pass - @test f(t0,[1,2],[1,1],v) ==[0,-1] # pass - # @test f(t0,2,[1],v) == -1 # fail - # @test f(t0,[2],1,v) == -1 # fail - # @test f(t0,[2],[1],v) == -1 # fail - @test f(t0,[2],[1],v) == [-1] # pass - # @test f(t0,[2,4],[1],v) == -1 # fail - # @test f(t0,x0[1:1],[1],v) == -1 # fail - # @test lag(t0,x0[1],1,v) == 5 # fail -end + @test f(t0, 2, 1, v) == -1 # pass + @test f(t0, [1, 2], [1, 1], v) == [0, -1] # pass + # @test f(t0,2,[1],v) == -1 # fail + # @test f(t0,[2],1,v) == -1 # fail + # @test f(t0,[2],[1],v) == -1 # fail + @test f(t0, [2], [1], v) == [-1] # pass + # @test f(t0,[2,4],[1],v) == -1 # fail + # @test f(t0,x0[1:1],[1],v) == -1 # fail + # @test lag(t0,x0[1],1,v) == 5 # fail +end -n=1 -m=1 -t0=0 -tf=1 -x0=[2] -xf=[0] +n = 1 +m = 1 +t0 = 0 +tf = 1 +x0 = [2] +xf = [0] ocp = Model() state!(ocp, n) control!(ocp, m) time!(ocp, [t0, tf]) constraint!(ocp, :initial, x0, :initial_constraint) constraint!(ocp, :control, -1, 1, :control_constraint) -dynamics!(ocp, (x, u) -> -x[1]+u[1]) +dynamics!(ocp, (x, u) -> -x[1] + u[1]) objective!(ocp, :lagrange, (x, u) -> x[1]^2 + u^2) f_vec = ocp.dynamics v = Real[] -f_vec(t0,x0,1,v) +f_vec(t0, x0, 1, v) lag_scal_vec = ocp.lagrange @testset verbose = true showtiming = true "x and u vector of dim 1 " begin - @test f_vec(t0,2,1,v) == -1 - @test f_vec(t0,2,[1],v) == -1 - @test f_vec(t0,[2],1,v) == -1 - @test f_vec(t0,[2],[1],v) == -1 - @test f_vec(t0,[2,4],[1],v) == -1 - @test f_vec(t0,x0[1:1],[1],v) == -1 - @test lag_scal_vec(t0,x0[1:1],1,v) == 5 -end - - - + @test f_vec(t0, 2, 1, v) == -1 + @test f_vec(t0, 2, [1], v) == -1 + @test f_vec(t0, [2], 1, v) == -1 + @test f_vec(t0, [2], [1], v) == -1 + @test f_vec(t0, [2, 4], [1], v) == -1 + @test f_vec(t0, x0[1:1], [1], v) == -1 + @test lag_scal_vec(t0, x0[1:1], 1, v) == 5 +end diff --git a/test/old/test_simple_integrator_energy.jl b/test/old/test_simple_integrator_energy.jl index 612b2e1d..f16a0e7c 100644 --- a/test/old/test_simple_integrator_energy.jl +++ b/test/old/test_simple_integrator_energy.jl @@ -1,18 +1,18 @@ # exponential, dim1, energy println(":exponential, :dim1, :energy problem") -prob = Problem(:exponential, :dim1, :energy) +prob = Problem(:exponential, :dim1, :energy) ocp = prob.model u_sol(t) = prob.solution.control(t)[1] # solve println("Is solvable ? ", CTDirect.is_solvable(ocp)) -sol = solve(ocp, grid_size=100, print_level=0) +sol = solve(ocp, grid_size = 100, print_level = 0) # check solution u = t -> sol.control(t)[1] T = sol.times -dT = T[2:end]-T[1:end-1] +dT = T[2:end] - T[1:end-1] N = length(T) -@test sum(dT .* [ abs(u(T[i])-u_sol(T[i])) for i ∈ 1:N-1] ) ≈ 0 atol=1e-1 -@test sol.objective ≈ prob.solution.objective atol=1e-2 +@test sum(dT .* [abs(u(T[i]) - u_sol(T[i])) for i ∈ 1:N-1]) ≈ 0 atol = 1e-1 +@test sol.objective ≈ prob.solution.objective atol = 1e-2 # @test constraints_violation(sol) < 1e-6 # ceci n'existe pas dans la OptimalControlSolution pour le moment diff --git a/test/old/test_utils.jl b/test/old/test_utils.jl index 6e1ca7f5..7f1cb61c 100644 --- a/test/old/test_utils.jl +++ b/test/old/test_utils.jl @@ -4,17 +4,17 @@ using LinearAlgebra Compute the infty norm : Max_{t_i}(||x(t_i)||_2) """ -function distance_infty(x::Function, x_sol::Function, T::TimesDisc) +function distance_infty(x::Function, x_sol::Function, T::TimesDisc) N = length(T) - return maximum([ norm(x(T[i])-x_sol(T[i])) for i in 1:N] ) + return maximum([norm(x(T[i]) - x_sol(T[i])) for i = 1:N]) end """ Compute the L2 norm : Int_{t_0}^{t_f}||u(t)|| dt """ -function distance_L2(u:: Function, u_sol:: Function, T::TimesDisc) +function distance_L2(u::Function, u_sol::Function, T::TimesDisc) N = length(T) - dT = T[2:end]-T[1:end-1] - return sum(dT .* [ abs(u(T[i])-u_sol(T[i])) for i ∈ 1:N-1] ) + dT = T[2:end] - T[1:end-1] + return sum(dT .* [abs(u(T[i]) - u_sol(T[i])) for i ∈ 1:N-1]) end diff --git a/test/problems/beam.jl b/test/problems/beam.jl index ade2c6ab..ad2b10e7 100644 --- a/test/problems/beam.jl +++ b/test/problems/beam.jl @@ -3,16 +3,16 @@ function beam() @def beam begin - t ∈ [ 0, 1 ], time + t ∈ [0, 1], time x ∈ R², state u ∈ R, control - x(0) == [0,1] - x(1) == [0,-1] + x(0) == [0, 1] + x(1) == [0, -1] ẋ(t) == [x₂(t), u(t)] - 0 ≤ x₁(t) ≤ 0.1 + 0 ≤ x₁(t) ≤ 0.1 -10 ≤ u(t) ≤ 10 ∫(u(t)^2) → min end - return ((ocp=beam, obj=8.898598, name="beam", init=nothing)) -end \ No newline at end of file + return ((ocp = beam, obj = 8.898598, name = "beam", init = nothing)) +end diff --git a/test/problems/bioreactor.jl b/test/problems/bioreactor.jl index 1e33c4f2..79007fa0 100644 --- a/test/problems/bioreactor.jl +++ b/test/problems/bioreactor.jl @@ -7,16 +7,16 @@ # growth model MONOD function growth(s, mu2m, Ks) - return mu2m * s / (s+Ks) + return mu2m * s / (s + Ks) end # light model: max^2 (0,sin) * mubar # DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi] function light(time, halfperiod) pi = 3.141592653589793 - days = time / (halfperiod*2) - tau = (days - floor(days)) * 2*pi - return max(0,sin(tau))^2 + days = time / (halfperiod * 2) + tau = (days - floor(days)) * 2 * pi + return max(0, sin(tau))^2 end # 1 day periodic problem @@ -34,7 +34,7 @@ function bioreactor_1day() T = halfperiod * 2 # ocp - t ∈ [ 0, T ], time + t ∈ [0, T], time x ∈ R³, state u ∈ R, control y = x[1] @@ -42,18 +42,20 @@ function bioreactor_1day() b = x[3] mu = light(t, halfperiod) * mubar mu2 = growth(s(t), mu2m, Ks) - [0,0,0.001] ≤ x(t) ≤ [Inf, Inf, Inf] + [0, 0, 0.001] ≤ x(t) ≤ [Inf, Inf, Inf] 0 ≤ u(t) ≤ 1 1 ≤ y(0) ≤ Inf 1 ≤ b(0) ≤ Inf - x(0)- x(T) == [0,0,0] - ẋ(t) == [mu*y(t)/(1+y(t))-(r+u(t))*y(t), - -mu2*b(t) + u(t)*beta*(gamma*y(t)-s(t)), - (mu2-u(t)*beta)*b(t)] - ∫(mu2*b(t)/(beta+c)) → max + x(0) - x(T) == [0, 0, 0] + ẋ(t) == [ + mu * y(t) / (1 + y(t)) - (r + u(t)) * y(t), + -mu2 * b(t) + u(t) * beta * (gamma * y(t) - s(t)), + (mu2 - u(t) * beta) * b(t), + ] + ∫(mu2 * b(t) / (beta + c)) → max end - return ((ocp=bioreactor_1, obj=0.614134, name="bioreactor_1day", init=nothing)) + return ((ocp = bioreactor_1, obj = 0.614134, name = "bioreactor_1day", init = nothing)) end @@ -73,7 +75,7 @@ function bioreactor_Ndays() T = 10 * N # ocp - t ∈ [ 0, T ], time + t ∈ [0, T], time x ∈ R³, state u ∈ R, control y = x[1] @@ -81,16 +83,23 @@ function bioreactor_Ndays() b = x[3] mu = light(t, halfperiod) * mubar mu2 = growth(s(t), mu2m, Ks) - [0,0,0.001] ≤ x(t) ≤ [Inf, Inf, Inf] + [0, 0, 0.001] ≤ x(t) ≤ [Inf, Inf, Inf] 0 ≤ u(t) ≤ 1 0.05 ≤ y(0) ≤ 0.25 0.5 ≤ s(0) ≤ 5 0.5 ≤ b(0) ≤ 3 - ẋ(t) == [mu*y(t)/(1+y(t))-(r+u(t))*y(t), - -mu2*b(t) + u(t)*beta*(gamma*y(t)-s(t)), - (mu2-u(t)*beta)*b(t)] - ∫(mu2*b(t)/(beta+c)) → max + ẋ(t) == [ + mu * y(t) / (1 + y(t)) - (r + u(t)) * y(t), + -mu2 * b(t) + u(t) * beta * (gamma * y(t) - s(t)), + (mu2 - u(t) * beta) * b(t), + ] + ∫(mu2 * b(t) / (beta + c)) → max end - return ((ocp=bioreactor_N, obj=19.0745, init=(state=[50,50,50],), name="bioreactor_Ndays")) -end \ No newline at end of file + return (( + ocp = bioreactor_N, + obj = 19.0745, + init = (state = [50, 50, 50],), + name = "bioreactor_Ndays", + )) +end diff --git a/test/problems/bolza.jl b/test/problems/bolza.jl index deacb853..3cf06edd 100644 --- a/test/problems/bolza.jl +++ b/test/problems/bolza.jl @@ -12,6 +12,5 @@ function bolza_freetf() tf + 0.5∫(u(t)^2) → min end - return ((ocp=ocp, obj=1.476, name="bolza_freetf", init=nothing)) + return ((ocp = ocp, obj = 1.476, name = "bolza_freetf", init = nothing)) end - diff --git a/test/problems/double_integrator.jl b/test/problems/double_integrator.jl index 4b5f49a5..2a8bc683 100644 --- a/test/problems/double_integrator.jl +++ b/test/problems/double_integrator.jl @@ -1,27 +1,27 @@ # double integrator function double_integrator_a() - + @def ocp begin tf ∈ R, variable - t ∈ [ 0, tf ], time + t ∈ [0, tf], time x ∈ R², state u ∈ R, control -1 ≤ u(t) ≤ 1 - x(0) == [ 0, 0 ] - x(tf) == [ 1, 0 ] - 0.05 ≤ tf ≤ Inf - ẋ(t) == [ x₂(t), u(t) ] + x(0) == [0, 0] + x(tf) == [1, 0] + 0.05 ≤ tf ≤ Inf + ẋ(t) == [x₂(t), u(t)] tf → min end - - return ((ocp=ocp, obj=2.0, name="double_integrator_a", init=nothing)) + + return ((ocp = ocp, obj = 2.0, name = "double_integrator_a", init = nothing)) end function double_integrator_T(T) @def ocp begin - t ∈ [ 0, T ], time + t ∈ [0, T], time x ∈ R², state u ∈ R, control q = x₁ @@ -30,26 +30,26 @@ function double_integrator_T(T) v(0) == 0 q(T) == 1 v(T) == 0 - ẋ(t) == [ v(t), u(t) ] + ẋ(t) == [v(t), u(t)] ∫(u(t)^2) → min end - return ((ocp=ocp, obj=nothing, name="double_integrator_T", init=nothing)) + return ((ocp = ocp, obj = nothing, name = "double_integrator_T", init = nothing)) end # min tf -function double_integrator_mintf(;lagrange=false) +function double_integrator_mintf(; lagrange = false) - ocp = Model(variable=true) + ocp = Model(variable = true) state!(ocp, 2) control!(ocp, 1) variable!(ocp, 1) - time!(ocp, t0=0, indf=1) - constraint!(ocp, :initial, val=[0,0]) - constraint!(ocp, :final, val=[1,0]) - constraint!(ocp, :control, lb=-1, ub=1) - constraint!(ocp, :variable, lb=0.1, ub=10) - dynamics!(ocp, (x, u, v) -> [x[2], u]) + time!(ocp, t0 = 0, indf = 1) + constraint!(ocp, :initial, val = [0, 0]) + constraint!(ocp, :final, val = [1, 0]) + constraint!(ocp, :control, lb = -1, ub = 1) + constraint!(ocp, :variable, lb = 0.1, ub = 10) + dynamics!(ocp, (x, u, v) -> [x[2], u]) if lagrange objective!(ocp, :lagrange, (x, u, v) -> 1) name = "double_integrator_lagrange" @@ -58,24 +58,24 @@ function double_integrator_mintf(;lagrange=false) name = "double_integrator_mayer" end - return ((ocp=ocp, obj=2.0, name=name, init=nothing)) + return ((ocp = ocp, obj = 2.0, name = name, init = nothing)) end # max t0 with free t0,tf -function double_integrator_freet0tf(lagrange=false) +function double_integrator_freet0tf(lagrange = false) - ocp = Model(variable=true) + ocp = Model(variable = true) state!(ocp, 2) control!(ocp, 1) variable!(ocp, 2) - time!(ocp, ind0=1, indf=2) - constraint!(ocp, :initial, val=[0,0]) - constraint!(ocp, :final, val=[1,0]) - constraint!(ocp, :control, lb=-1, ub=1) - constraint!(ocp, :variable, lb=[0.1,0.1], ub=[10,10]) - constraint!(ocp, :variable, f=v->v[2]-v[1], lb=0.1, ub=Inf) - dynamics!(ocp, (x, u, v) -> [x[2], u]) + time!(ocp, ind0 = 1, indf = 2) + constraint!(ocp, :initial, val = [0, 0]) + constraint!(ocp, :final, val = [1, 0]) + constraint!(ocp, :control, lb = -1, ub = 1) + constraint!(ocp, :variable, lb = [0.1, 0.1], ub = [10, 10]) + constraint!(ocp, :variable, f = v -> v[2] - v[1], lb = 0.1, ub = Inf) + dynamics!(ocp, (x, u, v) -> [x[2], u]) objective!(ocp, :mayer, (x0, xf, v) -> v[1], :max) - return ((ocp=ocp, obj=8.0, name="double_integrator_freet0tf", init=nothing)) -end \ No newline at end of file + return ((ocp = ocp, obj = 8.0, name = "double_integrator_freet0tf", init = nothing)) +end diff --git a/test/problems/fuller.jl b/test/problems/fuller.jl index 1740f118..e5892363 100644 --- a/test/problems/fuller.jl +++ b/test/problems/fuller.jl @@ -2,15 +2,15 @@ function fuller() @def fuller begin - t ∈ [ 0, 3.5 ], time + t ∈ [0, 3.5], time x ∈ R², state u ∈ R, control -1 ≤ u(t) ≤ 1 - x(0) == [ 0, 1 ] - x(3.5) == [ 0, 0 ] - ẋ(t) == [ x₂(t), u(t) ] + x(0) == [0, 1] + x(3.5) == [0, 0] + ẋ(t) == [x₂(t), u(t)] ∫(x₁(t)^2) → min end - return((ocp=fuller, obj=2.683944e-1, name="fuller", init=nothing)) -end \ No newline at end of file + return ((ocp = fuller, obj = 2.683944e-1, name = "fuller", init = nothing)) +end diff --git a/test/problems/goddard.jl b/test/problems/goddard.jl index 386f3b5d..18e44cac 100644 --- a/test/problems/goddard.jl +++ b/test/problems/goddard.jl @@ -5,15 +5,15 @@ # aux functions function F0(x, Cd, beta) r, v, m = x - D = Cd * v^2 * exp(-beta*(r - 1)) - return [ v, -D/m - 1/r^2, 0 ] + D = Cd * v^2 * exp(-beta * (r - 1)) + return [v, -D / m - 1 / r^2, 0] end function F1(x, Tmax, b) r, v, m = x - return [ 0, Tmax/m, -b*Tmax ] + return [0, Tmax / m, -b * Tmax] end -function goddard(;vmax=0.1, Tmax=3.5, functional_constraints=false) +function goddard(; vmax = 0.1, Tmax = 3.5, functional_constraints = false) # constants Cd = 310 beta = 500 @@ -22,30 +22,41 @@ function goddard(;vmax=0.1, Tmax=3.5, functional_constraints=false) v0 = 0 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] #ocp - goddard = Model(variable=true) + goddard = Model(variable = true) state!(goddard, 3) control!(goddard, 1) variable!(goddard, 1) - time!(goddard, t0=0, indf=1) - constraint!(goddard, :initial, val=x0) - constraint!(goddard, :final, rg=3, val=mf) + time!(goddard, t0 = 0, indf = 1) + constraint!(goddard, :initial, val = x0) + constraint!(goddard, :final, rg = 3, val = mf) if functional_constraints # note: the equations do not handle r<1 well # without the box constraint on x, the default init (0.1) is not suitable - constraint!(goddard, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0]) - constraint!(goddard, :control, f=(u,v)->u, lb=0, ub=1) + constraint!( + goddard, + :state, + f = (x, v) -> x, + lb = [r0, v0, mf], + ub = [r0 + 0.2, vmax, m0], + ) + constraint!(goddard, :control, f = (u, v) -> u, lb = 0, ub = 1) else - constraint!(goddard, :state, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0]) - constraint!(goddard, :control, lb=0, ub=1) + constraint!(goddard, :state, lb = [r0, v0, mf], ub = [r0 + 0.2, vmax, m0]) + constraint!(goddard, :control, lb = 0, ub = 1) end - constraint!(goddard, :variable, lb=0.01, ub=Inf) - objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) - dynamics!(goddard, (x, u, v) -> F0(x, Cd, beta) + u*F1(x, Tmax, b) ) + constraint!(goddard, :variable, lb = 0.01, ub = Inf) + objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) + dynamics!(goddard, (x, u, v) -> F0(x, Cd, beta) + u * F1(x, Tmax, b)) - return ((ocp=goddard, obj=1.01257, name="goddard", init=(state=[1.01,0.05,0.8],))) + return (( + ocp = goddard, + obj = 1.01257, + name = "goddard", + init = (state = [1.01, 0.05, 0.8],), + )) end @@ -59,42 +70,47 @@ function goddard_all() v0 = 0 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] vmax = 0.1 Tmax = 3.5 #ocp - goddard = Model(variable=true) + goddard = Model(variable = true) state!(goddard, 3) control!(goddard, 1) variable!(goddard, 1) - time!(goddard, t0=0, indf=1) + time!(goddard, t0 = 0, indf = 1) # initial constraint - constraint!(goddard, :initial, val=x0) + constraint!(goddard, :initial, val = x0) # final constraint - constraint!(goddard, :final, rg=3, val=mf) + constraint!(goddard, :final, rg = 3, val = mf) # state box (active at t0 and tf) - constraint!(goddard, :state, lb=[r0,v0,0], ub=[Inf,Inf,m0]) + constraint!(goddard, :state, lb = [r0, v0, 0], ub = [Inf, Inf, m0]) # control box (active on last bang arc) - constraint!(goddard, :control, lb=0, ub=Inf) + constraint!(goddard, :control, lb = 0, ub = Inf) # variable box (inactive) - constraint!(goddard, :variable, lb=0.01, ub=Inf) + constraint!(goddard, :variable, lb = 0.01, ub = Inf) # state constraint (active on constrained arc) - constraint!(goddard, :state, f=(x,v)->x[2], lb=-Inf, ub=vmax) + constraint!(goddard, :state, f = (x, v) -> x[2], lb = -Inf, ub = vmax) # control constraint (active on first bang arc) - constraint!(goddard, :control, f=(u,v)->u, lb=-Inf, ub=1) + constraint!(goddard, :control, f = (u, v) -> u, lb = -Inf, ub = 1) # 'mixed' constraint (active at tf) - constraint!(goddard, :mixed, f=(x,u,v)->x[3], lb=mf, ub=Inf) - objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) - dynamics!(goddard, (x, u, v) -> F0(x, Cd, beta) + u*F1(x, Tmax, b) ) + constraint!(goddard, :mixed, f = (x, u, v) -> x[3], lb = mf, ub = Inf) + objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) + dynamics!(goddard, (x, u, v) -> F0(x, Cd, beta) + u * F1(x, Tmax, b)) - return ((ocp=goddard, obj=1.01257, name="goddard_all", init=(state=[1.01,0.05,0.8],))) + return (( + ocp = goddard, + obj = 1.01257, + name = "goddard_all", + init = (state = [1.01, 0.05, 0.8],), + )) end # abstract definition -function goddard_a(;vmax=0.1, Tmax=3.5) +function goddard_a(; vmax = 0.1, Tmax = 3.5) # constants Cd = 310 beta = 500 @@ -103,11 +119,11 @@ function goddard_a(;vmax=0.1, Tmax=3.5) v0 = 0 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] @def goddard_a begin tf ∈ R, variable - t ∈ [ 0, tf ], time + t ∈ [0, tf], time x ∈ R^3, state u ∈ R, control 0.1 ≤ tf ≤ Inf @@ -120,9 +136,14 @@ function goddard_a(;vmax=0.1, Tmax=3.5) v0 ≤ v(t) ≤ 0.1 mf ≤ m(t) ≤ m0 0 ≤ u(t) ≤ 1 - ẋ(t) == F0(x(t), Cd, beta) + u(t)*F1(x(t), Tmax, b) + ẋ(t) == F0(x(t), Cd, beta) + u(t) * F1(x(t), Tmax, b) r(tf) → max end - return ((ocp=goddard_a, obj=1.01257, name="goddard_a", init=(state=[1.01,0.05,0.8],))) -end \ No newline at end of file + return (( + ocp = goddard_a, + obj = 1.01257, + name = "goddard_a", + init = (state = [1.01, 0.05, 0.8],), + )) +end diff --git a/test/problems/insurance.jl b/test/problems/insurance.jl index 4bbcac0d..c9628cd9 100644 --- a/test/problems/insurance.jl +++ b/test/problems/insurance.jl @@ -2,7 +2,7 @@ function insurance() @def insurance begin - + # constants gamma = 0.2 lambda = 0.25 @@ -14,7 +14,7 @@ function insurance() alpha = 4 tf = 10 - t ∈ [ 0, tf ], time + t ∈ [0, tf], time x ∈ R³, state u ∈ R⁵, control P ∈ R, variable @@ -27,7 +27,7 @@ function insurance() dUdR = u[5] 0 ≤ I(t) ≤ 1.5 - 0 ≤ m(t) ≤ 1.5 + 0 ≤ m(t) ≤ 1.5 0 ≤ h(t) ≤ 25 0 ≤ R(t) ≤ Inf 0 ≤ H(t) ≤ Inf @@ -40,21 +40,19 @@ function insurance() epsilon = k * t / (tf - t + 1) # illness distribution - fx = lambda*exp(-lambda*t) + exp(-lambda*tf)/tf + fx = lambda * exp(-lambda * t) + exp(-lambda * tf) / tf # expense effect - v = m(t)^(alpha/2) / (1+m(t)^(alpha/2)) - vprime = alpha/2 * m(t)^(alpha/2-1) / (1+m(t)^(alpha/2))^2 + v = m(t)^(alpha / 2) / (1 + m(t)^(alpha / 2)) + vprime = alpha / 2 * m(t)^(alpha / 2 - 1) / (1 + m(t)^(alpha / 2))^2 R(t) - (w - P + I(t) - m(t) - epsilon) == 0 H(t) - (h0 - gamma * t * (1 - v)) == 0 - U(t) - (1 - exp( - s * R(t))+ H(t)) == 0 - dUdR(t) - (s * exp( - s * R(t))) == 0 + U(t) - (1 - exp(-s * R(t)) + H(t)) == 0 + dUdR(t) - (s * exp(-s * R(t))) == 0 - ẋ(t) == [(1-gamma*t*vprime/dUdR(t))*h(t), - h(t), - (1+sigma)*I(t)*fx] - ∫(U(t)*fx) → max + ẋ(t) == [(1 - gamma * t * vprime / dUdR(t)) * h(t), h(t), (1 + sigma) * I(t) * fx] + ∫(U(t) * fx) → max end - return((ocp=insurance, obj=2.059511, name="insurance", init=nothing)) -end \ No newline at end of file + return ((ocp = insurance, obj = 2.059511, name = "insurance", init = nothing)) +end diff --git a/test/problems/jackson.jl b/test/problems/jackson.jl index 10421eb8..18e7f28f 100644 --- a/test/problems/jackson.jl +++ b/test/problems/jackson.jl @@ -7,19 +7,21 @@ function jackson() k2 = 10 k3 = 1 - t ∈ [ 0, 4 ], time + t ∈ [0, 4], time x ∈ R³, state u ∈ R, control - [0, 0, 0] ≤ x(t) ≤ [1.1, 1.1, 1.1] + [0, 0, 0] ≤ x(t) ≤ [1.1, 1.1, 1.1] 0 ≤ u(t) ≤ 1 - x(0) == [ 1, 0, 0 ] + x(0) == [1, 0, 0] a = x[1] b = x[2] - ẋ(t) == [-u(t)*(k1*a(t)-k2*b(t)), - u(t)*(k1*a(t)-k2*b(t)) - (1-u(t))*k3*b(t), - (1-u(t))*k3*b(t)] + ẋ(t) == [ + -u(t) * (k1 * a(t) - k2 * b(t)), + u(t) * (k1 * a(t) - k2 * b(t)) - (1 - u(t)) * k3 * b(t), + (1 - u(t)) * k3 * b(t), + ] x[3](4) → max end - return ((ocp=jackson, obj=0.192011, name="jackson", init=nothing)) -end \ No newline at end of file + return ((ocp = jackson, obj = 0.192011, name = "jackson", init = nothing)) +end diff --git a/test/problems/parametric.jl b/test/problems/parametric.jl index 3349f1df..43e2e4ae 100644 --- a/test/problems/parametric.jl +++ b/test/problems/parametric.jl @@ -1,27 +1,27 @@ # Parametric problem (name ??) function parametric(ρ) - + relu(x) = max(0, x) μ = 10 - p_relu(x) = log(abs(1 + exp(μ*x)))/μ - f(x) = 1-x - m(x) = (p_relu∘f)(x) + p_relu(x) = log(abs(1 + exp(μ * x))) / μ + f(x) = 1 - x + m(x) = (p_relu ∘ f)(x) T = 2 @def param begin τ ∈ R, variable - s ∈ [ 0, 1 ], time + s ∈ [0, 1], time x ∈ R², state u ∈ R², control x₁(0) == 0 x₂(0) == 1 x₁(1) == 1 - ẋ(s) == [τ*(u₁(s)+2), (T-τ)*u₂(s)] + ẋ(s) == [τ * (u₁(s) + 2), (T - τ) * u₂(s)] -1 ≤ u₁(s) ≤ 1 -1 ≤ u₂(s) ≤ 1 0 ≤ τ ≤ T - -(x₂(1)-2)^3 - ∫( ρ * ( τ*m(x₁(s))^2 + (T-τ)*m(x₂(s))^2 ) ) → min + -(x₂(1) - 2)^3 - ∫(ρ * (τ * m(x₁(s))^2 + (T - τ) * m(x₂(s))^2)) → min end - return ((ocp=param, obj=nothing, name="parametric", init=nothing)) -end \ No newline at end of file + return ((ocp = param, obj = nothing, name = "parametric", init = nothing)) +end diff --git a/test/problems/robbins.jl b/test/problems/robbins.jl index 73cd18ea..6c53191d 100644 --- a/test/problems/robbins.jl +++ b/test/problems/robbins.jl @@ -6,15 +6,15 @@ function robbins() beta = 0 gamma = 0.5 - t ∈ [ 0, 10 ], time + t ∈ [0, 10], time x ∈ R³, state u ∈ R, control 0 ≤ x[1](t) ≤ Inf - x(0) == [ 1, -2, 0 ] - x(10) == [ 0, 0, 0 ] - ẋ(t) == [ x[2](t), x[3](t), u(t) ] - ∫(alpha*x[1](t) + beta*x[1](t)^2 + gamma*u(t)^2) → min + x(0) == [1, -2, 0] + x(10) == [0, 0, 0] + ẋ(t) == [x[2](t), x[3](t), u(t)] + ∫(alpha * x[1](t) + beta * x[1](t)^2 + gamma * u(t)^2) → min end - return ((ocp=robbins, obj=20., name="robbins", init=nothing)) -end \ No newline at end of file + return ((ocp = robbins, obj = 20.0, name = "robbins", init = nothing)) +end diff --git a/test/problems/simple_integrator.jl b/test/problems/simple_integrator.jl index 065d18df..426961b5 100644 --- a/test/problems/simple_integrator.jl +++ b/test/problems/simple_integrator.jl @@ -6,12 +6,12 @@ function simple_integrator() ocp = Model() state!(ocp, 1) control!(ocp, 2) - time!(ocp, t0=0, tf=1) - constraint!(ocp, :initial, val=-1) - constraint!(ocp, :final, val=0) - constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) + time!(ocp, t0 = 0, tf = 1) + constraint!(ocp, :initial, val = -1) + constraint!(ocp, :final, val = 0) + constraint!(ocp, :control, lb = [0, 0], ub = [Inf, Inf]) dynamics!(ocp, (x, u) -> -x - u[1] + u[2]) - objective!(ocp, :lagrange, (x, u) -> (u[1]+u[2])^2) + objective!(ocp, :lagrange, (x, u) -> (u[1] + u[2])^2) - return ((ocp=ocp, obj=nothing, name="simple_integrator", init=nothing)) -end \ No newline at end of file + return ((ocp = ocp, obj = nothing, name = "simple_integrator", init = nothing)) +end diff --git a/test/problems/swimmer.jl b/test/problems/swimmer.jl index 3f7fc1a6..8fae1718 100644 --- a/test/problems/swimmer.jl +++ b/test/problems/swimmer.jl @@ -3,19 +3,19 @@ # +++ make 2 versions: 1 stroke periodic and free N strokes function swimmer() - + @def swimmer begin tf = 25 - t ∈ [ 0, tf ], time + t ∈ [0, tf], time x ∈ R^5, state u ∈ R^2, control # bounds -3.15 ≤ x[3](t) ≤ 3.15 [-1.5, -1.5] ≤ x[4:5](t) ≤ [1.5, 1.5] - [-1, -1] ≤ u(t) ≤ [1, 1] - + [-1, -1] ≤ u(t) ≤ [1, 1] + # initial conditions x[1:2](0) == [0, 0] -3.15 ≤ x[3](0) ≤ 0 @@ -36,28 +36,107 @@ function swimmer() u1 = u[1](t) u2 = u[2](t) - aux = 543 + 186*cos(b1) + 37*cos(2*b1) + 12*cos(b1 - 2*b3) + 30*cos(b1 - b3) + 2*cos(2*(b1 - b3)) + 12*cos(2*b1 - b3) + 186*cos(b3) + 37*cos(2*b3) - 6*cos(b1 + b3) - 3*cos(2*(b1 + b3)) - 6*cos(2*b1 + b3) - 6*cos(b1 + 2*b3) - - g11 = (-42*sin(b1 - th) - 2*sin(2*b1 - th) - 24*sin(th) - 300*sin(b1 + th) - 12*sin(2*b1 + th) - 6*sin(b1 - th - 2*b3) - sin(2*b1 - th - 2*b3) + 4*sin(th - 2*b3) - 12*sin(b1 + th - 2*b3) - sin(2*b1 + th - 2*b3) + 18*sin(b1 - th - b3) + 8*sin(th - b3) - 54*sin(b1 + th - b3) - 2*sin(2*b1 + th - b3) - 18*sin(b1 - th + b3) - 38*sin(th + b3) - 90*sin(b1 + th + b3) - 6*sin(b1 - th + 2*b3) - 18*sin(th + 2*b3) - 30*sin(b1 + th + 2*b3)) / (4*aux) - - g12 = (-42*cos(b1 - th) - 2*cos(2*b1 - th) + 24*cos(th) + 300*cos(b1 + th) + 12*cos(2*b1 + th) - 6*cos(b1 - th - 2*b3) - cos(2*b1 - th - 2*b3) - 4*cos(th - 2*b3) + 12*cos(b1 + th - 2*b3) + cos(2*b1 + th - 2*b3) + 18*cos(b1 - th - b3) - 8*cos(th - b3) + 54*cos(b1 + th - b3)+ 2*cos(2*b1 + th - b3) - 18*cos(b1 - th + b3) + 38*cos(th + b3) + 90*cos(b1 + th + b3) - 6*cos(b1 - th + 2*b3) + 18*cos(th + 2*b3) + 30*cos(b1 + th + 2*b3)) / (4*aux) - - g13 = -(105 + 186*cos(b1) + 2*cos(2*b1) + 12*cos(b1 - 2*b3) + 30*cos(b1 - b3) + cos(2*(b1 - b3)) - 4*cos(2*b3) - 6*cos(b1 + b3) - 6*cos(b1 + 2*b3)) / (2*aux) - - g21 = (8*sin(b1 - th) + 4*sin(2*b1 - th) + 24*sin(th) + 38*sin(b1 + th) + 18*sin(2*b1 + th) - 2*sin(b1 - th - 2*b3) - sin(2*b1 - th - 2*b3) - 2*sin(th - 2*b3) - sin(2*b1 + th - 2*b3) - 54*sin(b1 - th - b3) - 12*sin(2*b1 - th - b3) - 42*sin(th - b3) + 18*sin(b1 + th - b3) - 6*sin(2*b1 + th - b3) + 18*sin(b1 - th + b3) + 6*sin(2*b1 - th + b3) + 300*sin(th + b3) + 90*sin(b1 + th + b3) + 30*sin(2*b1 + th + b3) + 12*sin(th + 2*b3)) / (4*aux) - - g22 = (8*cos(b1 - th) + 4*cos(2*b1 - th) - 24*cos(th) - 38*cos(b1 + th) - 18*cos(2*b1 + th) - 2*cos(b1 - th - 2*b3) - cos(2*b1 - th - 2*b3) + 2*cos(th - 2*b3) + cos(2*b1 + th - 2*b3) - 54*cos(b1 - th - b3) - 12*cos(2*b1 - th - b3) + 42*cos(th - b3) - 18*cos(b1 + th - b3) + 6*cos(2*b1 + th - b3) + 18*cos(b1 - th + b3) + 6*cos(2*b1 - th + b3) - 300*cos(th + b3) - 90*cos(b1 + th + b3) - 30*cos(2*b1 + th + b3) - 12*cos(th + 2*b3)) / (4*aux) - - g23 = -(105 - 4*cos(2*b1) + 30*cos(b1 - b3) + cos(2*(b1 - b3)) + 12*cos(2*b1 - b3) + 186*cos(b3) + 2*cos(2*b3) - 6*cos(b1 + b3) - 6*cos(2*b1 + b3)) / (2*aux) - - ẋ(t) == [g11*u1 + g21*u2, - g12*u1 + g22*u2, - g13*u1 + g23*u2, - u1, - u2] + aux = + 543 + + 186 * cos(b1) + + 37 * cos(2 * b1) + + 12 * cos(b1 - 2 * b3) + + 30 * cos(b1 - b3) + + 2 * cos(2 * (b1 - b3)) + + 12 * cos(2 * b1 - b3) + + 186 * cos(b3) + + 37 * cos(2 * b3) - 6 * cos(b1 + b3) - 3 * cos(2 * (b1 + b3)) - + 6 * cos(2 * b1 + b3) - 6 * cos(b1 + 2 * b3) + + g11 = + ( + -42 * sin(b1 - th) - 2 * sin(2 * b1 - th) - 24 * sin(th) - + 300 * sin(b1 + th) - 12 * sin(2 * b1 + th) - 6 * sin(b1 - th - 2 * b3) - + sin(2 * b1 - th - 2 * b3) + 4 * sin(th - 2 * b3) - + 12 * sin(b1 + th - 2 * b3) - sin(2 * b1 + th - 2 * b3) + + 18 * sin(b1 - th - b3) + + 8 * sin(th - b3) - 54 * sin(b1 + th - b3) - 2 * sin(2 * b1 + th - b3) - + 18 * sin(b1 - th + b3) - 38 * sin(th + b3) - 90 * sin(b1 + th + b3) - + 6 * sin(b1 - th + 2 * b3) - 18 * sin(th + 2 * b3) - + 30 * sin(b1 + th + 2 * b3) + ) / (4 * aux) + + g12 = + ( + -42 * cos(b1 - th) - 2 * cos(2 * b1 - th) + + 24 * cos(th) + + 300 * cos(b1 + th) + + 12 * cos(2 * b1 + th) - 6 * cos(b1 - th - 2 * b3) - + cos(2 * b1 - th - 2 * b3) - 4 * cos(th - 2 * b3) + + 12 * cos(b1 + th - 2 * b3) + + cos(2 * b1 + th - 2 * b3) + + 18 * cos(b1 - th - b3) - 8 * cos(th - b3) + + 54 * cos(b1 + th - b3) + + 2 * cos(2 * b1 + th - b3) - 18 * cos(b1 - th + b3) + + 38 * cos(th + b3) + + 90 * cos(b1 + th + b3) - 6 * cos(b1 - th + 2 * b3) + + 18 * cos(th + 2 * b3) + + 30 * cos(b1 + th + 2 * b3) + ) / (4 * aux) + + g13 = + -( + 105 + + 186 * cos(b1) + + 2 * cos(2 * b1) + + 12 * cos(b1 - 2 * b3) + + 30 * cos(b1 - b3) + + cos(2 * (b1 - b3)) - 4 * cos(2 * b3) - 6 * cos(b1 + b3) - + 6 * cos(b1 + 2 * b3) + ) / (2 * aux) + + g21 = + ( + 8 * sin(b1 - th) + + 4 * sin(2 * b1 - th) + + 24 * sin(th) + + 38 * sin(b1 + th) + + 18 * sin(2 * b1 + th) - 2 * sin(b1 - th - 2 * b3) - + sin(2 * b1 - th - 2 * b3) - 2 * sin(th - 2 * b3) - + sin(2 * b1 + th - 2 * b3) - 54 * sin(b1 - th - b3) - + 12 * sin(2 * b1 - th - b3) - 42 * sin(th - b3) + 18 * sin(b1 + th - b3) - + 6 * sin(2 * b1 + th - b3) + + 18 * sin(b1 - th + b3) + + 6 * sin(2 * b1 - th + b3) + + 300 * sin(th + b3) + + 90 * sin(b1 + th + b3) + + 30 * sin(2 * b1 + th + b3) + + 12 * sin(th + 2 * b3) + ) / (4 * aux) + + g22 = + ( + 8 * cos(b1 - th) + 4 * cos(2 * b1 - th) - 24 * cos(th) - 38 * cos(b1 + th) - + 18 * cos(2 * b1 + th) - 2 * cos(b1 - th - 2 * b3) - + cos(2 * b1 - th - 2 * b3) + + 2 * cos(th - 2 * b3) + + cos(2 * b1 + th - 2 * b3) - 54 * cos(b1 - th - b3) - + 12 * cos(2 * b1 - th - b3) + 42 * cos(th - b3) - 18 * cos(b1 + th - b3) + + 6 * cos(2 * b1 + th - b3) + + 18 * cos(b1 - th + b3) + + 6 * cos(2 * b1 - th + b3) - 300 * cos(th + b3) - 90 * cos(b1 + th + b3) - + 30 * cos(2 * b1 + th + b3) - 12 * cos(th + 2 * b3) + ) / (4 * aux) + + g23 = + -( + 105 - 4 * cos(2 * b1) + + 30 * cos(b1 - b3) + + cos(2 * (b1 - b3)) + + 12 * cos(2 * b1 - b3) + + 186 * cos(b3) + + 2 * cos(2 * b3) - 6 * cos(b1 + b3) - 6 * cos(2 * b1 + b3) + ) / (2 * aux) + + ẋ(t) == [g11 * u1 + g21 * u2, g12 * u1 + g22 * u2, g13 * u1 + g23 * u2, u1, u2] x[1](tf) → max #∫(u1^2 + u2^2) → min end - return((ocp=swimmer, obj=0.984273, name="swimmer", init=nothing)) -end \ No newline at end of file + return ((ocp = swimmer, obj = 0.984273, name = "swimmer", init = nothing)) +end diff --git a/test/problems/vanderpol.jl b/test/problems/vanderpol.jl index 1f82facd..5f07a383 100644 --- a/test/problems/vanderpol.jl +++ b/test/problems/vanderpol.jl @@ -6,14 +6,14 @@ function vanderpol() omega = 1 epsilon = 1 - t ∈ [ 0, 2 ], time + t ∈ [0, 2], time x ∈ R², state u ∈ R, control - x(0) == [ 1, 0 ] - ẋ(t) == [ x[2](t), - epsilon*omega*(1-x[1](t)^2)*x[2](t) - omega^2*x[1](t) + u(t) ] - ∫(0.5*(x[1](t)^2 + x[2](t)^2 + u(t)^2)) → min + x(0) == [1, 0] + ẋ(t) == + [x[2](t), epsilon * omega * (1 - x[1](t)^2) * x[2](t) - omega^2 * x[1](t) + u(t)] + ∫(0.5 * (x[1](t)^2 + x[2](t)^2 + u(t)^2)) → min end - return((ocp=vanderpol, obj=1.047921, name="vanderpol", init=nothing)) -end \ No newline at end of file + return ((ocp = vanderpol, obj = 1.047921, name = "vanderpol", init = nothing)) +end diff --git a/test/refine_grid.jl b/test/refine_grid.jl index 2113d753..1045ddb7 100644 --- a/test/refine_grid.jl +++ b/test/refine_grid.jl @@ -17,26 +17,54 @@ if test4 # basic solve for comparison N = N_target println("\nBasic solve") - @time begin sol = solve(goddard, time_grid=LinRange(0,1,N+1), display=false) - @printf("steps %4d, objective %9.6f, iterations %4d\n", N, sol.objective, sol.iterations) + @time begin + sol = solve(goddard, time_grid = LinRange(0, 1, N + 1), display = false) + @printf( + "steps %4d, objective %9.6f, iterations %4d\n", + N, + sol.objective, + sol.iterations + ) end # init println("\nStep continuation") @time begin - N = 10 - sol = solve(goddard, time_grid=LinRange(0,1,N+1), display=false, max_iter=5) - @printf("steps %4d, objective %9.6f, iterations set to %4d\n", N, sol.objective, sol.iterations) - # loop on steps - while N < N_target - global N = min(N * 2, N_target) - global sol = solve(goddard, time_grid=LinRange(0,1,N+1), display=false, init=sol, max_iter=5) - @printf("steps %4d, objective %9.6f, iterations set to %4d\n", N, sol.objective, sol.iterations) - end - # final solve - N = N_target - sol = solve(goddard, time_grid=LinRange(0,1,N+1), display=false, init=sol) - @printf("steps %4d, objective %9.6f, iterations %4d\n", N, sol.objective, sol.iterations) + N = 10 + sol = + solve(goddard, time_grid = LinRange(0, 1, N + 1), display = false, max_iter = 5) + @printf( + "steps %4d, objective %9.6f, iterations set to %4d\n", + N, + sol.objective, + sol.iterations + ) + # loop on steps + while N < N_target + global N = min(N * 2, N_target) + global sol = solve( + goddard, + time_grid = LinRange(0, 1, N + 1), + display = false, + init = sol, + max_iter = 5, + ) + @printf( + "steps %4d, objective %9.6f, iterations set to %4d\n", + N, + sol.objective, + sol.iterations + ) + end + # final solve + N = N_target + sol = solve(goddard, time_grid = LinRange(0, 1, N + 1), display = false, init = sol) + @printf( + "steps %4d, objective %9.6f, iterations %4d\n", + N, + sol.objective, + sol.iterations + ) end end @@ -44,7 +72,7 @@ end # NB The function below is a dirty hack. # How is this handled by AD ? -function dt(t,v) +function dt(t, v) if t == 0 dt = v[1] elseif t == N_vars @@ -59,16 +87,16 @@ if test5 # number of time steps to be optimized N_vars = 10 - ocp = Model(variable=true, autonomous=false) + ocp = Model(variable = true, autonomous = false) state!(ocp, 2) control!(ocp, 1) variable!(ocp, N_vars) - time!(ocp, t0=0, tf=N_vars) - constraint!(ocp, :initial, lb=[0,0], ub=[0,0]) - constraint!(ocp, :final, lb=[1,0], ub=[1,0]) - constraint!(ocp, :control, lb=-1, ub=1) - constraint!(ocp, :variable, lb=0.01*ones(N_vars), ub=10*ones(N_vars)) - dynamics!(ocp, (t,x,u,v)-> [x[2], u] * dt(t,v)) + time!(ocp, t0 = 0, tf = N_vars) + constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0]) + constraint!(ocp, :final, lb = [1, 0], ub = [1, 0]) + constraint!(ocp, :control, lb = -1, ub = 1) + constraint!(ocp, :variable, lb = 0.01 * ones(N_vars), ub = 10 * ones(N_vars)) + dynamics!(ocp, (t, x, u, v) -> [x[2], u] * dt(t, v)) # min tf objective!(ocp, :mayer, (x0, xf, v) -> sum(v)) #= # min energy fixed tf @@ -77,24 +105,24 @@ if test5 objective!(ocp, :lagrange, (t, x, u, v) -> u^2 * v[Int(t)]) =# - sol = solve(ocp, grid_size=N_vars) + sol = solve(ocp, grid_size = N_vars) # actual time grid v = sol.variable - T_opt = zeros(N_vars+1) - for i in 1:N_vars + T_opt = zeros(N_vars + 1) + for i = 1:N_vars T_opt[i+1] = T_opt[i] + v[i] end println("Optimized time steps ", T_opt) println("And tf: ", sum(sol.variable)) - + #plot(sol) - U_opt = zeros(N_vars+1) + U_opt = zeros(N_vars + 1) # ffs julia - for i in 1:N_vars+1 + for i = 1:N_vars+1 U_opt[i] = sol.control(sol.times[i]) end - p=plot(T_opt, U_opt, markershape=:circle, show=true) + p = plot(T_opt, U_opt, markershape = :circle, show = true) end @@ -105,7 +133,7 @@ end # this feature needs a proper implementation anyway if test6 N_vars = 30 - goddard = Model(variable=true, autonomous=false) + goddard = Model(variable = true, autonomous = false) Cd = 310 Tmax = 3.5 β = 500 @@ -119,36 +147,42 @@ if test6 state!(goddard, 3) control!(goddard, 1) variable!(goddard, N_vars) - time!(goddard, t0=0, tf=N_vars) - constraint!(goddard, :initial, lb=x0, ub=x0) - constraint!(goddard, :final, rg=3, lb=mf, ub=mf) - constraint!(goddard, :state, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0]) - constraint!(goddard, :control, lb=0, ub=1) - constraint!(goddard, :variable, lb=0.05/N_vars*ones(N_vars), ub=Inf*ones(N_vars)) - objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) + time!(goddard, t0 = 0, tf = N_vars) + constraint!(goddard, :initial, lb = x0, ub = x0) + constraint!(goddard, :final, rg = 3, lb = mf, ub = mf) + constraint!(goddard, :state, lb = [r0, v0, mf], ub = [r0 + 0.2, vmax, m0]) + constraint!(goddard, :control, lb = 0, ub = 1) + constraint!( + goddard, + :variable, + lb = 0.05 / N_vars * ones(N_vars), + ub = Inf * ones(N_vars), + ) + objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) function F0(x) r, v, m = x - D = Cd * v^2 * exp(-β*(r - 1)) - return [ v, -D/m - 1/r^2, 0 ] + D = Cd * v^2 * exp(-β * (r - 1)) + return [v, -D / m - 1 / r^2, 0] end function F1(x) r, v, m = x - return [ 0, Tmax/m, -b*Tmax ] + return [0, Tmax / m, -b * Tmax] end - dynamics!(goddard, (t, x, u, v) -> (F0(x) + u*F1(x)) * dt(t,v) ) + dynamics!(goddard, (t, x, u, v) -> (F0(x) + u * F1(x)) * dt(t, v)) # start will small time steps ? - sol = solve(goddard, grid_size=N_vars, tol=1e-12, init=(variable=zeros(N_vars),)) + sol = + solve(goddard, grid_size = N_vars, tol = 1e-12, init = (variable = zeros(N_vars),)) # actual time grid v = sol.variable - T_opt = zeros(N_vars+1) - for i in 1:N_vars + T_opt = zeros(N_vars + 1) + for i = 1:N_vars T_opt[i+1] = T_opt[i] + v[i] end println("Optimized time steps ", T_opt) println("And tf: ", sum(sol.variable)) - + U_opt = sol.control.(sol.times) - p=plot(T_opt, U_opt, markershape=:circle, show=true) -end \ No newline at end of file + p = plot(T_opt, U_opt, markershape = :circle, show = true) +end diff --git a/test/runtests.jl b/test/runtests.jl index d563dcd2..e8c397c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,5 +5,5 @@ include("deps.jl") # +++ add explicit list (cf OC) (keep suite folder) @testset verbose = true showtiming = true "Test suite" begin # run all scripts in subfolder suite/ - include.(filter(contains(r".jl$"), readdir("./suite"; join=true))) + include.(filter(contains(r".jl$"), readdir("./suite"; join = true))) end diff --git a/test/scalar_vector.jl b/test/scalar_vector.jl index c73c7c11..5bf9abb7 100644 --- a/test/scalar_vector.jl +++ b/test/scalar_vector.jl @@ -2,7 +2,7 @@ n = 3 m = 1 N = 10 -xu = ones((N+1)*n+N*m) +xu = ones((N + 1) * n + N * m) function getx(xu, step, n) if n == 1 @@ -16,7 +16,7 @@ function getu(xu, step, n, m) #if m == 1 # return xu[(N+1)*n + step*m + 1] #else - return xu[(N+1)*n + step*m + 1 : (N+1)*n + step*m + m] + return xu[(N+1)*n+step*m+1:(N+1)*n+step*m+m] #end end @@ -24,7 +24,7 @@ function dynamics_scalar(x, u, n) f = zeros(n) f[1] = x[2] f[2] = u - f[3] = u*u + f[3] = u * u return f end @@ -32,13 +32,13 @@ function dynamics_vector(x, u, n) f = zeros(n) f[1] = x[2] f[2] = u[1] - f[3] = u[1]*u[1] + f[3] = u[1] * u[1] return f end # ici u semble bien de type scalaire des la premiere evaluation -@code_warntype dynamics_scalar(getx(xu,5,3), getu(xu,5,3,1)) -@code_warntype dynamics_vector(getx(xu,5,3), getu(xu,5,3,1)) +@code_warntype dynamics_scalar(getx(xu, 5, 3), getu(xu, 5, 3, 1)) +@code_warntype dynamics_vector(getx(xu, 5, 3), getu(xu, 5, 3, 1)) # nb le code dynamics_scalar passe aussi avec un u vectoriel -# au niveau de CTBase qu'est ce qui coince exactement si on attend un argument scalaire et qu'on recupere un vecteur de taille 1 ? \ No newline at end of file +# au niveau de CTBase qu'est ce qui coince exactement si on attend un argument scalaire et qu'on recupere un vecteur de taille 1 ? diff --git a/test/suite/test_abstract_ocp.jl b/test/suite/test_abstract_ocp.jl index f106f18d..a88d1505 100644 --- a/test/suite/test_abstract_ocp.jl +++ b/test/suite/test_abstract_ocp.jl @@ -7,8 +7,8 @@ end @testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin prob = double_integrator_a() - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end if !isdefined(Main, :goddard_a) @@ -16,6 +16,6 @@ if !isdefined(Main, :goddard_a) end @testset verbose = true showtiming = true ":goddard :max_rf :abstract :constr" begin prob = goddard_a() - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end diff --git a/test/suite/test_constraints.jl b/test/suite/test_constraints.jl index 1c0287be..405ff9e5 100644 --- a/test/suite/test_constraints.jl +++ b/test/suite/test_constraints.jl @@ -15,7 +15,7 @@ if check_constraint_mult pyplot() ocp = goddard_all() - sol = direct_solve(ocp.ocp, display=false, init=ocp.init) + sol = direct_solve(ocp.ocp, display = false, init = ocp.init) # plot state, control and costate psol = plot(sol) @@ -31,7 +31,12 @@ if check_constraint_mult # variable constraints / box println("\nVariable constraints: ", sol.variable_constraints) println("multipliers: ", sol.mult_variable_constraints) - println("\nVariable box multipliers LB: ", sol.mult_variable_box_lower, " UB ", sol.mult_variable_box_upper) + println( + "\nVariable box multipliers LB: ", + sol.mult_variable_box_lower, + " UB ", + sol.mult_variable_box_upper, + ) # PATH CONSTRAINTS # state box @@ -39,27 +44,34 @@ if check_constraint_mult x_box_lb = invert(sol.mult_state_box_lower.(T)) x_box_ub = invert(sol.mult_state_box_upper.(T)) - p1_a = plot(T, x[1], label="r") - p1_b = plot(T, [x_box_lb[1] x_box_ub[1]], label=["LB" "UB"]) - p1 = plot(p1_a, p1_b, layout=(2,1)) + p1_a = plot(T, x[1], label = "r") + p1_b = plot(T, [x_box_lb[1] x_box_ub[1]], label = ["LB" "UB"]) + p1 = plot(p1_a, p1_b, layout = (2, 1)) - p2_a = plot(T, x[2], label="v") - p2_b = plot(T, [x_box_lb[2] x_box_ub[2]], label=["LB" "UB"]) - p2 = plot(p2_a, p2_b, layout=(2,1)) + p2_a = plot(T, x[2], label = "v") + p2_b = plot(T, [x_box_lb[2] x_box_ub[2]], label = ["LB" "UB"]) + p2 = plot(p2_a, p2_b, layout = (2, 1)) - p3_a = plot(T, x[3], label="m") - p3_b = plot(T, [x_box_lb[3] x_box_ub[3]], label=["LB" "UB"]) - p3 = plot(p3_a, p3_b, layout=(2,1)) + p3_a = plot(T, x[3], label = "m") + p3_b = plot(T, [x_box_lb[3] x_box_ub[3]], label = ["LB" "UB"]) + p3 = plot(p3_a, p3_b, layout = (2, 1)) # control box u = sol.control.(T) u_box_lb = flatten(sol.mult_control_box_lower.(T)) u_box_ub = flatten(sol.mult_control_box_upper.(T)) - p4_a = plot(T, u, label="u") - p4_b = plot(T, [u_box_lb u_box_ub], label=["LB" "UB"]) - p4 = plot(p4_a, p4_b, layout=(2,1)) - - p_box = plot(p1, p2, p3, p4, layout=(2,2), title=["r box" "" "v box" "" "m box" "" "u box" ""]) + p4_a = plot(T, u, label = "u") + p4_b = plot(T, [u_box_lb u_box_ub], label = ["LB" "UB"]) + p4 = plot(p4_a, p4_b, layout = (2, 1)) + + p_box = plot( + p1, + p2, + p3, + p4, + layout = (2, 2), + title = ["r box" "" "v box" "" "m box" "" "u box" ""], + ) display(p_box) readline() #ffs julia, fix your damn plots @@ -67,47 +79,53 @@ if check_constraint_mult # control constraints c_u = flatten(sol.control_constraints.(T)) m_c_u = flatten(sol.mult_control_constraints.(T)) - p5_a = plot(T, c_u, label="c_u") - p5_b = plot(T, m_c_u, label="mul") - p5 = plot(p5_a, p5_b, layout=(2,1)) + p5_a = plot(T, c_u, label = "c_u") + p5_b = plot(T, m_c_u, label = "mul") + p5 = plot(p5_a, p5_b, layout = (2, 1)) # state constraints c_x = flatten(sol.state_constraints.(T)) m_c_x = flatten(sol.mult_state_constraints.(T)) - p6_a = plot(T, c_x, label="c_x") - p6_b = plot(T, m_c_x, label="mul") - p6 = plot(p6_a, p6_b, layout=(2,1)) + p6_a = plot(T, c_x, label = "c_x") + p6_b = plot(T, m_c_x, label = "mul") + p6 = plot(p6_a, p6_b, layout = (2, 1)) # mixed constraints c_xu = flatten(sol.mixed_constraints.(T)) m_c_xu = flatten(sol.mult_mixed_constraints.(T)) - p7_a = plot(T, c_xu, label="c_xu") - p7_b = plot(T, m_c_xu, label="mul") - p7 = plot(p7_a, p7_b, layout=(2,1)) - - p_cons = plot(p5, p6, p7, layout=(1,3), title=["control cons" "" "state cons" "" "mixed cons" ""]) + p7_a = plot(T, c_xu, label = "c_xu") + p7_b = plot(T, m_c_xu, label = "mul") + p7 = plot(p7_a, p7_b, layout = (2, 1)) + + p_cons = plot( + p5, + p6, + p7, + layout = (1, 3), + title = ["control cons" "" "state cons" "" "mixed cons" ""], + ) else -# box constraints -@testset verbose = true showtiming = true ":goddard :box_constraints" begin - ocp = goddard() - sol = direct_solve(ocp.ocp, display=false) - @test sol.objective ≈ ocp.obj rtol=1e-2 -end + # box constraints + @testset verbose = true showtiming = true ":goddard :box_constraints" begin + ocp = goddard() + sol = direct_solve(ocp.ocp, display = false) + @test sol.objective ≈ ocp.obj rtol = 1e-2 + end + + # functional constraints + @testset verbose = true showtiming = true ":goddard :functional_constraints" begin + ocp = goddard(functional_constraints = true) + sol = direct_solve(ocp.ocp, display = false, init = ocp.init) + @test sol.objective ≈ ocp.obj rtol = 1e-2 + end + + # all constraints + @testset verbose = true showtiming = true ":goddard :all_constraints" begin + ocp = goddard_all() + sol = direct_solve(ocp.ocp, display = false, init = ocp.init) + @test sol.objective ≈ ocp.obj rtol = 1e-2 + end -# functional constraints -@testset verbose = true showtiming = true ":goddard :functional_constraints" begin - ocp = goddard(functional_constraints=true) - sol = direct_solve(ocp.ocp, display=false, init=ocp.init) - @test sol.objective ≈ ocp.obj rtol=1e-2 end - -# all constraints -@testset verbose = true showtiming = true ":goddard :all_constraints" begin - ocp = goddard_all() - sol = direct_solve(ocp.ocp, display=false, init=ocp.init) - @test sol.objective ≈ ocp.obj rtol=1e-2 -end - -end \ No newline at end of file diff --git a/test/suite/test_continuation.jl b/test/suite/test_continuation.jl index c79c7522..a3d89175 100644 --- a/test/suite/test_continuation.jl +++ b/test/suite/test_continuation.jl @@ -13,13 +13,13 @@ if test1 @testset verbose = true showtiming = true ":continuation :double_integrator" begin init = () obj_list = [] - for T=1:5 - ocp = double_integrator_T(T).ocp - sol = direct_solve(ocp, display=false, init=init) + for T = 1:5 + ocp = double_integrator_T(T).ocp + sol = direct_solve(ocp, display = false, init = init) init = sol push!(obj_list, sol.objective) end - @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol=1e-2 + @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol = 1e-2 end end @@ -35,11 +35,11 @@ if test2 obj_list = [] for ρ in [0.1, 5, 10, 30, 100] ocp = parametric(ρ).ocp - sol = direct_solve(ocp, display=false, init=init) + sol = direct_solve(ocp, display = false, init = init) init = sol push!(obj_list, sol.objective) end - @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol=1e-2 + @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol = 1e-2 end end @@ -49,27 +49,34 @@ if test3 if !isdefined(Main, :goddard) include("../problems/goddard.jl") end - sol0 = direct_solve(goddard().ocp, display=false) + sol0 = direct_solve(goddard().ocp, display = false) @testset verbose = true showtiming = true ":continuation :goddard" begin sol = sol0 Tmax_list = [] obj_list = [] - for Tmax=3.5:-0.5:1 - sol = direct_solve(goddard(Tmax=Tmax).ocp, display=false, init=sol) + for Tmax = 3.5:-0.5:1 + sol = direct_solve(goddard(Tmax = Tmax).ocp, display = false, init = sol) push!(Tmax_list, Tmax) push!(obj_list, sol.objective) end - @test obj_list ≈ [1.0125, 1.0124, 1.0120, 1.0112, 1.0092, 1.0036] rtol=1e-2 + @test obj_list ≈ [1.0125, 1.0124, 1.0120, 1.0112, 1.0092, 1.0036] rtol = 1e-2 if draw_plot using Plots # plot obj(vmax) - pobj = plot(Tmax_list, obj_list, label="r(tf)", xlabel="Maximal thrust (Tmax)", ylabel="Maximal altitude r(tf)",seriestype=:scatter) + pobj = plot( + Tmax_list, + obj_list, + label = "r(tf)", + xlabel = "Maximal thrust (Tmax)", + ylabel = "Maximal altitude r(tf)", + seriestype = :scatter, + ) # plot multiple solutions plot(sol0) p = plot!(sol) - display(plot(pobj, p, layout=2, reuse=false, size=(1000,500))) + display(plot(pobj, p, layout = 2, reuse = false, size = (1000, 500))) end end end diff --git a/test/suite/test_grid.jl b/test/suite/test_grid.jl index d2edb53e..09591a83 100644 --- a/test/suite/test_grid.jl +++ b/test/suite/test_grid.jl @@ -5,19 +5,19 @@ if !isdefined(Main, :simple_integrator) include("../problems/simple_integrator.jl") end ocp = simple_integrator().ocp -sol0 = direct_solve(ocp, display=false) +sol0 = direct_solve(ocp, display = false) # solve with explicit and non uniform time grid @testset verbose = true showtiming = true ":explicit_grid" begin - time_grid = LinRange(0,1,CTDirect.__grid_size()+1) - sol = direct_solve(ocp, time_grid=time_grid, display=false) + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1) + sol = direct_solve(ocp, time_grid = time_grid, display = false) @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) end @testset verbose = true showtiming = true ":non_uniform_grid" begin - time_grid = [0,0.1,0.3,0.6,0.98,0.99,1] - sol = direct_solve(ocp, time_grid=time_grid, display=false) - @test sol.objective ≈ 0.309 rtol=1e-2 + time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1] + sol = direct_solve(ocp, time_grid = time_grid, display = false) + @test sol.objective ≈ 0.309 rtol = 1e-2 end @@ -26,16 +26,20 @@ if !isdefined(Main, :double_integrator_freet0tf) include("../problems/double_integrator.jl") end ocp = double_integrator_freet0tf().ocp -sol0 = direct_solve(ocp, display=false) +sol0 = direct_solve(ocp, display = false) @testset verbose = true showtiming = true ":explicit_grid" begin - sol = direct_solve(ocp, time_grid=LinRange(0,1,CTDirect.__grid_size()+1), display=false) + sol = direct_solve( + ocp, + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), + display = false, + ) @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) end @testset verbose = true showtiming = true ":max_t0 :non_uniform_grid" begin - sol = direct_solve(ocp, time_grid=[0,0.1,0.6,0.95,1], display=false) - @test sol.objective ≈ 7.48 rtol=1e-2 + sol = direct_solve(ocp, time_grid = [0, 0.1, 0.6, 0.95, 1], display = false) + @test sol.objective ≈ 7.48 rtol = 1e-2 end @@ -44,14 +48,18 @@ if !isdefined(Main, :double_integrator_T) include("../problems/double_integrator.jl") end ocp = double_integrator_T(2).ocp -sol0 = direct_solve(ocp, display=false) +sol0 = direct_solve(ocp, display = false) @testset verbose = true showtiming = true ":explicit_grid" begin - sol = direct_solve(ocp, time_grid=LinRange(0,1,CTDirect.__grid_size()+1), display=false) + sol = direct_solve( + ocp, + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), + display = false, + ) @test (sol.objective == sol0.objective) && (sol.iterations == sol0.iterations) end @testset verbose = true showtiming = true ":non_uniform_grid" begin - sol = direct_solve(ocp, time_grid=[0,0.3,1,1.9,2], display=false) - @test sol.objective ≈ 2.43 rtol=1e-2 + sol = direct_solve(ocp, time_grid = [0, 0.3, 1, 1.9, 2], display = false) + @test sol.objective ≈ 2.43 rtol = 1e-2 end diff --git a/test/suite/test_initial_guess.jl b/test/suite/test_initial_guess.jl index 0e36d2cf..e25fe135 100644 --- a/test/suite/test_initial_guess.jl +++ b/test/suite/test_initial_guess.jl @@ -17,7 +17,7 @@ end # reference solution prob = double_integrator_a() ocp = prob.ocp -sol0 = direct_solve(ocp, display=false) +sol0 = direct_solve(ocp, display = false) # constant initial guess x_const = [0.5, 0.2] @@ -25,13 +25,13 @@ u_const = 0.5 v_const = 0.15 # functional initial guess -x_func = t->[t^2, sqrt(t)] -u_func = t->(cos(10*t)+1)*0.5 +x_func = t -> [t^2, sqrt(t)] +u_func = t -> (cos(10 * t) + 1) * 0.5 # interpolated initial guess x_vec = [[0, 0], [1, 2], [5, -1]] x_matrix = [0 0; 1 2; 5 -1] -u_vec = [0, 0.3, .1] +u_vec = [0, 0.3, 0.1] ################################################# @@ -39,88 +39,146 @@ u_vec = [0, 0.3, .1] # 1.a default initial guess @testset verbose = true showtiming = true ":default_init_no_arg" begin - sol = direct_solve(ocp, display=false, max_iter=maxiter) + sol = direct_solve(ocp, display = false, max_iter = maxiter) @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) end @testset verbose = true showtiming = true ":default_init_()" begin - sol = direct_solve(ocp, display=false, init=(), max_iter=maxiter) + sol = direct_solve(ocp, display = false, init = (), max_iter = maxiter) @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) end @testset verbose = true showtiming = true ":default_init_nothing" begin - sol = direct_solve(ocp, display=false, init=nothing,max_iter=maxiter) + sol = direct_solve(ocp, display = false, init = nothing, max_iter = maxiter) @test(check_xf(sol, [0.1, 0.1]) && check_uf(sol, 0.1) && check_v(sol, 0.1)) end # 1.b constant initial guess @testset verbose = true showtiming = true ":constant_x" begin - sol = direct_solve(ocp, display=false, init=(state=x_const,), max_iter=maxiter) + sol = direct_solve(ocp, display = false, init = (state = x_const,), max_iter = maxiter) @test(check_xf(sol, x_const)) end @testset verbose = true showtiming = true ":constant_u" begin - sol = direct_solve(ocp, display=false, init=(control=u_const,), max_iter=maxiter) + sol = + direct_solve(ocp, display = false, init = (control = u_const,), max_iter = maxiter) @test(check_uf(sol, u_const)) end @testset verbose = true showtiming = true ":constant_v" begin - sol = direct_solve(ocp, display=false, init=(variable=v_const,), max_iter=maxiter) + sol = + direct_solve(ocp, display = false, init = (variable = v_const,), max_iter = maxiter) @test(check_v(sol, v_const)) end @testset verbose = true showtiming = true ":constant_xu" begin - sol = direct_solve(ocp, display=false, init=(state=x_const, control=u_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (state = x_const, control = u_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_const) && check_uf(sol, u_const)) end @testset verbose = true showtiming = true ":constant_xv" begin - sol = direct_solve(ocp, display=false, init=(state=x_const, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (state = x_const, variable = v_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_const) && check_v(sol, v_const)) end @testset verbose = true showtiming = true ":constant_uv" begin - sol = direct_solve(ocp, display=false, init=(control=u_const, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (control = u_const, variable = v_const), + max_iter = maxiter, + ) @test(check_uf(sol, u_const) && check_v(sol, v_const)) end @testset verbose = true showtiming = true ":constant_xuv" begin - sol = direct_solve(ocp, display=false, init=(state=x_const, control=u_const, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (state = x_const, control = u_const, variable = v_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_const) && check_uf(sol, u_const) && check_v(sol, v_const)) -end +end # 1. functional initial guess @testset verbose = true showtiming = true ":functional_x" begin - sol = direct_solve(ocp, display=false, init=(state=x_func,), max_iter=maxiter) + sol = direct_solve(ocp, display = false, init = (state = x_func,), max_iter = maxiter) @test(check_xf(sol, x_func(sol.time_grid[end]))) end @testset verbose = true showtiming = true ":functional_u" begin - sol = direct_solve(ocp, display=false, init=(control=u_func,), max_iter=maxiter) + sol = direct_solve(ocp, display = false, init = (control = u_func,), max_iter = maxiter) @test(check_uf(sol, u_func(sol.time_grid[end]))) end @testset verbose = true showtiming = true ":functional_xu" begin - sol = direct_solve(ocp, display=false, init=(state=x_func, control=u_func), max_iter=maxiter) - @test(check_xf(sol, x_func(sol.time_grid[end])) && check_uf(sol, u_func(sol.time_grid[end]))) + sol = direct_solve( + ocp, + display = false, + init = (state = x_func, control = u_func), + max_iter = maxiter, + ) + @test( + check_xf(sol, x_func(sol.time_grid[end])) && + check_uf(sol, u_func(sol.time_grid[end])) + ) end # 1.d interpolated initial guess -t_vec = [0, .1, v_const] +t_vec = [0, 0.1, v_const] @testset verbose = true showtiming = true ":vector_txu :constant_v" begin - sol = direct_solve(ocp, display=false, init=(time=t_vec, state=x_vec, control=u_vec, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (time = t_vec, state = x_vec, control = u_vec, variable = v_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) end -t_matrix = [0 .1 v_const] +t_matrix = [0 0.1 v_const] @testset verbose = true showtiming = true ":matrix_t :vector_xu :constant_v" begin - sol = direct_solve(ocp, display=false, init=(time=t_matrix, state=x_vec, control=u_vec, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (time = t_matrix, state = x_vec, control = u_vec, variable = v_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) end @testset verbose = true showtiming = true ":matrix_x :vector_tu :constant_v" begin - sol = direct_solve(ocp, display=false, init=(time=t_vec, state=x_matrix, control=u_vec, variable=v_const), max_iter=maxiter) + sol = direct_solve( + ocp, + display = false, + init = (time = t_vec, state = x_matrix, control = u_vec, variable = v_const), + max_iter = maxiter, + ) @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_vec[end]) && check_v(sol, v_const)) end # 1.e mixed initial guess @testset verbose = true showtiming = true ":vector_tx :functional_u :constant_v" begin - sol = direct_solve(ocp, display=false, init=(time=t_vec, state=x_vec, control=u_func, variable=v_const), max_iter=maxiter) - @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.time_grid[end])) && check_v(sol, v_const)) + sol = direct_solve( + ocp, + display = false, + init = (time = t_vec, state = x_vec, control = u_func, variable = v_const), + max_iter = maxiter, + ) + @test( + check_xf(sol, x_vec[end]) && + check_uf(sol, u_func(sol.time_grid[end])) && + check_v(sol, v_const) + ) end # 1.f warm start @testset verbose = true showtiming = true ":warm_start" begin - sol = direct_solve(ocp, display=false, init=sol0, max_iter=maxiter) - @test(check_xf(sol, sol.state(sol.time_grid[end])) && check_uf(sol, sol.control(sol.time_grid[end])) && check_v(sol, sol.variable)) + sol = direct_solve(ocp, display = false, init = sol0, max_iter = maxiter) + @test( + check_xf(sol, sol.state(sol.time_grid[end])) && + check_uf(sol, sol.control(sol.time_grid[end])) && + check_v(sol, sol.variable) + ) end ################################################# @@ -129,15 +187,27 @@ docp, nlp = direct_transcription(ocp) tag = CTDirect.IpoptTag() # mixed init @testset verbose = true showtiming = true ":docp_mixed_init" begin - set_initial_guess(docp, nlp, (time=t_vec, state=x_vec, control=u_func, variable=v_const)) - dsol = CTDirect.solve_docp(tag, docp, nlp, display=false, max_iter=maxiter) + set_initial_guess( + docp, + nlp, + (time = t_vec, state = x_vec, control = u_func, variable = v_const), + ) + dsol = CTDirect.solve_docp(tag, docp, nlp, display = false, max_iter = maxiter) sol = OptimalControlSolution(docp, dsol) - @test(check_xf(sol, x_vec[end]) && check_uf(sol, u_func(sol.time_grid[end])) && check_v(sol, v_const)) + @test( + check_xf(sol, x_vec[end]) && + check_uf(sol, u_func(sol.time_grid[end])) && + check_v(sol, v_const) + ) end # warm start @testset verbose = true showtiming = true ":docp_warm_start" begin set_initial_guess(docp, nlp, sol0) - dsol = CTDirect.solve_docp(tag, docp, nlp, display=false, max_iter=maxiter) + dsol = CTDirect.solve_docp(tag, docp, nlp, display = false, max_iter = maxiter) sol = OptimalControlSolution(docp, dsol) - @test(check_xf(sol, sol.state(sol.time_grid[end])) && check_uf(sol, sol.control(sol.time_grid[end])) && check_v(sol, sol.variable)) + @test( + check_xf(sol, sol.state(sol.time_grid[end])) && + check_uf(sol, sol.control(sol.time_grid[end])) && + check_v(sol, sol.variable) + ) end diff --git a/test/suite/test_misc.jl b/test/suite/test_misc.jl index 582834a3..b430b639 100644 --- a/test/suite/test_misc.jl +++ b/test/suite/test_misc.jl @@ -20,11 +20,11 @@ if !isdefined(Main, :simple_integrator) include("../problems/simple_integrator.jl") end ocp = simple_integrator().ocp -sol0 = direct_solve(ocp, display=false) +sol0 = direct_solve(ocp, display = false) # test save / load solution in JLD2 format @testset verbose = true showtiming = true ":save_load :JLD2" begin - save(sol0, filename_prefix="solution_test") + save(sol0, filename_prefix = "solution_test") sol_reloaded = load("solution_test") @test sol0.objective == sol_reloaded.objective end @@ -32,8 +32,8 @@ end #= # test export / read solution in JSON format @testset verbose = true showtiming = true ":export_read :JSON" begin - export_ocp_solution(sol0, filename_prefix="solution_test") + export_ocp_solution(sol0, filename_prefix = "solution_test") sol_reloaded = import_ocp_solution("solution_test") @test sol0.objective == sol_reloaded.objective end -=# \ No newline at end of file +=# diff --git a/test/suite/test_nlp.jl b/test/suite/test_nlp.jl index 4c4382b5..884aa806 100644 --- a/test/suite/test_nlp.jl +++ b/test/suite/test_nlp.jl @@ -7,32 +7,32 @@ ocp = simple_integrator().ocp @testset verbose = true showtiming = true ":control_dim_2" begin @test is_solvable(ocp) - @test (:adnlp,:ipopt) in available_methods() - @test (:adnlp,:madnlp) in available_methods() + @test (:adnlp, :ipopt) in available_methods() + @test (:adnlp, :madnlp) in available_methods() end @testset verbose = true showtiming = true ":solve_docp :ipopt" begin docp, nlp = direct_transcription(ocp) tag = CTDirect.IpoptTag() - dsol = CTDirect.solve_docp(tag, docp, nlp, display=false) + dsol = CTDirect.solve_docp(tag, docp, nlp, display = false) sol = OptimalControlSolution(docp, dsol) - @test sol.objective ≈ 0.313 rtol=1e-2 - sol = OptimalControlSolution(docp, primal=dsol.solution) - @test sol.objective ≈ 0.313 rtol=1e-2 - sol = OptimalControlSolution(docp, primal=dsol.solution, dual=dsol.multipliers) - @test sol.objective ≈ 0.313 rtol=1e-2 + @test sol.objective ≈ 0.313 rtol = 1e-2 + sol = OptimalControlSolution(docp, primal = dsol.solution) + @test sol.objective ≈ 0.313 rtol = 1e-2 + sol = OptimalControlSolution(docp, primal = dsol.solution, dual = dsol.multipliers) + @test sol.objective ≈ 0.313 rtol = 1e-2 end @testset verbose = true showtiming = true ":solve_docp :madnlp" begin docp, nlp = direct_transcription(ocp) tag = CTDirect.MadNLPTag() - dsol = CTDirect.solve_docp(tag, docp, nlp, display=false) + dsol = CTDirect.solve_docp(tag, docp, nlp, display = false) sol = OptimalControlSolution(docp, dsol) - @test sol.objective ≈ 0.313 rtol=1e-2 - sol = OptimalControlSolution(docp, primal=dsol.solution) - @test sol.objective ≈ 0.313 rtol=1e-2 - sol = OptimalControlSolution(docp, primal=dsol.solution, dual=dsol.multipliers) - @test sol.objective ≈ 0.313 rtol=1e-2 + @test sol.objective ≈ 0.313 rtol = 1e-2 + sol = OptimalControlSolution(docp, primal = dsol.solution) + @test sol.objective ≈ 0.313 rtol = 1e-2 + sol = OptimalControlSolution(docp, primal = dsol.solution, dual = dsol.multipliers) + @test sol.objective ≈ 0.313 rtol = 1e-2 end # check solution building @@ -40,22 +40,22 @@ if !isdefined(Main, :double_integrator_T) include("../problems/double_integrator.jl") end ocp = double_integrator_T(1).ocp -x_opt = t -> [6*(t^2 / 2 - t^3 / 3), 6*(t - t^2)] -u_opt = t -> 6 - 12*t -p_opt = t -> [24 , 12 - 24*t] +x_opt = t -> [6 * (t^2 / 2 - t^3 / 3), 6 * (t - t^2)] +u_opt = t -> 6 - 12 * t +p_opt = t -> [24, 12 - 24 * t] @testset verbose = true showtiming = true ":analytic_solution :ipopt" begin - sol = direct_solve(ocp, display=false) + sol = direct_solve(ocp, display = false) T = sol.time_grid - @test isapprox(x_opt.(T), sol.state.(T), rtol=1e-2) - @test isapprox(u_opt.(T), sol.control.(T), rtol=1e-2) - @test isapprox(p_opt.(T), sol.costate.(T), rtol=1e-2) + @test isapprox(x_opt.(T), sol.state.(T), rtol = 1e-2) + @test isapprox(u_opt.(T), sol.control.(T), rtol = 1e-2) + @test isapprox(p_opt.(T), sol.costate.(T), rtol = 1e-2) end @testset verbose = true showtiming = true ":analytic_solution :madnlp" begin - sol = direct_solve(ocp, :madnlp, display=false) + sol = direct_solve(ocp, :madnlp, display = false) T = sol.time_grid - @test isapprox(x_opt.(T), sol.state.(T), rtol=1e-2) - @test isapprox(u_opt.(T), sol.control.(T), rtol=1e-2) - @test isapprox(p_opt.(T), sol.costate.(T), rtol=1e-2) + @test isapprox(x_opt.(T), sol.state.(T), rtol = 1e-2) + @test isapprox(u_opt.(T), sol.control.(T), rtol = 1e-2) + @test isapprox(p_opt.(T), sol.costate.(T), rtol = 1e-2) end diff --git a/test/suite/test_objective.jl b/test/suite/test_objective.jl index 189b7057..72a4874a 100644 --- a/test/suite/test_objective.jl +++ b/test/suite/test_objective.jl @@ -7,22 +7,22 @@ end @testset verbose = true showtiming = true ":min_tf :mayer" begin prob = double_integrator_mintf() - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end # min tf (lagrange) @testset verbose = true showtiming = true ":min_tf :lagrange" begin - prob = double_integrator_mintf(lagrange=true) - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + prob = double_integrator_mintf(lagrange = true) + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end # max t0 (free t0 and tf) @testset verbose = true showtiming = true ":max_t0" begin prob = double_integrator_freet0tf() - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end # bolza, non-autonomous mayer term, tf in dynamics @@ -31,8 +31,8 @@ if !isdefined(Main, :bolza_freetf) end prob = bolza_freetf() @testset verbose = true showtiming = true ":bolza :tf_in_dyn_and_cost" begin - sol = direct_solve(prob.ocp, display=false) - @test sol.objective ≈ prob.obj rtol=1e-2 + sol = direct_solve(prob.ocp, display = false) + @test sol.objective ≈ prob.obj rtol = 1e-2 end #= @@ -65,4 +65,4 @@ end ∫(u(s)^2) → min end sol = direct_solve(ocp2) -=# \ No newline at end of file +=#