diff --git a/Project.toml b/Project.toml index cc055a7..94861e6 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "0.1.0" [deps] CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" [weakdeps] CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5" @@ -18,6 +17,5 @@ OptimalControlModels = "CTDirect" [compat] CTBase = "0.13" CTDirect = "0.12" -DataFrames = "1" JuMP = "1" julia = "1.10" diff --git a/docs/src/list_of_problems.md b/docs/src/list_of_problems.md index ba11457..0db57ed 100644 --- a/docs/src/list_of_problems.md +++ b/docs/src/list_of_problems.md @@ -11,9 +11,14 @@ For each problem, we provide the following data in [MetaData](https://github.com - `minimize::Bool`: true if optimize == minimize To get the list of metadata, you can use the following code: -```julia +```@example metadata using OptimalControlProblems OptimalControlProblems.metadata +nothing # hide +``` +To access the metadata of a specific problem, you can execute the following command: +```@example metadata +OptimalControlProblems.metadata[:chain] ``` ## List of Problems @@ -21,19 +26,23 @@ The table below summarizes the status of the each problem: | Problem | With JuMP | With OptimalControl | | --- | --- | --- | +| `beam` | ✅ | ✅| +| `bioreactor` | ✅ | ✅| +| `cart_pendulum` | ✅ | ✅| | `chain` | ✅ | ✅| +| `dielectrophoretic_particle` | ✅ | ✅| +| `double_oscillator` | ✅ | ✅| +| `ducted_fan` | ✅ | ✅ +| `electrical_vehicle` | ✅ | ✅| | `glider` | ✅ | ✅ | +| `insurance` | ✅ | ✅| +| `jackson` | ✅ | ✅| +| `moonlander` | ✅ | ❌| +| `quadrotor` | ✅ | ❌| +| `robbins` | ✅ | ✅| | `robot` | ✅ | ✅| | `rocket` | ✅ | ✅| -| `steering` | ✅ | ✅| | `space_shuttle` | ✅ | ❌| -| `cart_pendulum` | ✅ | ✅| -| `moonlander` | ✅ | ❌| +| `steering` | ✅ | ✅| | `truck_trailer` | ❌ | ❌| -| `quadrotor` | ✅ | ❌| -| `dielectrophoretic_particle` | ✅ | ✅| -| `ducted_fan` | ✅ | ✅| -| `double_oscillator` | ✅ | ✅| -| `electrical_vehicle` | ✅ | ✅| | `vanderpol` | ✅ | ✅| - diff --git a/ext/JuMPModels/beam.jl b/ext/JuMPModels/beam.jl new file mode 100644 index 0000000..55f22e3 --- /dev/null +++ b/ext/JuMPModels/beam.jl @@ -0,0 +1,37 @@ +""" +The Beam Problem: + The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.beam(::JuMPBackend;nh::Int=100) + # Parameters + tf = 1 + step = tf/nh + + # Model + model = JuMP.Model() + + @variables(model, begin + 0.0 <= x1[0:nh] <= 0.1, (start = 0.0) + x2[0:nh] , (start = 0.0) + -10.0 <= u[0:nh] <= 10.0 , (start = 0.0) + end) + + # Boundary constraints + @constraints(model, begin + x1[0] == 0 + x2[0] == 1 + x1[nh] == 0 + x2[nh] == -1 + end) + + # Dynamics + @constraints(model, begin + con_x1[t=1:nh], x1[t] == x1[t-1] + 0.5 * step * (x2[t] + x2[t-1]) + con_x2[t=1:nh], x2[t] == x2[t-1] + 0.5 * step * (u[t] + u[t-1]) + end) + + + @objective(model, Min, sum(u[t]^2 for t in 0:nh)) + + return model +end \ No newline at end of file diff --git a/ext/JuMPModels/bioreactor.jl b/ext/JuMPModels/bioreactor.jl new file mode 100644 index 0000000..e343c3c --- /dev/null +++ b/ext/JuMPModels/bioreactor.jl @@ -0,0 +1,57 @@ +""" +The Bioreactor Problem: + The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.bioreactor(::JuMPBackend;nh::Int=100, N::Int=30) + # Parameters + beta = 1 + c = 2 + gamma = 1 + halfperiod = 5 + Ks = 0.05 + mu2m = 0.1 + mubar = 1 + r = 0.005 + T = 10 * N + + # Model + model = JuMP.Model() + + @variables(model, begin + y[0:nh] >= 0, (start = 50) + s[0:nh] >= 0, (start = 50) + b[0:nh] >= 0.001, (start = 50) + 0 <= u[0:nh] <= 1 , (start = 0.5) + end) + + # Boundary constraints + @constraints(model, begin + 0.05 <= y[0] <= 0.25 + 0.5 <= s[0] <= 5 + 0.5 <= b[0] <= 3 + end) + + # Dynamics + @expressions(model, begin + step, T/nh + # intermediate variables + mu2[t=0:nh], mu2m * s[t] / (s[t] + Ks) + days[t=0:nh], (t * step) / (halfperiod * 2) + tau[t=0:nh], (days[t] - floor(days[t])) * 2 * pi + light[t=0:nh], max(0, sin(tau[t]))^2 + mu[t=0:nh], light[t] * mubar + # dynamics + dy[t=0:nh], mu[t] * y[t] / (1 + y[t]) - (r + u[t]) * y[t] + ds[t=0:nh], -mu2[t] * b[t] + u[t] * beta * (gamma * y[t] - s[t]) + db[t=0:nh], (mu2[t] - u[t] * beta) * b[t] + end) + @constraints(model, begin + con_y[t=1:nh], y[t] == y[t - 1] + 0.5 * step * (dy[t] + dy[t - 1]) + con_s[t=1:nh], s[t] == s[t - 1] + 0.5 * step * (ds[t] + ds[t - 1]) + con_b[t=1:nh], b[t] == b[t - 1] + 0.5 * step * (db[t] + db[t - 1]) + end) + + @objective(model, Max, sum(b[t] / (beta + c) for t in 0:nh)) + + return model +end \ No newline at end of file diff --git a/ext/JuMPModels/insurance.jl b/ext/JuMPModels/insurance.jl new file mode 100644 index 0000000..34f2153 --- /dev/null +++ b/ext/JuMPModels/insurance.jl @@ -0,0 +1,70 @@ +""" +The Insurance Problem: + The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.insurance(::JuMPBackend;nh::Int=100, N::Int=30) + # constants + gamma = 0.2 + lambda = 0.25 + h0 = 1.5 + w = 1 + s = 10 + k = 0 + sigma = 0 + alpha = 4 + tf = 10 + + # Model + model = JuMP.Model() + + @variables(model, begin + 0 <= I[0:nh] <= 1.5, (start = 0.1) + 0 <= m[0:nh] <= 1.5, (start = 0.1) + 0 <= x3[0:nh] <= 1.5, (start = 0.1) + 0 <= h[0:nh] <= 25, (start = 0.1) + 0 <= R[0:nh], (start = 0.1) + 0 <= H[0:nh], (start = 0.1) + 0 <= U[0:nh], (start = 0.1) + 0.001 <= dUdR[0:nh], (start = 0.1) + P >= 0, (start = 0.1) + end) + + # Boundary constraints + @constraints(model, begin + I[0] == 0 + m[0] == 0.001 + x3[0] == 0 + P - x3[nh] == 0 + end) + + @expressions(model, begin + step, tf/nh + epsilon[t=0:nh], k * t * step / (tf - t * step + 1) + fx[t=0:nh], lambda * exp(-lambda * t * step) + exp(-lambda * tf) / tf + v[t=0:nh], m[t]^(alpha / 2) / (1 + m[t]^(alpha / 2)) + vprime[t=0:nh], alpha / 2 * m[t]^(alpha / 2 - 1) / (1 + m[t]^(alpha / 2))^2 + end) + @constraints(model, begin + cond1[t=0:nh], R[t] - (w - P + I[t] - m[t] - epsilon[t]) == 0 + cond2[t=0:nh], H[t] - (h0 - gamma * step * t *(1 - v[t])) == 0 + cond3[t=0:nh], U[t] - (1 - exp(-s * R[t]) + H[t]) == 0 + cond4[t=0:nh], dUdR[t] - (s * exp(-s * R[t])) == 0 + end) + + # Dynamics + @expressions(model, begin + dI[t=0:nh], (1 - gamma * t * step * vprime[t] / dUdR[t]) * h[t] + dm[t=0:nh], h[t] + dx3[t=0:nh], (1 + sigma) * I[t] * fx[t] + end) + @constraints(model, begin + con_dI[t=1:nh], I[t] == I[t - 1] + 0.5 * step * (dI[t] + dI[t - 1]) + con_dm[t=1:nh], m[t] == m[t - 1] + 0.5 * step * (dm[t] + dm[t - 1]) + con_dx3[t=1:nh], x3[t] == x3[t - 1] + 0.5 * step * (dx3[t] + dx3[t - 1]) + end) + + # Objective + @objective(model, Max, sum(U[t] * fx[t] for t in 0:nh)) + + return model +end \ No newline at end of file diff --git a/ext/JuMPModels/jackson.jl b/ext/JuMPModels/jackson.jl new file mode 100644 index 0000000..eaed342 --- /dev/null +++ b/ext/JuMPModels/jackson.jl @@ -0,0 +1,47 @@ +""" +The Jackson Problem: + The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.jackson(::JuMPBackend;nh::Int=100, N::Int=30) + # constants + k1 = 1 + k2 = 10 + k3 = 1 + tf = 4 + step = tf/nh + + # Model + model = JuMP.Model() + + @variables(model, begin + 0 <= a[0:nh] <= 1.1, (start = 0.1) + 0 <= b[0:nh] <= 1.1, (start = 0.1) + 0 <= x3[0:nh] <= 1.1, (start = 0.1) + 0 <= u[0:nh] <= 1, (start = 0.1) + end) + + # Boundary constraints + @constraints(model, begin + a[0] == 1 + b[0] == 0 + x3[0] == 0 + end) + + # Dynamics + @expressions(model, begin + da[t=0:nh], -u[t] * (k1 * a[t] - k2 * b[t]) + db[t=0:nh], u[t] * (k1 * a[t] - k2 * b[t]) - (1 - u[t]) * k3 * b[t] + dx3[t=0:nh], (1 - u[t]) * k3 * b[t] + end) + @constraints(model, begin + con_da[t=1:nh], a[t] == a[t - 1] + 0.5 * step * (da[t] + da[t - 1]) + con_db[t=1:nh], b[t] == b[t - 1] + 0.5 * step * (db[t] + db[t - 1]) + con_dx3[t=1:nh], x3[t] == x3[t - 1] + 0.5 * step * (dx3[t] + dx3[t - 1]) + end) + + # Objective + @objective(model, Max, x3[nh]) + + return model + +end \ No newline at end of file diff --git a/ext/JuMPModels/robbins.jl b/ext/JuMPModels/robbins.jl new file mode 100644 index 0000000..d1413d1 --- /dev/null +++ b/ext/JuMPModels/robbins.jl @@ -0,0 +1,45 @@ +""" +The Robbins Problem: + The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.robbins(::JuMPBackend;nh::Int=100, N::Int=30) + # constants + alpha = 3 + beta = 0 + gamma = 0.5 + tf = 10 + step = tf/nh + + # Model + model = JuMP.Model() + + @variables(model, begin + 0 <= x1[0:nh], (start = 0.1) + x2[0:nh], (start = 0.1) + x3[0:nh], (start = 0.1) + u[0:nh], (start = 0.1) + end) + + # Boundary constraints + @constraints(model, begin + x1[0] == 1 + x2[0] == -2 + x3[0] == 0 + x1[nh] == 0 + x2[nh] == 0 + x3[nh] == 0 + end) + + # Dynamics + @constraints(model, begin + con_x1[t=1:nh], x1[t] == x1[t - 1] + 0.5 * step * (x2[t] + x2[t - 1]) + con_x2[t=1:nh], x2[t] == x2[t - 1] + 0.5 * step * (x3[t] + x3[t - 1]) + con_dx3[t=1:nh], x3[t] == x3[t - 1] + 0.5 * step * (u[t] + u[t - 1]) + end) + + # Objective + @objective(model, Min, sum(alpha * x1[t] + beta * x1[t]^2 + gamma * u[t]^2 for t in 0:nh)) + + return model + +end \ No newline at end of file diff --git a/ext/MetaData/beam.jl b/ext/MetaData/beam.jl new file mode 100644 index 0000000..056532f --- /dev/null +++ b/ext/MetaData/beam.jl @@ -0,0 +1 @@ +beam_meta = Dict( :name => "beam", :nh => 100, :nvar => 404, :ncon => 305, :minimize => true,) diff --git a/ext/MetaData/bioreactor.jl b/ext/MetaData/bioreactor.jl new file mode 100644 index 0000000..de73bed --- /dev/null +++ b/ext/MetaData/bioreactor.jl @@ -0,0 +1 @@ +bioreactor_meta = Dict(:name => "bioreactor", :nh => 100, :nvar => 505, :ncon => 404, :minimize => false) diff --git a/ext/MetaData/insurance.jl b/ext/MetaData/insurance.jl new file mode 100644 index 0000000..ba6fa5c --- /dev/null +++ b/ext/MetaData/insurance.jl @@ -0,0 +1 @@ +insurance_meta = Dict(:name => "insurance", :nh => 100, :nvar => 910, :ncon => 809, :minimize => false) \ No newline at end of file diff --git a/ext/MetaData/jackson.jl b/ext/MetaData/jackson.jl new file mode 100644 index 0000000..15f97b2 --- /dev/null +++ b/ext/MetaData/jackson.jl @@ -0,0 +1 @@ +jackson_meta = Dict(:name => "jackson", :nh => 100, :nvar => 404, :ncon => 303, :minimize => false) \ No newline at end of file diff --git a/ext/MetaData/robbins.jl b/ext/MetaData/robbins.jl new file mode 100644 index 0000000..b1b4202 --- /dev/null +++ b/ext/MetaData/robbins.jl @@ -0,0 +1 @@ +robbins_meta = Dict(:name => "robbins", :nh => 100, :nvar => 505, :ncon => 407, :minimize => true) \ No newline at end of file diff --git a/ext/OptimalControlModels/beam.jl b/ext/OptimalControlModels/beam.jl new file mode 100644 index 0000000..f438a97 --- /dev/null +++ b/ext/OptimalControlModels/beam.jl @@ -0,0 +1,27 @@ +""" +The Beam Problem: + The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.beam(::OptimalControlBackend;nh::Int=100) + # Model + @def ocp 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 + + # Initial guess + init = (state = [0.0, 0.0], control = 0.0,) + + # NLPModel + nlp = direct_transcription(ocp ,init = init, grid_size = nh)[2] + + return nlp + +end \ No newline at end of file diff --git a/ext/OptimalControlModels/bioreactor.jl b/ext/OptimalControlModels/bioreactor.jl new file mode 100644 index 0000000..2e52b2f --- /dev/null +++ b/ext/OptimalControlModels/bioreactor.jl @@ -0,0 +1,73 @@ +""" +The Bioreactor Problem: + The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.bioreactor(::OptimalControlBackend; nh::Int=100, N::Int=30) + # constants + beta = 1 + c = 2 + gamma = 1 + halfperiod = 5 + Ks = 0.05 + mu2m = 0.1 + mubar = 1 + r = 0.005 + T = 10 * N + + # Model + @def ocp begin + # constants + beta = 1 + c = 2 + gamma = 1 + halfperiod = 5 + Ks = 0.05 + mu2m = 0.1 + mubar = 1 + r = 0.005 + T = 10 * N + + # ocp + t ∈ [0, T], time + x ∈ R³, state + u ∈ R, control + y = x[1] + s = x[2] + b = x[3] + mu2 = mu2m * s(t) / (s(t) + 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 + + # dynamics + ẋ(t) == dynamics(t, x(t), u(t)) + + # objective + ∫(b(t) / (beta + c)) → max + end + + # dynamics + function dynamics(t, x, u) + y, s, b = x + pi = 3.141592653589793 + days = t / (halfperiod * 2) + tau = (days - floor(days)) * 2 * pi + light = max(0, sin(tau))^2 + mu = light * mubar + mu2 = mu2m * s / (s + Ks) + return [mu * y / (1 + y) - (r + u) * y, + -mu2 * b + u * beta * (gamma * y - s), + (mu2 - u * beta) * b] + end + + # Initial guess + init = (state = [50, 50, 50], control = 0.5) + + # NLPModel + nlp = direct_transcription(ocp ,init = init, grid_size = nh)[2] + + return nlp + +end \ No newline at end of file diff --git a/ext/OptimalControlModels/insurance.jl b/ext/OptimalControlModels/insurance.jl new file mode 100644 index 0000000..4c8c6d8 --- /dev/null +++ b/ext/OptimalControlModels/insurance.jl @@ -0,0 +1,81 @@ +""" +The Insurance Problem: + The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.insurance(::OptimalControlBackend; nh::Int=100) + # constants + gamma = 0.2 + lambda = 0.25 + h0 = 1.5 + w = 1 + s = 10 + k = 0 + sigma = 0 + alpha = 4 + tf = 10 + + # Model + @def ocp 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 + + # dynamics + ẋ(t) == [(1 - gamma * t * vprime / dUdR(t)) * h(t), h(t), (1 + sigma) * I(t) * fx] + + # objective + ∫(U(t) * fx) → max + end + + # Initial guess + init = () + + # NLPModel + nlp = direct_transcription(ocp ,init = init, grid_size = nh)[2] + + return nlp + +end \ No newline at end of file diff --git a/ext/OptimalControlModels/jackson.jl b/ext/OptimalControlModels/jackson.jl new file mode 100644 index 0000000..a4274a1 --- /dev/null +++ b/ext/OptimalControlModels/jackson.jl @@ -0,0 +1,35 @@ +""" +The Jackson Problem: + The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.jackson(::OptimalControlBackend; nh::Int=100, N::Int=30) + @def ocp 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 + + # Initial guess + init = () + + # NLPModel + nlp = direct_transcription(ocp ,init = init, grid_size = nh)[2] + + return nlp +end \ No newline at end of file diff --git a/ext/OptimalControlModels/robbins.jl b/ext/OptimalControlModels/robbins.jl new file mode 100644 index 0000000..78f3957 --- /dev/null +++ b/ext/OptimalControlModels/robbins.jl @@ -0,0 +1,30 @@ +""" +The Robbins Problem: + The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) +""" +function OptimalControlProblems.robbins(::OptimalControlBackend; nh::Int=100, N::Int=30) + @def ocp 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 + + # Initial guess + init = () + + # NLPModel + nlp = direct_transcription(ocp ,init = init, grid_size = nh)[2] + + return nlp + +end \ No newline at end of file diff --git a/src/OptimalControlProblems.jl b/src/OptimalControlProblems.jl index 07c89f1..8296019 100644 --- a/src/OptimalControlProblems.jl +++ b/src/OptimalControlProblems.jl @@ -1,7 +1,6 @@ module OptimalControlProblems using CTBase -using DataFrames abstract type AbstractModelBackend end struct JuMPBackend <: AbstractModelBackend end @@ -58,10 +57,18 @@ The following keys are valid: - `ncon::Int`: number of general constraints - `minimize::Bool`: true if optimize == minimize """ -const metadata = DataFrame(infos .=> [Array{T}(undef, number_of_problems) for T in types]) +const metadata = Dict() -for data in infos, i in 1:number_of_problems - metadata[!, data][i] = eval(Meta.parse("$(split(files[i], ".")[1])_meta"))[data] +for i in 1:number_of_problems + file_key = Symbol(split(files[i], ".")[1]) + metadata[file_key] = Dict() + for (data, T) in zip(infos, types) + value = eval(Meta.parse("$(file_key)_meta"))[data] + if !(value isa T) + error("Type mismatch: Expected $(T) for $(data), but got $(typeof(value))") + end + metadata[file_key][data] = value + end end export JuMPBackend, OptimalControlBackend