From c535c8379b847f30d8906fdf4e3f299e50e8d9b8 Mon Sep 17 00:00:00 2001 From: Pierre Martinon Date: Fri, 26 Jul 2024 12:58:58 +0200 Subject: [PATCH] problem library now defines functions (#186) * problem library now defines functions that return named tuples for problems and relevant data, including init * updated tests --- benchmark/benchmark.jl | 26 +++--- problems/beam.jl | 27 +++--- problems/bioreactor.jl | 96 +++++++++++++++++++++ problems/bioreactor_1day_periodic.jl | 56 ------------- problems/bioreactor_Ndays.jl | 54 ------------ problems/fuller.jl | 27 +++--- problems/goddard.jl | 100 ++++++++++++++++------ problems/insurance.jl | 119 +++++++++++++-------------- problems/jackson.jl | 43 +++++----- problems/robbins.jl | 36 ++++---- problems/swimmer.jl | 90 ++++++++++---------- problems/vanderpol.jl | 31 ++++--- test/suite/test_abstract_ocp.jl | 23 +----- test/suite/test_constraints.jl | 22 ++--- test/suite/test_continuation.jl | 10 +-- 15 files changed, 380 insertions(+), 380 deletions(-) create mode 100644 problems/bioreactor.jl delete mode 100644 problems/bioreactor_1day_periodic.jl delete mode 100644 problems/bioreactor_Ndays.jl diff --git a/benchmark/benchmark.jl b/benchmark/benchmark.jl index 2a5c80b3..ff78aa13 100644 --- a/benchmark/benchmark.jl +++ b/benchmark/benchmark.jl @@ -17,23 +17,25 @@ blas_config = LinearAlgebra.BLAS.lbt_get_config() ####################################################### # set parameters tol = 1e-8 -grid_size = 5000 +grid_size = 500 precompile = true @printf("Settings: tol=%g grid_size=%d precompile=%s\n\n", tol, grid_size, precompile) ####################################################### -# load examples -# +++ could add :init field in the problems ! +# load examples library +problem_path = pwd()*"/problems" +for problem_file in filter(contains(r".jl$"), readdir(problem_path; join=true)) + include(problem_file) +end + +# load problems for benchmark print("Loading problems: ") names_list = ["beam", "bioreactor_1day_periodic", "fuller", "goddard", "insurance", "jackson", "swimmer", "vanderpol"] problem_list = [] -problem_path = pwd()*"/problems" -for problem_file in filter(contains(r".jl$"), readdir(problem_path; join=true)) - ocp_data = include(problem_file) - if ocp_data.name in names_list - @printf("%s ", ocp_data.name) - push!(problem_list,ocp_data) - end +for problem_name in names_list + ocp_data = getfield(Main, Symbol(problem_name))() + @printf("%s ", ocp_data.name) + push!(problem_list,ocp_data) end println("") @@ -55,8 +57,8 @@ end t_list = [] println("Benchmark step") for problem in problem_list - t = @elapsed local sol = solve(problem[:ocp], display=false, linear_solver=linear_solver, grid_size=grid_size, tol=tol) - if !isapprox(sol.objective, problem[:obj], rtol=5e-2) + t = @elapsed local sol = solve(problem[:ocp], init=problem[:init], display=false, 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 @printf("%-30s completed in %6.2f s after %4d iterations\n",problem[:name],t,sol.iterations) diff --git a/problems/beam.jl b/problems/beam.jl index 413f77b4..ade2c6ab 100644 --- a/problems/beam.jl +++ b/problems/beam.jl @@ -1,15 +1,18 @@ # Beam example from bocop -@def beam begin - t ∈ [ 0, 1 ], time - x ∈ R², state - u ∈ R, control - x(0) == [0,1] - x(1) == [0,-1] - ẋ(t) == [x₂(t), u(t)] - 0 ≤ x₁(t) ≤ 0.1 - -10 ≤ u(t) ≤ 10 - ∫(u(t)^2) → min -end +function beam() -return ((ocp=beam, obj=8.898598, name="beam")) + @def beam begin + t ∈ [ 0, 1 ], time + x ∈ R², state + u ∈ R, control + x(0) == [0,1] + x(1) == [0,-1] + ẋ(t) == [x₂(t), u(t)] + 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 diff --git a/problems/bioreactor.jl b/problems/bioreactor.jl new file mode 100644 index 00000000..0c9c5239 --- /dev/null +++ b/problems/bioreactor.jl @@ -0,0 +1,96 @@ +# Bioreactor example from bocop + +# METHANE PROBLEM +# mu2 according to growth model +# mu according to light model +# time scale is [0,10] for 24h (day then night) + +# growth model MONOD +function growth(s, mu2m, 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 +end + +# 1 day periodic problem +function bioreactor_1day_periodic() + @def bioreactor_1 begin + # constants + beta = 1 + c = 2 + gamma = 1 + Ks = 0.05 + mu2m = 0.1 + mubar = 1 + r = 0.005 + halfperiod = 5 + T = halfperiod * 2 + + # ocp + t ∈ [ 0, T ], time + x ∈ R³, state + u ∈ R, control + y = x[1] + s = x[2] + b = x[3] + mu = light(t, halfperiod) * mubar + mu2 = growth(s(t), mu2m, Ks) + [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 + end + + return ((ocp=bioreactor_1, obj=0.614134, name="bioreactor_1day_periodic", init=nothing)) +end + + +# N days (non periodic) +function bioreactor_Ndays() + @def bioreactor_N begin + # constants + beta = 1 + c = 2 + gamma = 1 + halfperiod = 5 + Ks = 0.05 + mu2m = 0.1 + mubar = 1 + r = 0.005 + N = 30 + T = 10 * N + + # ocp + t ∈ [ 0, T ], time + x ∈ R³, state + u ∈ R, control + y = x[1] + s = x[2] + b = x[3] + mu = light(t, halfperiod) * mubar + mu2 = growth(s(t), mu2m, Ks) + [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 + end + + return ((ocp=bioreactor_N, obj=19.0745, init=(state=[50,50,50],), name="bioreactor_Ndays")) +end \ No newline at end of file diff --git a/problems/bioreactor_1day_periodic.jl b/problems/bioreactor_1day_periodic.jl deleted file mode 100644 index 22e86c3f..00000000 --- a/problems/bioreactor_1day_periodic.jl +++ /dev/null @@ -1,56 +0,0 @@ -# Bioreactor example from bocop - - - -function growth(s, mu2m, Ks) - # MONOD - return mu2m * s / (s+Ks) -end - -function daynightgrowth(time, halfperiod, mubar) - # light model: max^2 (0,sin) * mubar - # DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi] - pi = 3.141592653589793 - days = time / (halfperiod*2) - tau = (days - floor(days)) * 2*pi - mu = max(0,sin(tau))^2 * mubar - return mu -end - -#= METHANE PROBLEM -mu2 according to growth model, mu according to light model -time scale is [0,10] for 24h (day then night)=# -@def bioreactor_1 begin - # constants - beta = 1 - c = 2 - gamma = 1 - Ks = 0.05 - mu2m = 0.1 - mubar = 1 - r = 0.005 - T = 10 - - t ∈ [ 0, T ], time - x ∈ R³, state - u ∈ R, control - y = x[1] - s = x[2] - b = x[3] - mu = daynightgrowth(t, T/2, mubar) - mu2 = growth(s(t), mu2m, Ks) - [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 -end - -#sol=solve(bioreactor_1) -#plot(sol) - -return ((ocp=bioreactor_1, obj=0.614134, name="bioreactor_1day_periodic")) \ No newline at end of file diff --git a/problems/bioreactor_Ndays.jl b/problems/bioreactor_Ndays.jl deleted file mode 100644 index 51f05f87..00000000 --- a/problems/bioreactor_Ndays.jl +++ /dev/null @@ -1,54 +0,0 @@ -# Bioreactor example from bocop - -# constants -beta = 1 -c = 2 -gamma = 1 -halfperiod = 5 -Ks = 0.05 -mu2m = 0.1 -mubar = 1 -r = 0.005 -pi = 3.141592653589793 -T = 300 - -function growth(s) - # MONOD - return mu2m * s / (s+Ks) -end - -function daynightgrowth(time) - # light model: max^2 (0,sin) * mubar - # DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi] - days = time / (halfperiod*2) - tau = (days - floor(days)) * 2*pi - mu = max(0,sin(tau))^2 * mubar - return mu -end - -#= METHANE PROBLEM -mu2 according to growth model, mu according to light model -time scale is [0,10] for 24h (day then night)=# -@def bioreactor_N begin - t ∈ [ 0, T ], time - x ∈ R³, state - u ∈ R, control - y = x[1] - s = x[2] - b = x[3] - mu = daynightgrowth(t) - mu2 = growth(s(t)) - [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 -end - -#sol=solve(bioreactor_N, grid_size=1000, init=(state=[50,50,50],)) -#plot(sol) -return ((ocp=bioreactor_N, obj=19.0745, name="bioreactor_Ndays")) \ No newline at end of file diff --git a/problems/fuller.jl b/problems/fuller.jl index ea440351..1740f118 100644 --- a/problems/fuller.jl +++ b/problems/fuller.jl @@ -1,17 +1,16 @@ # Fuller example +function fuller() -@def fuller begin - 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₁(t)^2) → min -end + @def fuller begin + 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₁(t)^2) → min + end -#sol = solve(fuller, grid_size=1000) -#plot(sol) - -return((ocp=fuller, obj=2.683944e-1, name="fuller")) \ No newline at end of file + return((ocp=fuller, obj=2.683944e-1, name="fuller", init=nothing)) +end \ No newline at end of file diff --git a/problems/goddard.jl b/problems/goddard.jl index 7c7e8825..df0237de 100644 --- a/problems/goddard.jl +++ b/problems/goddard.jl @@ -1,37 +1,83 @@ # Goddard problem # free final time, fixed final mass, max altitude +# constraint on max speed -goddard = Model(variable=true) -Cd = 310 -Tmax = 3.5 -β = 500 -b = 2 -r0 = 1 -v0 = 0 -vmax = 0.1 -m0 = 1 -mf = 0.6 -x0 = [ r0, v0, m0 ] -state!(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=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.01, ub=Inf) -objective!(goddard, :mayer, (x0, xf, v) -> xf[1], :max) -function F0(x) +# aux functions +function F0(x, Cd, beta) r, v, m = x - D = Cd * v^2 * exp(-β*(r - 1)) + D = Cd * v^2 * exp(-beta*(r - 1)) return [ v, -D/m - 1/r^2, 0 ] end -function F1(x) +function F1(x, Tmax, b) r, v, m = x return [ 0, Tmax/m, -b*Tmax ] end -dynamics!(goddard, (x, u, v) -> F0(x) + u*F1(x) ) -# return problem and objective -return ((ocp=goddard, obj=1.01257, name="goddard")) \ No newline at end of file +function goddard(;vmax=0.1, Tmax=3.5, functional_constraints=false) + # constants + Cd = 310 + beta = 500 + b = 2 + r0 = 1 + v0 = 0 + m0 = 1 + mf = 0.6 + x0 = [ r0, v0, m0 ] + + #ocp + 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) + 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) + else + 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) ) + + 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) + # constants + Cd = 310 + beta = 500 + b = 2 + r0 = 1 + v0 = 0 + m0 = 1 + mf = 0.6 + x0 = [ r0, v0, m0 ] + + @def goddard_a begin + tf ∈ R, variable + t ∈ [ 0, tf ], time + x ∈ R^3, state + u ∈ R, control + 0.1 ≤ tf ≤ Inf + r = x[1] + v = x[2] + m = x[3] + x(0) == x0 + m(tf) == mf + r0 ≤ r(t) ≤ 1.1 + 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) + 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 diff --git a/problems/insurance.jl b/problems/insurance.jl index c1933e7a..4bbcac0d 100644 --- a/problems/insurance.jl +++ b/problems/insurance.jl @@ -1,61 +1,60 @@ # Insurance (non audit) example from Bocop - -@def insurance begin - - # constants - gamma = 0.2 - lambda = 0.25 - h0 = 1.5 - w = 1 - s = 10 - k = 0 - sigma = 0 - alpha = 4 - tf = 10 - - t ∈ [ 0, tf ], time - x ∈ R³, state - u ∈ R⁵, control - P ∈ R, variable - I = x[1] # Insurance - m = x[2] # Expense - h = u[1] - R = u[2] # Revenue - H = u[3] # Health - U = u[4] # Utility - dUdR = u[5] - - 0 ≤ I(t) ≤ 1.5 - 0 ≤ m(t) ≤ 1.5 - 0 ≤ h(t) ≤ 25 - 0 ≤ R(t) ≤ Inf - 0 ≤ H(t) ≤ Inf - 0 ≤ U(t) ≤ Inf - 0.001 ≤ dUdR(t) ≤ Inf - 0 ≤ P ≤ Inf - - x(0) == [0, 0.001, 0] - P - x[3](tf) == 0 - - epsilon = k * t / (tf - t + 1) - # illness distribution - 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 - - 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 - - ẋ(t) == [(1-gamma*t*vprime/dUdR(t))*h(t), - h(t), - (1+sigma)*I(t)*fx] - ∫(U(t)*fx) → max -end - -#sol = solve(insurance, grid_size=1000) -#plot(sol) - -return((ocp=insurance, obj=2.059511, name="insurance")) \ No newline at end of file +function insurance() + + @def insurance begin + + # constants + gamma = 0.2 + lambda = 0.25 + h0 = 1.5 + w = 1 + s = 10 + k = 0 + sigma = 0 + alpha = 4 + tf = 10 + + t ∈ [ 0, tf ], time + x ∈ R³, state + u ∈ R⁵, control + P ∈ R, variable + I = x[1] # Insurance + m = x[2] # Expense + h = u[1] + R = u[2] # Revenue + H = u[3] # Health + U = u[4] # Utility + dUdR = u[5] + + 0 ≤ I(t) ≤ 1.5 + 0 ≤ m(t) ≤ 1.5 + 0 ≤ h(t) ≤ 25 + 0 ≤ R(t) ≤ Inf + 0 ≤ H(t) ≤ Inf + 0 ≤ U(t) ≤ Inf + 0.001 ≤ dUdR(t) ≤ Inf + 0 ≤ P ≤ Inf + + x(0) == [0, 0.001, 0] + P - x[3](tf) == 0 + + epsilon = k * t / (tf - t + 1) + # illness distribution + 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 + + 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 + + ẋ(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 diff --git a/problems/jackson.jl b/problems/jackson.jl index 682a0915..10421eb8 100644 --- a/problems/jackson.jl +++ b/problems/jackson.jl @@ -1,26 +1,25 @@ # Jackson example from Bocop +function jackson() -@def jackson begin - # constants - k1 = 1 - k2 = 10 - k3 = 1 + @def jackson begin + # constants + k1 = 1 + k2 = 10 + k3 = 1 - t ∈ [ 0, 4 ], time - x ∈ R³, state - u ∈ R, control - [0, 0, 0] ≤ x(t) ≤ [1.1, 1.1, 1.1] - 0 ≤ u(t) ≤ 1 - 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)] - x[3](4) → max -end + t ∈ [ 0, 4 ], time + x ∈ R³, state + u ∈ R, control + [0, 0, 0] ≤ x(t) ≤ [1.1, 1.1, 1.1] + 0 ≤ u(t) ≤ 1 + 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)] + x[3](4) → max + end -#sol = solve(jackson) -#plot(sol) - -return ((ocp=jackson, obj=0.192011, name="jackson")) \ No newline at end of file + return ((ocp=jackson, obj=0.192011, name="jackson", init=nothing)) +end \ No newline at end of file diff --git a/problems/robbins.jl b/problems/robbins.jl index fb02b447..73cd18ea 100644 --- a/problems/robbins.jl +++ b/problems/robbins.jl @@ -1,22 +1,20 @@ # Robbins example from Bocop +function robbins() + @def robbins begin + # constants + alpha = 3 + beta = 0 + gamma = 0.5 -@def robbins begin - # constants - alpha = 3 - beta = 0 - gamma = 0.5 + 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 + end - 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 -end - -#sol = solve(robbins) -#plot(sol) - -return ((ocp=robbins, obj=20., name="robbins")) \ No newline at end of file + return ((ocp=robbins, obj=20., name="robbins", init=nothing)) +end \ No newline at end of file diff --git a/problems/swimmer.jl b/problems/swimmer.jl index e9a14ffc..3f7fc1a6 100644 --- a/problems/swimmer.jl +++ b/problems/swimmer.jl @@ -1,61 +1,63 @@ # Microswimmer example from Bocop -@def swimmer begin +# +++ make 2 versions: 1 stroke periodic and free N strokes - tf = 25 - 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] +function swimmer() - # initial conditions - x[1:2](0) == [0, 0] - -3.15 ≤ x[3](0) ≤ 0 - 0 ≤ x[4](0) ≤ Inf # to break symmetry + @def swimmer begin + + tf = 25 + t ∈ [ 0, tf ], time + x ∈ R^5, state + u ∈ R^2, control - # final conditions - #x[1](tf) ≤ -0.5 # target displacement for min energy problem - x[2](tf) == 0 + # 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] + + # initial conditions + x[1:2](0) == [0, 0] + -3.15 ≤ x[3](0) ≤ 0 + 0 ≤ x[4](0) ≤ Inf # to break symmetry - # periodicity - #x[3:5](tf) - x[3:5](0) == [0, 0, 0] - #x[3:5](tf) - x[3:5](0) == 0 NB. WITH JUST 0 PARSE OK BUT DIM ERROR LATER + # final conditions + #x[1](tf) ≤ -0.5 # target displacement for min energy problem + x[2](tf) == 0 - #coefficients for dynamics - th = x[3](t) - b1 = x[4](t) - b3 = x[5](t) - u1 = u[1](t) - u2 = u[2](t) + # periodicity + #x[3:5](tf) - x[3:5](0) == [0, 0, 0] + #x[3:5](tf) - x[3:5](0) == 0 NB. WITH JUST 0 PARSE OK BUT DIM ERROR LATER - 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) + #coefficients for dynamics + th = x[3](t) + b1 = x[4](t) + b3 = x[5](t) + u1 = u[1](t) + u2 = u[2](t) - 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) + 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) - 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) + 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) - 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) + 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) - 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) + 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) - 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) + 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) - 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) + 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) - ẋ(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 + 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) -#sol = solve(swimmer, grid_size=100) -#plot(sol) + ẋ(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")) \ No newline at end of file + return((ocp=swimmer, obj=0.984273, name="swimmer", init=nothing)) +end \ No newline at end of file diff --git a/problems/vanderpol.jl b/problems/vanderpol.jl index c716aeb9..1f82facd 100644 --- a/problems/vanderpol.jl +++ b/problems/vanderpol.jl @@ -1,20 +1,19 @@ # Van der Pol example from Bocop +function vanderpol() -@def vanderpol begin - # constants - omega = 1 - epsilon = 1 + @def vanderpol begin + # constants + omega = 1 + epsilon = 1 - 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 -end + 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 + end -#sol = solve(vanderpol) -#plot(sol) - -return((ocp=vanderpol, obj=1.047921, name="vanderpol")) \ No newline at end of file + return((ocp=vanderpol, obj=1.047921, name="vanderpol", init=nothing)) +end \ No newline at end of file diff --git a/test/suite/test_abstract_ocp.jl b/test/suite/test_abstract_ocp.jl index 3a82dafd..11aec27d 100644 --- a/test/suite/test_abstract_ocp.jl +++ b/test/suite/test_abstract_ocp.jl @@ -19,27 +19,8 @@ end @test sol1.objective ≈ 2.0 rtol=1e-2 end -# goddard -# NB. the ≤ is not the same as <= (parse error for <=) -@def ocp3 begin - tf ∈ R, variable - t ∈ [ 0, tf ], time - x ∈ R^3, state - u ∈ R, control - 0.1 ≤ tf ≤ Inf - r = x[1] - v = x[2] - m = x[3] - x(0) == [1, 0, 1] - m(tf) == 0.6 - 1 ≤ r(t) ≤ 1.1 - 0 ≤ v(t) ≤ vmax - 0 ≤ u(t) ≤ 1 - ẋ(t) == F0(x(t)) + u(t)*F1(x(t)) - r(tf) → max -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 \ No newline at end of file +end diff --git a/test/suite/test_constraints.jl b/test/suite/test_constraints.jl index 38c69a7b..3cb71ae5 100644 --- a/test/suite/test_constraints.jl +++ b/test/suite/test_constraints.jl @@ -1,27 +1,15 @@ println("Test: constraint types") +# box constraints @testset verbose = true showtiming = true ":goddard :box_constraints" begin - sol = solve(goddard, print_level=0, tol=1e-8) + ocp = goddard() + sol = solve(ocp.ocp, print_level=0, tol=1e-8) @test sol.objective ≈ 1.0125 rtol=1e-2 end # functional constraints -ocp1 = Model(variable=true) -state!(ocp1, 3) -control!(ocp1, 1) -variable!(ocp1, 1) -time!(ocp1, t0=0, indf=1) -constraint!(ocp1, :initial, lb=x0, ub=x0) -constraint!(ocp1, :final, rg=3, lb=mf, ub=Inf) -constraint!(ocp1, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0]) -constraint!(ocp1, :control, f=(u,v)->u, lb=0, ub=1) -constraint!(ocp1, :variable, f=v->v, lb=0.01, ub=Inf) -objective!(ocp1, :mayer, (x0, xf, v) -> xf[1], :max) -dynamics!(ocp1, (x, u, v) -> F0(x) + u*F1(x) ) - -# note: the equations do not handle r<1 well -# without the box constraint on x, the default init (0.1) is not suitable @testset verbose = true showtiming = true ":goddard :functional_constraints" begin - sol1 = solve(ocp1, print_level=0, tol=1e-8, init=(state=[1.01,0.05,0.75],)) + ocp1 = goddard(functional_constraints=true) + sol1 = solve(ocp1.ocp, print_level=0, tol=1e-8, init=ocp1.init) @test sol1.objective ≈ 1.0125 rtol=1e-2 end \ No newline at end of file diff --git a/test/suite/test_continuation.jl b/test/suite/test_continuation.jl index 3a66b959..3982e084 100644 --- a/test/suite/test_continuation.jl +++ b/test/suite/test_continuation.jl @@ -77,18 +77,16 @@ if test2 end -# global variable used in ocp +# parametric ocp definition if test3 - Tmax = 3.5 - sol0 = solve(goddard, print_level=0) + 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_local=3.5:-0.5:1 - global Tmax = Tmax_local - sol = solve(goddard, print_level=0, init=sol) + 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