11 dev add ctdirectjl test problems

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"

ext/JuMPModels/beam.jl

@@ -0,0 +1,39 @@
+"""
+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
+
+    # 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,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

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

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

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

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/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/robbins.jl

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

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

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