11 dev add ctdirectjl test problems

Project.toml

@@ -5,6 +5,7 @@ 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,4 +19,5 @@ OptimalControlModels = "CTDirect"
 CTBase = "0.13"
 CTDirect = "0.12"
 JuMP = "1"
+DataFrames = "1"
 julia = "1.10"

ext/JuMPModels/beam.jl

@@ -0,0 +1,39 @@
+"""
+The Beam Problem:
+    The problem is formulated as a JuMP modeln and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop)
+"""
+function OptimalControlProblems.beam(::JuMPBackend;nh::Int=100)
+    # Parameters
+    tf = 1
+
+    # 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
+    @expressions(model, begin
+        step, tf/nh
+    end)    
+    @constraints(model, begin
+        con_x1[t=1:nh], x1[t] == x1[t-1] + 0.5 * step[t] * (x2[t] + x2[t-1])
+        con_x2[t=1:nh], x2[t] == x2[t-1] + 0.5 * step[t] * (u[t] + u[t-1])
+    end)
+
+
+    @objective(model, Min, sum(u[t]^2 for t in 0:nh))
+
+    return model
+end

ext/JuMPModels/bioreactor.jl

@@ -0,0 +1,58 @@
+"""
+The Bioreactor Problem:
+    The problem is formulated as a JuMP modeln 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
+    pi = 3.141592653589793
+
+    # 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

ext/JuMPModels/insurance.jl

@@ -0,0 +1,70 @@
+"""
+The Insurance Problem:
+    The problem is formulated as a JuMP modeln 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

ext/JuMPModels/jackson.jl

@@ -0,0 +1,47 @@
+"""
+The Jackson Problem:
+    The problem is formulated as a JuMP modeln 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

ext/JuMPModels/robbins.jl

@@ -0,0 +1,45 @@
+"""
+The Robbins Problem:
+    The problem is formulated as a JuMP modeln 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

ext/JuMPModels/vanderpol.jl

@@ -1,6 +1,6 @@
 """
 The Van der Pol Problem:
-    The problem is formulated as an OptimalControl modeln and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop)
+    The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop)
 """
 function OptimalControlProblems.vanderpol(::JuMPBackend; nh::Int=100)
     # parameters

ext/MetaData/beam.jl

@@ -0,0 +1 @@
+beam_meta = Dict( :name => "beam", :nh => 100, :nvar => 404, :ncon => 305, :minimize => true,)

ext/MetaData/bioreactor.jl

@@ -0,0 +1 @@
+bioreactor_meta = Dict(:name => "bioreactor", :nh => 100, :nvar => 505, :ncon => 404, :minimize => false)

ext/MetaData/cart_pendulum.jl

@@ -1,3 +1 @@
-cart_pendulum_meta = Dict(
-    :name => "cart_pendulum", :nvar => 507, :ncon => 405, :minimize => true
-)
+cart_pendulum_meta = Dict( :name => "cart_pendulum", :nh => 100, :nvar => 507, :ncon => 405, :minimize => true)

ext/MetaData/chain.jl

@@ -1 +1 @@
-chain_meta = Dict(:name => "chain", :nvar => 404, :ncon => 305, :minimize => true)
+chain_meta = Dict(:name => "chain", :nh => 100, :nvar => 404, :ncon => 305, :minimize => true)

ext/MetaData/dielectrophoretic_particle.jl

@@ -1,3 +1 @@
-dielectrophoretic_particle_meta = Dict(
-    :name => "dielectrophoretic_particle", :nvar => 304, :ncon => 203, :minimize => true
-)
+dielectrophoretic_particle_meta = Dict(:name => "dielectrophoretic_particle", :nh => 300, :nvar => 304, :ncon => 203, :minimize => true)

ext/MetaData/double_oscillator.jl

@@ -1,3 +1 @@
-double_oscillator_meta = Dict(
-    :name => "double_oscillator", :nvar => nothing, :ncon => nothing, :minimize => true
-)
+double_oscillator_meta = Dict(:name => "double_oscillator", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true)

ext/MetaData/ducted_fan.jl

@@ -1,3 +1 @@
-ducted_fan_meta = Dict(
-    :name => "ducted_fan", :nvar => nothing, :ncon => nothing, :minimize => true
-)
+ducted_fan_meta = Dict(:name => "ducted_fan", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true)

ext/MetaData/electrical_vehicle.jl

@@ -1,3 +1 @@
-electrical_vehicle_meta = Dict(
-    :name => "electrical_vehicle", :nvar => nothing, :ncon => nothing, :minimize => true
-)
+electrical_vehicle_meta = Dict(:name => "electrical_vehicle", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true)

ext/MetaData/glider.jl

@@ -1 +1 @@
-glider_meta = Dict(:name => "glider", :nvar => 506, :ncon => 407, :minimize => false)
+glider_meta = Dict(:name => "glider", :nh => 100, :nvar => 506, :ncon => 407, :minimize => false)

ext/MetaData/insurance.jl

@@ -0,0 +1 @@
+insurance_meta = Dict(:name => "insurance", :nh => 100, :nvar => 910, :ncon => 809, :minimize => false)

ext/MetaData/jackson.jl

@@ -0,0 +1 @@
+jackson_meta = Dict(:name => "jackson", :nh => 100, :nvar => 404, :ncon => 303, :minimize => false)

ext/MetaData/moonlander.jl

@@ -1 +1 @@
-moonlander_meta = Dict(:name => "moonlander", :nvar => 610, :ncon => 809, :minimize => true)
+moonlander_meta = Dict(:name => "moonlander", :nh => 100, :nvar => 610, :ncon => 809, :minimize => true)

ext/MetaData/quadrotor.jl

@@ -1,3 +1 @@
-quadrotor_meta = Dict(
-    :name => "quadrotor", :nvar => nothing, :ncon => nothing, :minimize => true
-)
+quadrotor_meta = Dict(:name => "quadrotor", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true)

ext/MetaData/robbins.jl

@@ -0,0 +1 @@
+robbins_meta = Dict(:name => "robbins", :nh => 100, :nvar => 505, :ncon => 407, :minimize => true)

ext/MetaData/robot.jl

@@ -1 +1 @@
-robot_meta = Dict(:name => "robot", :nvar => 910, :ncon => 612, :minimize => true)
+robot_meta = Dict(:name => "robot", :nh => 100, :nvar => 910, :ncon => 612, :minimize => true)

ext/MetaData/rocket.jl

@@ -1 +1 @@
-rocket_meta = Dict(:name => "rocket", :nvar => 405, :ncon => 304, :minimize => false)
+rocket_meta = Dict(:name => "rocket", :nh => 100, :nvar => 405, :ncon => 304, :minimize => false)

ext/MetaData/space_shuttle.jl

@@ -1,3 +1 @@
-space_shuttle_meta = Dict(
-    :name => "space_shuttle", :nvar => nothing, :ncon => nothing, :minimize => false
-)
+space_shuttle_meta = Dict(:name => "space_shuttle", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => false)

ext/MetaData/steering.jl

@@ -1 +1 @@
-steering_meta = Dict(:name => "steering", :nvar => 506, :ncon => 408, :minimize => true)
+steering_meta = Dict(:name => "steering", :nh => 100, :nvar => 506, :ncon => 408, :minimize => true)

ext/MetaData/truck_trailer.jl

@@ -1 +1 @@
-truck_meta = Dict(:name => "truck", :nvar => nothing, :ncon => nothing, :minimize => true)
+truck_trailer_meta = Dict(:name => "truck", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true)

ext/MetaData/vanderpol.jl

@@ -1 +1 @@
-vanderpol_meta = Dict(:name => "vanderpol", :nvar => 404, :ncon => 303, :minimize => true)
+vanderpol_meta = Dict(:name => "vanderpol", :nh => 100, :nvar => 404, :ncon => 303, :minimize => true)