From 75dd4f40ca92a6f57bba2699608e6fb8ab3643d9 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Sun, 25 Aug 2024 11:19:15 +0200 Subject: [PATCH] formatter --- .github/workflows/Formatter.yml | 33 ++++ docs/make.jl | 46 ++--- src/OptimalControl.jl | 12 +- src/solve.jl | 11 +- test/Goddard.jl | 26 +-- test/problems/beam.jl | 12 +- test/problems/bioreactor.jl | 54 +++--- test/problems/fuller.jl | 12 +- test/problems/goddard.jl | 62 ++++--- test/problems/insurance.jl | 26 ++- test/problems/jackson.jl | 18 +- test/problems/parametric.jl | 16 +- test/problems/robbins.jl | 14 +- test/problems/swimmer.jl | 129 +++++++++++--- test/problems/vanderpol.jl | 14 +- test/runtests.jl | 6 +- test/test-nlp.jl | 22 ++- test/test_abstract_ocp.jl | 46 ++--- test/test_basic.jl | 16 +- test/test_continuation.jl | 181 ++++++++++--------- test/test_goddard_direct.jl | 14 +- test/test_goddard_indirect.jl | 42 ++--- test/test_grid.jl | 104 ++++++----- test/test_initial_guess.jl | 301 +++++++++++++++++++------------- test/test_objective.jl | 116 ++++++------ 25 files changed, 776 insertions(+), 557 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/docs/make.jl b/docs/make.jl index 36281906..82cf9043 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,8 +8,8 @@ using CTDirect # to add docstrings from external packages Modules = [CTBase, CTFlows, CTDirect, OptimalControl] for Module ∈ Modules - isnothing( DocMeta.getdocmeta(Module, :DocTestSetup) ) && - DocMeta.setdocmeta!(Module, :DocTestSetup, :(using $Module); recursive = true) + isnothing(DocMeta.getdocmeta(Module, :DocTestSetup)) && + DocMeta.setdocmeta!(Module, :DocTestSetup, :(using $Module); recursive = true) end # @@ -18,10 +18,7 @@ makedocs( sitename = "OptimalControl.jl", format = Documenter.HTML( prettyurls = false, - size_threshold_ignore = [ - "api-ctbase/types.md", - "tutorial-plot.md", - ], + size_threshold_ignore = ["api-ctbase/types.md", "tutorial-plot.md"], assets = [ asset("https://control-toolbox.org/assets/css/documentation.css"), asset("https://control-toolbox.org/assets/js/documentation.js"), @@ -30,47 +27,40 @@ makedocs( pages = [ "Introduction" => "index.md", "Basic examples" => [ - "Energy minimisation" => "tutorial-basic-example.md", + "Energy minimisation" => "tutorial-basic-example.md", "Time mininimisation" => "tutorial-double-integrator.md", ], "Manual" => [ - "Abstract syntax" => "tutorial-abstract.md", - "Initial guess" => "tutorial-initial-guess.md", - "Solve" => "tutorial-solve.md", - "Plot a solution" => "tutorial-plot.md", - "Flow" => "tutorial-flow.md", - "Functional syntax" => "tutorial-functional.md", - "Control-toolbox REPL" => "tutorial-repl.md", + "Abstract syntax" => "tutorial-abstract.md", + "Initial guess" => "tutorial-initial-guess.md", + "Solve" => "tutorial-solve.md", + "Plot a solution" => "tutorial-plot.md", + "Flow" => "tutorial-flow.md", + "Functional syntax" => "tutorial-functional.md", + "Control-toolbox REPL" => "tutorial-repl.md", ], "Tutorials" => [ "tutorial-continuation.md", - "Goddard: direct, indirect" => "tutorial-goddard.md", + "Goddard: direct, indirect" => "tutorial-goddard.md", "tutorial-iss.md", "Linear–quadratic regulator" => "tutorial-lqr-basic.md", "tutorial-nlp.md", ], "API" => [ "api-optimalcontrol.md", - "Subpackages" => [ - "api-ctbase.md", - "api-ctdirect.md", - "api-ctflows.md", - ], + "Subpackages" => ["api-ctbase.md", "api-ctdirect.md", "api-ctflows.md"], ], "Developers" => [ "OptimalControl.jl" => "dev-optimalcontrol.md", "Subpackages" => [ - "CTBase.jl" => "dev-ctbase.md", + "CTBase.jl" => "dev-ctbase.md", "CTDirect.jl" => "dev-ctdirect.md", - "CTFlows.jl" => "dev-ctflows.md", + "CTFlows.jl" => "dev-ctflows.md", ], ], - "JuliaCon 2024"=> "juliacon2024.md", - ] + "JuliaCon 2024" => "juliacon2024.md", + ], ) # -deploydocs( - repo = "github.com/control-toolbox/OptimalControl.jl.git", - devbranch = "main" -) +deploydocs(repo = "github.com/control-toolbox/OptimalControl.jl.git", devbranch = "main") diff --git a/src/OptimalControl.jl b/src/OptimalControl.jl index fdcb0778..35c7a731 100644 --- a/src/OptimalControl.jl +++ b/src/OptimalControl.jl @@ -23,7 +23,7 @@ using CTFlows import CommonSolve: solve, CommonSolve # declarations -const __display = CTBase.__display +const __display = CTBase.__display const __ocp_init = CTBase.__ocp_init include("solve.jl") @@ -53,7 +53,15 @@ export OptimalControlModel, OptimalControlSolution export Autonomous, NonAutonomous export NonFixed, Fixed export Model, __OCPModel -export variable!, time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!, constraint +export variable!, + time!, + constraint!, + dynamics!, + objective!, + state!, + control!, + remove_constraint!, + constraint #export is_time_independent, is_time_dependent, is_min, is_max, is_variable_dependent, is_variable_independent export Lie, @Lie, Poisson, Lift, ⋅, ∂ₜ export @def diff --git a/src/solve.jl b/src/solve.jl index 29056cff..13bcdd05 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -19,7 +19,7 @@ Remove from the description, the Symbol that are specific to [OptimalControl.jl] be passed. """ function clean(d::Description) - return remove(d, (:direct, )) + return remove(d, (:direct,)) end """ @@ -52,15 +52,12 @@ julia> sol = solve(ocp, init=(state=t->[-1+t, t*(t-1)], control=0.5)) julia> sol = solve(ocp, init=(state=t->[-1+t, t*(t-1)], control=t->6-12*t)) ``` """ -function CommonSolve.solve( - ocp::OptimalControlModel, - description::Symbol...; - kwargs...) +function CommonSolve.solve(ocp::OptimalControlModel, description::Symbol...; kwargs...) # get the full description method = getFullDescription(description, available_methods()) # solve the problem - :direct ∈ method && return CTDirect.direct_solve(ocp, clean(description)...; kwargs...) + :direct ∈ method && return CTDirect.direct_solve(ocp, clean(description)...; kwargs...) -end \ No newline at end of file +end diff --git a/test/Goddard.jl b/test/Goddard.jl index b9f1d50f..42c0746c 100644 --- a/test/Goddard.jl +++ b/test/Goddard.jl @@ -11,7 +11,7 @@ function Goddard() vmax = 0.1 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] # the model @def ocp begin @@ -26,11 +26,11 @@ function Goddard() vmax = 0.1 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] # variables tf ∈ R, variable - t ∈ [ t0, tf ], time + t ∈ [t0, tf], time x ∈ R³, state u ∈ R, control r = x₁ @@ -38,14 +38,14 @@ function Goddard() m = x₃ # constraints - 0 ≤ u(t) ≤ 1, (u_con) - r(t) ≥ r0, (x_con_rmin) - 0 ≤ v(t) ≤ vmax, (x_con_vmax) - x(t0) == x0, (initial_con) - m(tf) == mf, (final_con) + 0 ≤ u(t) ≤ 1, (u_con) + r(t) ≥ r0, (x_con_rmin) + 0 ≤ v(t) ≤ vmax, (x_con_vmax) + x(t0) == x0, (initial_con) + m(tf) == mf, (final_con) # dynamics - ẋ(t) == F0(x(t)) + u(t)*F1(x(t)) + ẋ(t) == F0(x(t)) + u(t) * F1(x(t)) # objective r(tf) → max @@ -54,12 +54,12 @@ function Goddard() # dynamics 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 # @@ -67,4 +67,4 @@ function Goddard() return ocp, objective -end \ No newline at end of file +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 0c9c5239..b5dbe243 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_periodic() T = halfperiod * 2 # ocp - t ∈ [ 0, T ], time + t ∈ [0, T], time x ∈ R³, state u ∈ R, control y = x[1] @@ -42,18 +42,25 @@ function bioreactor_1day_periodic() 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_periodic", init=nothing)) + return (( + ocp = bioreactor_1, + obj = 0.614134, + name = "bioreactor_1day_periodic", + init = nothing, + )) end @@ -73,7 +80,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 +88,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/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 df0237de..e3efae60 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,34 +22,45 @@ 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, lb=x0, ub=x0) - constraint!(goddard, :final, rg=3, lb=mf, ub=Inf) + time!(goddard, t0 = 0, indf = 1) + constraint!(goddard, :initial, lb = x0, ub = x0) + constraint!(goddard, :final, rg = 3, lb = mf, ub = Inf) 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 # abstratc definition -function goddard_a(;vmax=0.1, Tmax=3.5) +function goddard_a(; vmax = 0.1, Tmax = 3.5) # constants Cd = 310 beta = 500 @@ -58,11 +69,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 @@ -75,9 +86,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 f36c8b0c..d6e79537 100644 --- a/test/problems/parametric.jl +++ b/test/problems/parametric.jl @@ -1,29 +1,29 @@ # Parametric problem (name ??) function myocp(ρ) - + 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 param end # return problem and objective -return ((ocp=myocp(1.), obj=-0.335936, name="param")) \ No newline at end of file +return ((ocp = myocp(1.0), obj = -0.335936, name = "param")) 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/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/runtests.jl b/test/runtests.jl index 2039294b..0cd75a8c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,12 +21,12 @@ include("problems/goddard.jl") :grid, :initial_guess, :objective, - ) + ) @testset "$(name)" begin test_name = Symbol(:test_, name) - println("Testing: "*string(name)) + println("Testing: " * string(name)) include("$(test_name).jl") @eval $test_name() end end -end \ No newline at end of file +end diff --git a/test/test-nlp.jl b/test/test-nlp.jl index c7757579..ec322537 100644 --- a/test/test-nlp.jl +++ b/test/test-nlp.jl @@ -8,12 +8,12 @@ using Plots x ∈ R², state u ∈ R, control - x(0) == [ -1, 0 ] - x(1) == [ 0, 0 ] + x(0) == [-1, 0] + x(1) == [0, 0] - ẋ(t) == [ x₂(t), u(t) ] + ẋ(t) == [x₂(t), u(t)] - ∫( 0.5u(t)^2 ) → min + ∫(0.5u(t)^2) → min end @@ -21,8 +21,8 @@ docp, nlp = direct_transcription(ocp) #println(nlp.meta.x0) using NLPModelsIpopt -nlp_sol = ipopt(nlp; print_level=5, mu_strategy="adaptive", tol=1e-8, sb="yes") -sol = OptimalControlSolution(docp, primal=nlp_sol.solution, dual=nlp_sol.multipliers) +nlp_sol = ipopt(nlp; print_level = 5, mu_strategy = "adaptive", tol = 1e-8, sb = "yes") +sol = OptimalControlSolution(docp, primal = nlp_sol.solution, dual = nlp_sol.multipliers) plot(sol) using MadNLP @@ -32,10 +32,14 @@ set_initial_guess(docp, nlp, sol) #println(nlp.meta.x0) mad_nlp_sol = madnlp(nlp) -docp, nlp = direct_transcription(ocp, init=sol) +docp, nlp = direct_transcription(ocp, init = sol) mad_nlp_sol = madnlp(nlp) -mad_sol = OptimalControlSolution(docp, primal=madnlp_sol.solution, dual=madnlp_sol.multipliers) +mad_sol = OptimalControlSolution( + docp, + primal = madnlp_sol.solution, + dual = madnlp_sol.multipliers, +) plot(mad_sol) @@ -44,4 +48,4 @@ using Percival per_nlp_sol = percival(nlp, verbose = 1) per_sol = OptimalControlSolution(docp, primal=per_nlp_sol.solution, dual=per_nlp_sol.multipliers) plot(per_sol) -=# \ No newline at end of file +=# diff --git a/test/test_abstract_ocp.jl b/test/test_abstract_ocp.jl index 0a256342..1d63b9bd 100644 --- a/test/test_abstract_ocp.jl +++ b/test/test_abstract_ocp.jl @@ -1,29 +1,29 @@ # for parser testing function test_abstract_ocp() -# double integrator min tf, abstract definition -@def ocp1 begin - tf ∈ R, variable - 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) ] - tf → min -end + # double integrator min tf, abstract definition + @def ocp1 begin + tf ∈ R, variable + 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)] + tf → min + end -@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin - sol1 = solve(ocp1, print_level=0, tol=1e-12) - @test sol1.objective ≈ 2.0 rtol=1e-2 -end + @testset verbose = true showtiming = true ":double_integrator :min_tf :abstract" begin + sol1 = solve(ocp1, print_level = 0, tol = 1e-12) + @test sol1.objective ≈ 2.0 rtol = 1e-2 + end -@testset verbose = true showtiming = true ":goddard :max_rf :abstract :constr" begin - ocp3 = goddard_a().ocp - sol3 = solve(ocp3, print_level=0, tol=1e-12) - @test sol3.objective ≈ 1.0125 rtol=1e-2 -end + @testset verbose = true showtiming = true ":goddard :max_rf :abstract :constr" begin + ocp3 = goddard_a().ocp + sol3 = solve(ocp3, print_level = 0, tol = 1e-12) + @test sol3.objective ≈ 1.0125 rtol = 1e-2 + end -end \ No newline at end of file +end diff --git a/test/test_basic.jl b/test/test_basic.jl index 4e1d7f49..2cb689a8 100644 --- a/test/test_basic.jl +++ b/test/test_basic.jl @@ -1,16 +1,16 @@ function test_basic() @def ocp begin - t ∈ [ 0, 1 ], time + t ∈ [0, 1], time x ∈ R², state u ∈ R, control - x(0) == [ -1, 0 ] - x(1) == [ 0, 0 ] - ẋ(t) == [ x₂(t), u(t) ] - ∫( 0.5u(t)^2 ) → min + x(0) == [-1, 0] + x(1) == [0, 0] + ẋ(t) == [x₂(t), u(t)] + ∫(0.5u(t)^2) → min end - sol = solve(ocp; display=false) - @test sol.objective ≈ 6 atol=1e-2 + sol = solve(ocp; display = false) + @test sol.objective ≈ 6 atol = 1e-2 -end \ No newline at end of file +end diff --git a/test/test_continuation.jl b/test/test_continuation.jl index 6c9449be..416f1801 100644 --- a/test/test_continuation.jl +++ b/test/test_continuation.jl @@ -1,107 +1,114 @@ # discrete continuation function test_continuation() -test1 = true -test2 = true -test3 = true -draw_plot = false + test1 = true + test2 = true + test3 = true + draw_plot = false -# parametric ocp definition -if test1 - function ocp_T(T) - @def ocp begin - t ∈ [ 0, T ], time - x ∈ R², state - u ∈ R, control - q = x₁ - v = x₂ - q(0) == 0 - v(0) == 0 - q(T) == 1 - v(T) == 0 - ẋ(t) == [ v(t), u(t) ] - ∫(u(t)^2) → min + # parametric ocp definition + if test1 + function ocp_T(T) + @def ocp begin + t ∈ [0, T], time + x ∈ R², state + u ∈ R, control + q = x₁ + v = x₂ + q(0) == 0 + v(0) == 0 + q(T) == 1 + v(T) == 0 + ẋ(t) == [v(t), u(t)] + ∫(u(t)^2) → min + end + return ocp end - return ocp - end - @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin - init = () - obj_list = [] - for T=1:5 - ocp = ocp_T(T) - sol = solve(ocp, print_level=0, init=init) - init = sol - push!(obj_list, sol.objective) + @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin + init = () + obj_list = [] + for T = 1:5 + ocp = ocp_T(T) + sol = solve(ocp, print_level = 0, 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 end - @test obj_list ≈ [12, 1.5, 0.44, 0.19, 0.096] rtol=1e-2 end -end -# parametric ocp definition -if test2 - relu(x) = max(0, x) - μ = 10 - p_relu(x) = log(abs(1 + exp(μ*x)))/μ - f(x) = 1-x - m(x) = (p_relu∘f)(x) - T = 2 - function myocp(ρ) - @def ocp begin - τ ∈ R, variable - 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)] - -1 ≤ u₁(s) ≤ 1 - -1 ≤ u₂(s) ≤ 1 - 0 ≤ τ ≤ T - -(x₂(1)-2)^3 - ∫( ρ * ( τ*m(x₁(s))^2 + (T-τ)*m(x₂(s))^2 ) ) → min + # parametric ocp definition + if test2 + relu(x) = max(0, x) + μ = 10 + p_relu(x) = log(abs(1 + exp(μ * x))) / μ + f(x) = 1 - x + m(x) = (p_relu ∘ f)(x) + T = 2 + function myocp(ρ) + @def ocp begin + τ ∈ R, variable + 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)] + -1 ≤ u₁(s) ≤ 1 + -1 ≤ u₂(s) ≤ 1 + 0 ≤ τ ≤ T + -(x₂(1) - 2)^3 - ∫(ρ * (τ * m(x₁(s))^2 + (T - τ) * m(x₂(s))^2)) → min + end + return ocp end - return ocp - end - @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin - init = () - obj_list = [] - for ρ in [0.1, 5, 10, 30, 100] - ocp = myocp(ρ) - sol = solve(ocp, print_level=0, init=init) - init = sol - push!(obj_list, sol.objective) + @testset verbose = true showtiming = true ":parametric_ocp :warm_start" begin + init = () + obj_list = [] + for ρ in [0.1, 5, 10, 30, 100] + ocp = myocp(ρ) + sol = solve(ocp, print_level = 0, init = init) + init = sol + push!(obj_list, sol.objective) + end + @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol = 1e-2 end - @test obj_list ≈ [-0.034, -1.7, -6.2, -35, -148] rtol=1e-2 end -end -# parametric ocp definition -if test3 - sol0 = solve(goddard().ocp, print_level=0) + # parametric ocp definition + if test3 + sol0 = solve(goddard().ocp, print_level = 0) - @testset verbose = true showtiming = true ":global_variable :warm_start" begin - sol = sol0 - Tmax_list = [] - obj_list = [] - for Tmax=3.5:-0.5:1 - sol = solve(goddard(Tmax=Tmax).ocp, print_level=0, 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 + @testset verbose = true showtiming = true ":global_variable :warm_start" begin + sol = sol0 + Tmax_list = [] + obj_list = [] + for Tmax = 3.5:-0.5:1 + sol = solve(goddard(Tmax = Tmax).ocp, print_level = 0, 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 - if draw_plot - # plot obj(vmax) - 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))) + if draw_plot + # plot obj(vmax) + 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))) + end end end -end -end \ No newline at end of file +end diff --git a/test/test_goddard_direct.jl b/test/test_goddard_direct.jl index 21836f7c..d7f3f06b 100644 --- a/test/test_goddard_direct.jl +++ b/test/test_goddard_direct.jl @@ -1,11 +1,11 @@ function test_goddard_direct() -# goddard with state constraint - maximize altitude -ocp, objective = Goddard() + # goddard with state constraint - maximize altitude + ocp, objective = Goddard() -# initial guess (constant state and control functions) -init = (state=[1.01, 0.05, 0.8], control=0.1, variable=0.2) -sol = solve(ocp; grid_size=10, display=false, init=init) -@test sol.objective ≈ objective atol=5e-3 + # initial guess (constant state and control functions) + init = (state = [1.01, 0.05, 0.8], control = 0.1, variable = 0.2) + sol = solve(ocp; grid_size = 10, display = false, init = init) + @test sol.objective ≈ objective atol = 5e-3 -end \ No newline at end of file +end diff --git a/test/test_goddard_indirect.jl b/test/test_goddard_indirect.jl index e999b863..42d22fda 100644 --- a/test/test_goddard_indirect.jl +++ b/test/test_goddard_indirect.jl @@ -14,22 +14,22 @@ function test_goddard_indirect() vmax = 0.1 m0 = 1 mf = 0.6 - x0 = [ r0, v0, m0 ] + x0 = [r0, v0, m0] # remove_constraint!(ocp, :x_con_rmin) - g(x) = vmax-constraint(ocp, :x_con_vmax)(x, Real[]) # g(x, u) ≥ 0 (cf. nonnegative multiplier) - final_mass_cons(xf) = constraint(ocp, :final_con)(x0, xf, Real[])-mf + g(x) = vmax - constraint(ocp, :x_con_vmax)(x, Real[]) # g(x, u) ≥ 0 (cf. nonnegative multiplier) + final_mass_cons(xf) = constraint(ocp, :final_con)(x0, xf, Real[]) - mf 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 # bang controls @@ -39,13 +39,13 @@ function test_goddard_indirect() # singular control H0 = Lift(F0) H1 = Lift(F1) - H01 = @Lie {H0, H1} + H01 = @Lie {H0, H1} H001 = @Lie {H0, H01} H101 = @Lie {H1, H01} us(x, p) = -H001(x, p) / H101(x, p) # boundary control - ub(x) = -Lie(F0, g)(x) / Lie(F1, g)(x) + ub(x) = -Lie(F0, g)(x) / Lie(F1, g)(x) μ(x, p) = H01(x, p) / Lie(F1, g)(x) # flows @@ -62,7 +62,7 @@ function test_goddard_indirect() x3, p3 = fb(t2, x2, p2, t3) xf, pf = f0(t3, x3, p3, tf) s[1] = final_mass_cons(xf) - s[2:3] = pf[1:2] - [ 1, 0 ] + s[2:3] = pf[1:2] - [1, 0] s[4] = H1(x1, p1) s[5] = H01(x1, p1) s[6] = g(x2) @@ -80,17 +80,19 @@ function test_goddard_indirect() # test shooting function s = zeros(eltype(p0), 7) shoot!(s, p0, t1, t2, t3, tf) - s_guess_sol = [-0.02456074767656735, - -0.05699760226157302, - 0.0018629693253921868, - -0.027013078908634858, - -0.21558816838342798, - -0.0121146739026253, - 0.015713236406057297] - @test s ≈ s_guess_sol atol=1e-6 + s_guess_sol = [ + -0.02456074767656735, + -0.05699760226157302, + 0.0018629693253921868, + -0.027013078908634858, + -0.21558816838342798, + -0.0121146739026253, + 0.015713236406057297, + ] + @test s ≈ s_guess_sol atol = 1e-6 # - ξ0 = [ p0 ; t1 ; t2 ; t3 ; tf ] + ξ0 = [p0; t1; t2; t3; tf] # foo!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7]) @@ -103,7 +105,7 @@ function test_goddard_indirect() t1 = sol.u[4] t2 = sol.u[5] t3 = sol.u[6] - tf = sol.u[7]; + tf = sol.u[7] shoot!(s, p0, t1, t2, t3, tf) @@ -111,4 +113,4 @@ function test_goddard_indirect() @test SciMLBase.successful_retcode(sol) @test norm(s) < 1e-6 -end \ No newline at end of file +end diff --git a/test/test_grid.jl b/test/test_grid.jl index ae70c96f..67286226 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -5,75 +5,79 @@ function test_grid() ocp = Model() state!(ocp, 1) control!(ocp, 2) - time!(ocp, t0=0, tf=1) - constraint!(ocp, :initial, lb=-1, ub=-1) - constraint!(ocp, :final, lb=0, ub=0) - constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) + time!(ocp, t0 = 0, tf = 1) + constraint!(ocp, :initial, lb = -1, ub = -1) + constraint!(ocp, :final, lb = 0, ub = 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) - time_grid = [0,0.1,0.3,0.6,0.98,0.99,1] - sol = 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 = solve(ocp; time_grid = time_grid, display = false) + @test sol.objective ≈ 0.309 rtol = 1e-2 # 1. simple integrator min energy (dual control for test) ocp = Model() state!(ocp, 1) control!(ocp, 2) - time!(ocp, t0=0, tf=1) - constraint!(ocp, :initial, lb=-1, ub=-1) - constraint!(ocp, :final, lb=0, ub=0) - constraint!(ocp, :control, lb=[0,0], ub=[Inf, Inf]) + time!(ocp, t0 = 0, tf = 1) + constraint!(ocp, :initial, lb = -1, ub = -1) + constraint!(ocp, :final, lb = 0, ub = 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) - sol0 = solve(ocp, print_level=0) - + objective!(ocp, :lagrange, (x, u) -> (u[1] + u[2])^2) + sol0 = solve(ocp, print_level = 0) + # 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 = solve(ocp, time_grid=time_grid, print_level=0) + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1) + sol = solve(ocp, time_grid = time_grid, print_level = 0) @test sol.objective == sol0.objective @test 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 = solve(ocp, time_grid=time_grid, print_level=0) - @test sol.objective ≈ 0.309 rtol=1e-2 + time_grid = [0, 0.1, 0.3, 0.6, 0.98, 0.99, 1] + sol = solve(ocp, time_grid = time_grid, print_level = 0) + @test sol.objective ≈ 0.309 rtol = 1e-2 end - - + + # 2. integrator free times - 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, 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.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, 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.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) - sol0 = solve(ocp, print_level=0) - + sol0 = solve(ocp, print_level = 0) + @testset verbose = true showtiming = true ":explicit_grid" begin - sol = solve(ocp, time_grid=LinRange(0, 1, CTDirect.__grid_size()+1), print_level=0) - @test sol.objective == sol0.objective + sol = solve( + ocp, + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), + print_level = 0, + ) + @test sol.objective == sol0.objective @test sol.iterations == sol0.iterations end - + @testset verbose = true showtiming = true ":non_uniform_grid" begin - sol = solve(ocp, time_grid=[0,0.1,0.3,0.5,0.6,0.8,0.95,1], print_level=0) + sol = solve(ocp, time_grid = [0, 0.1, 0.3, 0.5, 0.6, 0.8, 0.95, 1], print_level = 0) #plot(sol, show=true) - @test sol.objective ≈ 7.96 rtol=1e-2 + @test sol.objective ≈ 7.96 rtol = 1e-2 end - + # 3. parametric ocp (T=2) with explicit / non-uniform grid @def ocpT2 begin - t ∈ [ 0, 2 ], time + t ∈ [0, 2], time x ∈ R², state u ∈ R, control q = x₁ @@ -82,21 +86,25 @@ function test_grid() v(0) == 0 q(2) == 1 v(2) == 0 - ẋ(t) == [ v(t), u(t) ] + ẋ(t) == [v(t), u(t)] ∫(u(t)^2) → min end - sol0 = solve(ocpT2, print_level=0) - + sol0 = solve(ocpT2, print_level = 0) + @testset verbose = true showtiming = true ":explicit_grid" begin - sol = solve(ocpT2, time_grid=LinRange(0, 1, CTDirect.__grid_size()+1), print_level=0) + sol = solve( + ocpT2, + time_grid = LinRange(0, 1, CTDirect.__grid_size() + 1), + print_level = 0, + ) @test sol.objective == sol0.objective @test sol.iterations == sol0.iterations end - + @testset verbose = true showtiming = true ":non_uniform_grid" begin - sol = solve(ocpT2, time_grid=[0,0.3,1,1.9,2],print_level=0) - @test sol.objective ≈ 2.43 rtol=1e-2 + sol = solve(ocpT2, time_grid = [0, 0.3, 1, 1.9, 2], print_level = 0) + @test sol.objective ≈ 2.43 rtol = 1e-2 end -end \ No newline at end of file +end diff --git a/test/test_initial_guess.jl b/test/test_initial_guess.jl index 92af44a7..c313a350 100644 --- a/test/test_initial_guess.jl +++ b/test/test_initial_guess.jl @@ -1,137 +1,198 @@ # test initial guess options function test_initial_guess() -# use 0 iterations to check initial guess, >0 to check cv -maxiter = 0 + # use 0 iterations to check initial guess, >0 to check cv + maxiter = 0 -# test functions -function check_xf(sol, xf) - return xf == sol.state(sol.times[end]) -end -function check_uf(sol, uf) - return uf == sol.control(sol.times[end]) -end -function check_v(sol, v) - return v == sol.variable -end + # test functions + function check_xf(sol, xf) + return xf == sol.state(sol.times[end]) + end + function check_uf(sol, uf) + return uf == sol.control(sol.times[end]) + end + function check_v(sol, v) + return v == sol.variable + end -# reference solution -@def ocp begin - tf ∈ R, variable - 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) ] - tf → min -end -sol0 = solve(ocp, print_level=0) + # reference solution + @def ocp begin + tf ∈ R, variable + 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)] + tf → min + end + sol0 = solve(ocp, print_level = 0) -# constant initial guess -x_const = [0.5, 0.2] -u_const = 0.5 -v_const = 0.15 + # constant initial guess + x_const = [0.5, 0.2] + 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 + # functional initial guess + 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] + # interpolated initial guess + x_vec = [[0, 0], [1, 2], [5, -1]] + x_matrix = [0 0; 1 2; 5 -1] + u_vec = [0, 0.3, 0.1] -################################################# -# 1 Pass initial guess to all-in-one solve call + ################################################# + # 1 Pass initial guess to all-in-one solve call -# 1.a default initial guess -@testset verbose = true showtiming = true ":default_init_no_arg" begin - sol = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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.a default initial guess + @testset verbose = true showtiming = true ":default_init_no_arg" begin + sol = solve(ocp, print_level = 0, 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 = solve(ocp, print_level = 0, 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 = solve(ocp, print_level = 0, 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 = solve(ocp, print_level=0, init=(state=x_const,), max_iter=maxiter) - @test(check_xf(sol, x_const)) -end -@testset verbose = true showtiming = true ":constant_u" begin - sol = solve(ocp, print_level=0, init=(control=u_const,), max_iter=maxiter) - @test(check_uf(sol, u_const)) -end -@testset verbose = true showtiming = true ":constant_v" begin - sol = solve(ocp, print_level=0, init=(variable=v_const,), max_iter=maxiter) - @test(check_v(sol, v_const)) -end -@testset verbose = true showtiming = true ":constant_xu" begin - sol = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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 + # 1.b constant initial guess + @testset verbose = true showtiming = true ":constant_x" begin + sol = solve(ocp, print_level = 0, init = (state = x_const,), max_iter = maxiter) + @test(check_xf(sol, x_const)) + end + @testset verbose = true showtiming = true ":constant_u" begin + sol = solve(ocp, print_level = 0, init = (control = u_const,), max_iter = maxiter) + @test(check_uf(sol, u_const)) + end + @testset verbose = true showtiming = true ":constant_v" begin + sol = solve(ocp, print_level = 0, init = (variable = v_const,), max_iter = maxiter) + @test(check_v(sol, v_const)) + end + @testset verbose = true showtiming = true ":constant_xu" begin + sol = solve( + ocp, + print_level = 0, + 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 = solve( + ocp, + print_level = 0, + 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 = solve( + ocp, + print_level = 0, + 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 = solve( + ocp, + print_level = 0, + 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 -# 1. functional initial guess -@testset verbose = true showtiming = true ":functional_x" begin - sol = solve(ocp, print_level=0, init=(state=x_func,), max_iter=maxiter) - @test(check_xf(sol, x_func(sol.times[end]))) -end -@testset verbose = true showtiming = true ":functional_u" begin - sol = solve(ocp, print_level=0, init=(control=u_func,), max_iter=maxiter) - @test(check_uf(sol, u_func(sol.times[end]))) -end -@testset verbose = true showtiming = true ":functional_xu" begin - sol = solve(ocp, print_level=0, init=(state=x_func, control=u_func), max_iter=maxiter) - @test(check_xf(sol, x_func(sol.times[end])) && check_uf(sol, u_func(sol.times[end]))) -end + # 1. functional initial guess + @testset verbose = true showtiming = true ":functional_x" begin + sol = solve(ocp, print_level = 0, init = (state = x_func,), max_iter = maxiter) + @test(check_xf(sol, x_func(sol.times[end]))) + end + @testset verbose = true showtiming = true ":functional_u" begin + sol = solve(ocp, print_level = 0, init = (control = u_func,), max_iter = maxiter) + @test(check_uf(sol, u_func(sol.times[end]))) + end + @testset verbose = true showtiming = true ":functional_xu" begin + sol = solve( + ocp, + print_level = 0, + init = (state = x_func, control = u_func), + max_iter = maxiter, + ) + @test( + check_xf(sol, x_func(sol.times[end])) && check_uf(sol, u_func(sol.times[end])) + ) + end -# 1.d interpolated initial guess -t_vec = [0, .1, v_const] -@testset verbose = true showtiming = true ":vector_txu :constant_v" begin - sol = solve(ocp, print_level=0, 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] -@testset verbose = true showtiming = true ":matrix_t :vector_xu :constant_v" begin - sol = solve(ocp, print_level=0, 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 = solve(ocp, print_level=0, 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.d interpolated initial guess + t_vec = [0, 0.1, v_const] + @testset verbose = true showtiming = true ":vector_txu :constant_v" begin + sol = solve( + ocp, + print_level = 0, + 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 0.1 v_const] + @testset verbose = true showtiming = true ":matrix_t :vector_xu :constant_v" begin + sol = solve( + ocp, + print_level = 0, + 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 = solve( + ocp, + print_level = 0, + 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 = solve(ocp, print_level=0, 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.times[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 = solve( + ocp, + print_level = 0, + 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.times[end])) && + check_v(sol, v_const) + ) + end -# 1.f warm start -@testset verbose = true showtiming = true ":warm_start" begin - sol = solve(ocp, print_level=0, init=sol0, max_iter=maxiter) - @test(check_xf(sol, sol.state(sol.times[end])) && check_uf(sol, sol.control(sol.times[end])) && check_v(sol, sol.variable)) -end + # 1.f warm start + @testset verbose = true showtiming = true ":warm_start" begin + sol = solve(ocp, print_level = 0, init = sol0, max_iter = maxiter) + @test( + check_xf(sol, sol.state(sol.times[end])) && + check_uf(sol, sol.control(sol.times[end])) && + check_v(sol, sol.variable) + ) + end -end \ No newline at end of file +end diff --git a/test/test_objective.jl b/test/test_objective.jl index 777a2d13..0ffa9a09 100644 --- a/test/test_objective.jl +++ b/test/test_objective.jl @@ -1,71 +1,71 @@ # test some objective options, variable tf function test_objective() -# min tf -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 1) -time!(ocp, t0=0, indf=1) -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.1, ub=10) -dynamics!(ocp, (x, u, v) -> [x[2], u]) -objective!(ocp, :mayer, (x0, xf, v) -> v) + # min tf + ocp = Model(variable = true) + state!(ocp, 2) + control!(ocp, 1) + variable!(ocp, 1) + time!(ocp, t0 = 0, indf = 1) + 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.1, ub = 10) + dynamics!(ocp, (x, u, v) -> [x[2], u]) + objective!(ocp, :mayer, (x0, xf, v) -> v) -@testset verbose = true showtiming = true ":min_tf :mayer" begin - sol = solve(ocp, print_level=0, tol=1e-12) - @test sol.objective ≈ 2.0 rtol=1e-2 -end + @testset verbose = true showtiming = true ":min_tf :mayer" begin + sol = solve(ocp, print_level = 0, tol = 1e-12) + @test sol.objective ≈ 2.0 rtol = 1e-2 + end -# min tf (lagrange) -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 1) -time!(ocp, t0=0, indf=1) -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.1, ub=10) -dynamics!(ocp, (x, u, v) -> [x[2], u]) -objective!(ocp, :lagrange, (x, u, v) -> 1) + # min tf (lagrange) + ocp = Model(variable = true) + state!(ocp, 2) + control!(ocp, 1) + variable!(ocp, 1) + time!(ocp, t0 = 0, indf = 1) + 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.1, ub = 10) + dynamics!(ocp, (x, u, v) -> [x[2], u]) + objective!(ocp, :lagrange, (x, u, v) -> 1) -@testset verbose = true showtiming = true ":min_tf :lagrange" begin - sol = solve(ocp, print_level=0, tol=1e-12) - @test sol.objective ≈ 2.0 rtol=1e-2 -end + @testset verbose = true showtiming = true ":min_tf :lagrange" begin + sol = solve(ocp, print_level = 0, tol = 1e-12) + @test sol.objective ≈ 2.0 rtol = 1e-2 + end -# max t0 (free t0 and tf) -ocp = Model(variable=true) -state!(ocp, 2) -control!(ocp, 1) -variable!(ocp, 2) -time!(ocp, ind0=1, indf=2) -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.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) + # max t0 (free t0 and tf) + ocp = Model(variable = true) + state!(ocp, 2) + control!(ocp, 1) + variable!(ocp, 2) + time!(ocp, ind0 = 1, indf = 2) + 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.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) -@testset verbose = true showtiming = true ":max_t0" begin - sol = solve(ocp, print_level=0, tol=1e-12) - @test sol.objective ≈ 8.0 rtol=1e-2 -end + @testset verbose = true showtiming = true ":max_t0" begin + sol = solve(ocp, print_level = 0, tol = 1e-12) + @test sol.objective ≈ 8.0 rtol = 1e-2 + end -@testset verbose = true showtiming = true ":max_t0 :explicit_grid" begin - sol = solve(ocp, time_grid=LinRange(0,1,101), print_level=0, tol=1e-12) - @test sol.objective ≈ 8.0 rtol=1e-2 -end + @testset verbose = true showtiming = true ":max_t0 :explicit_grid" begin + sol = solve(ocp, time_grid = LinRange(0, 1, 101), print_level = 0, tol = 1e-12) + @test sol.objective ≈ 8.0 rtol = 1e-2 + end -@testset verbose = true showtiming = true ":max_t0 :non_uniform_grid" begin - sol = solve(ocp, time_grid=[0,0.1,0.6,0.95,1], print_level=0, tol=1e-12) - @test sol.objective ≈ 7.48 rtol=1e-2 -end + @testset verbose = true showtiming = true ":max_t0 :non_uniform_grid" begin + sol = solve(ocp, time_grid = [0, 0.1, 0.6, 0.95, 1], print_level = 0, tol = 1e-12) + @test sol.objective ≈ 7.48 rtol = 1e-2 + end -end \ No newline at end of file +end