test new inits

ext/JuMPModels/bioreactor.jl

@@ -79,20 +79,20 @@ function OptimalControlProblems.bioreactor(
     @variables(
         model,
         begin
-            y[0:N] ≥ y_l, (start = 50)
-            s[0:N] ≥ s_l, (start = 50)
-            b[0:N] ≥ b_l, (start = 50)
-            u_l ≤ u[0:N] ≤ u_u, (start = 0.5)
+            y[0:N] ≥ 0, (start = 0.15)     
+            s[0:N] ≥ 0, (start = 2.75)     
+            b[0:N] ≥ 0.001, (start = 1.75) 
+            0 ≤ u[0:N] ≤ 1, (start = 0.5)
         end
     )
 
     # boundary constraints
     @constraints(
         model,
         begin
-            y_t0_l ≤ y[0] ≤ y_t0_u
-            s_t0_l ≤ s[0] ≤ s_t0_u
-            b_t0_l ≤ b[0] ≤ b_t0_u
+            y[0] == 0.05  
+            s[0] == 0.5   
+            b[0] == 0.5   
         end
     )
 

ext/JuMPModels/glider.jl

@@ -76,12 +76,12 @@ function OptimalControlProblems.glider(
     @variables(
         model,
         begin
-            tf ≥ tf_l, (start = 1)
-            x[k = 0:N] ≥ x_l, (start = x_t0 + vx_t0 * k / N)
+            tf ≥ tf_l, (start = 100.0)  ####
+            x[k = 0:N] ≥ x_l, (start = x_t0 + (k / N) * 1248.0) ####
             y[k = 0:N], (start = y_t0 + (k / N) * (y_tf - y_t0))
             vx[k = 0:N] ≥ vx_l, (start = vx_t0)
             vy[k = 0:N], (start = vy_t0)
-            cL_min ≤ cL[k = 0:N] ≤ cL_max, (start = cL_max / 2)
+            cL_min ≤ cL[k = 0:N] ≤ cL_max, (start = 1.0) ####
         end
     )
 

ext/JuMPModels/robot.jl

@@ -92,22 +92,53 @@ function OptimalControlProblems.robot(
     @variables(
         model,
         begin
-            ρ_l ≤ ρ[k = 0:N] ≤ ρ_u, (start = ρ_t0)
-            θ_l ≤ θ[k = 0:N] ≤ θ_u, (start = 2π/3 * (k / N)^2)
-            ϕ_l ≤ ϕ[k = 0:N] ≤ ϕ_u, (start = ϕ_t0)
+            ρ_l ≤ ρ[k = 0:N] ≤ ρ_u
+            θ_l ≤ θ[k = 0:N] ≤ θ_u
+            ϕ_l ≤ ϕ[k = 0:N] ≤ ϕ_u
 
-            dρ[k = 0:N], (start = 0)
-            dθ[k = 0:N], (start = 4π/3 * (k / N))
-            dϕ[k = 0:N], (start = 0)
+            dρ[k = 0:N]
+            dθ[k = 0:N]
+            dϕ[k = 0:N]
 
-            uρ_l ≤ uρ[0:N] ≤ uρ_u, (start = 0)
-            uθ_l ≤ uθ[0:N] ≤ uθ_u, (start = 0)
-            uϕ_l ≤ uϕ[0:N] ≤ uϕ_u, (start = 0)
+            uρ_l ≤ uρ[0:N] ≤ uρ_u
+            uθ_l ≤ uθ[0:N] ≤ uθ_u
+            uϕ_l ≤ uϕ[0:N] ≤ uϕ_u
 
-            tf ≥ tf_l, (start = 1)
+            tf ≥ tf_l
         end
     )
 
+    #INITIAL GUESS
+    set_start_value(tf, 9.1)
+    for k in 0:grid_size
+        # Coefficient d'interpolation entre 0 et 1
+        alpha = k / grid_size
+
+        # Interpolation linéaire entre t0 et tf
+        ρ_val = ρ_t0 + alpha * (ρ_tf - ρ_t0)
+        θ_val = θ_t0 + alpha * (θ_tf - θ_t0)
+        ϕ_val = ϕ_t0 + alpha * (ϕ_tf - ϕ_t0)
+
+        # SECURITÉ ANTI-CRASH (Division par zéro)
+        if abs(ϕ_val) < 1e-6
+            ϕ_val = 1e-6
+        end
+
+        # Application des valeurs aux variables JuMP
+        set_start_value(ρ[k], ρ_val)
+        set_start_value(θ[k], θ_val)
+        set_start_value(ϕ[k], ϕ_val)
+
+        # Vitesses et contrôles à 0 par défaut
+        set_start_value(dρ[k], 0.0)
+        set_start_value(dθ[k], 0.0)
+        set_start_value(dϕ[k], 0.0)
+
+        set_start_value(uρ[k], 0.0)
+        set_start_value(uθ[k], 0.0)
+        set_start_value(uϕ[k], 0.0)
+    end
+
     # Boundary condition
     @constraints(
         model,
@@ -137,8 +168,7 @@ function OptimalControlProblems.robot(
             #
             Δt, (tf - t0) / N
 
-            #
-            I_θ[i = 0:N], ((L - ρ[i])^3 + ρ[i]^3) * (sin(ϕ[i]))^2
+            I_θ[i = 0:N], ((L - ρ[i])^3 + ρ[i]^3) * ((sin(ϕ[i]))^2 + 1e-9)
             I_ϕ[i = 0:N], (L - ρ[i])^3 + ρ[i]^3
 
             #

ext/JuMPModels/space_shuttle.jl

@@ -184,7 +184,7 @@ function OptimalControlProblems.space_shuttle(
     set_start_value.(model[:ψ], vec(initial_guess[:, 6]))
     set_start_value.(model[:α], vec(initial_guess[:, 7]))
     set_start_value.(model[:β], vec(initial_guess[:, 8]))
-    set_start_value.(model[:tf], (tf_l+tf_u)/2)
+    set_start_value.(model[:tf], 2000.0) #######
 
     # Functions to restore `h` and `v` to their true scale
     @expression(model, h[j = 0:N], scaled_h[j] * scaling_h)

ext/JuMPModels/steering.jl

@@ -51,7 +51,7 @@ function OptimalControlProblems.steering(
     x₄_tf = params[:x₄_tf]
 
     #
-    tf_start = 1
+    tf_start = 0.6
     function gen_x0(k, i)
         if i == 1 || i == 4
             return 0.0

ext/OptimalControlModels/bioreactor.jl

@@ -27,82 +27,68 @@ julia> docp = OptimalControlProblems.bioreactor(OptimalControlBackend(); N=100);
 
 - BOCOP repository: https://github.com/control-toolbox/bocop/tree/main/bocop
 """
+
 function OptimalControlProblems.bioreactor(
     ::OptimalControlBackend,
     description::Symbol...;
     grid_size::Int=grid_size_data(:bioreactor),
     parameters::Union{Nothing,NamedTuple}=nothing,
     kwargs...,
 )
-
-    # parameters
+    # --- 1. Parameters ---
     params = parameters_data(:bioreactor, parameters)
-    t0 = params[:t0]
-    tf = params[:tf]
-    β = params[:β]
-    c = params[:c]
-    γ = params[:γ]
+    t0, tf = params[:t0], params[:tf]
+    β, c, γ = params[:β], params[:c], params[:γ]
     halfperiod = params[:halfperiod]
-    Ks = params[:Ks]
-    μ2m = params[:μ2m]
-    μbar = params[:μbar]
-    r = params[:r]
-    y_l = params[:y_l]
-    s_l = params[:s_l]
-    b_l = params[:b_l]
-    u_l = params[:u_l]
-    u_u = params[:u_u]
-    y_t0_l = params[:y_t0_l]
-    y_t0_u = params[:y_t0_u]
-    s_t0_l = params[:s_t0_l]
-    s_t0_u = params[:s_t0_u]
-    b_t0_l = params[:b_t0_l]
-    b_t0_u = params[:b_t0_u]
-
-    # Model
+    Ks, μ2m, μbar, r = params[:Ks], params[:μ2m], params[:μbar], params[:r]
+    y_l, s_l, b_l = params[:y_l], params[:s_l], params[:b_l]
+    u_l, u_u = params[:u_l], params[:u_u]
+    y_t0_l, y_t0_u = params[:y_t0_l], params[:y_t0_u]
+    s_t0_l, s_t0_u = params[:s_t0_l], params[:s_t0_u]
+    b_t0_l, b_t0_u = params[:b_t0_l], params[:b_t0_u]
+
+    # --- 2. Model ---
     ocp = @def begin
         t ∈ [t0, tf], time
-        x = (y, s, b) ∈ R³, state
+        x = (y, s, b, k) ∈ R⁴, state
         u ∈ R, control
 
-        x(t) ≥ [y_l, s_l, b_l]
-        u_l ≤ u(t) ≤ u_u
-        [y_t0_l, s_t0_l, b_t0_l] ≤ x(t0) ≤ [y_t0_u, s_t0_u, b_t0_u]
-
-        μ = light(t, halfperiod) * μbar
-        μ2 = growth(s(t), μ2m, Ks)
-
-        ẋ(t) == [
-            μ * y(t) / (1 + y(t)) - (r + u(t)) * y(t),
-            -μ2 * b(t) + u(t) * β * (γ * y(t) - s(t)),
-            (μ2 - u(t) * β) * b(t),
-        ]
+
+        y(t) ≥ 0
+        s(t) ≥ 0
+        b(t) ≥ 1e-3       
+        k(t) ≥ t0        
+        k(t) ≤ tf         
 
-        -∫(μ2 * b(t) / (β + c)) → min
-    end
-
-    # METHANE PROBLEM
-    # μ2 according to growth model
-    # μ according to light model
-    # time scale is [0,10] for 24h (day then night)
-
-    # growth model MONOD
-    function growth(s, μ2m, Ks)
-        return μ2m * s / (s + Ks)
-    end
-
-    # light model: max^2 (0,sin) * μbar
-    # DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi]
-    function light(time, halfperiod)
-        days = time / (halfperiod * 2)
-        tau = (days - floor(days)) * 2π
-        return max(0, sin(tau))^2
+
+        u_l ≤ u(t) ≤ u_u
+
+
+        y_t0_l ≤ y(t0) ≤ y_t0_u
+        s_t0_l ≤ s(t0) ≤ s_t0_u
+        b_t0_l ≤ b(t0) ≤ b_t0_u
+        k(t0) == t0       
+
+
+
+        # dy/dt
+        ẋ[1](t) == (μbar * max(0, sin(k(t) * π / halfperiod))^2) * y(t) / (1 + y(t)) - (r + u(t)) * y(t)
+
+        # ds/dt
+        ẋ[2](t) == -(μ2m * s(t) / (s(t) + Ks)) * b(t) + u(t) * β * (γ * y(t) - s(t))
+
+        # db/dt
+        ẋ[3](t) == ((μ2m * s(t) / (s(t) + Ks)) - u(t) * β) * b(t)
+
+        ẋ[4](t) == 1
+
+        # Objectif
+        -∫((μ2m * s(t) / (s(t) + Ks)) * b(t) / (β + c)) → min
     end
 
-    # initial guess
-    init = (state=[50, 50, 50], control=0.5)
+    # --- 3. Transcription ---
+    init = (state=[0.15, 2.75, 1.75, t0], control=0.5)
 
-    # discretise the optimal control problem
     docp = direct_transcription(
         ocp,
         description...;
@@ -114,4 +100,4 @@ function OptimalControlProblems.bioreactor(
     )
 
     return docp
-end
+end

ext/OptimalControlModels/glider.jl

@@ -110,11 +110,15 @@ function OptimalControlProblems.glider(
     end
 
     # initial guess
-    tfinit = 1
-    xinit =
-        t -> [x_t0 + vx_t0 * t / tfinit, y_t0 + t / tfinit * (y_tf - y_t0), vx_t0, vy_t0]
-    uinit = cL_max / 2
-    init = (state=xinit, control=uinit, variable=tfinit)
+    tfinit = 100.0 
+    xinit = t -> [
+        x_t0 + (t / tfinit) * 1248.0,                         
+        y_t0 + (t / tfinit) * (y_tf - y_t0),         
+        vx_t0 + (t / tfinit) * (vx_tf - vx_t0),       
+        vy_t0 + (t / tfinit) * (vy_tf - vy_t0)        
+    ]
+    uinit = 1.0
+    init = (state=xinit, control=uinit, variable=tfinit) 
 
     # discretise the optimal control problem
     docp = direct_transcription(

ext/OptimalControlModels/robot.jl

@@ -105,30 +105,42 @@ function OptimalControlProblems.robot(
         dθ(tf) == dθ_tf, (dθ_tf)
         dϕ(tf) == dϕ_tf, (dϕ_tf)
 
-        # aliases
-        I_θ = ((L - ρ(t))^3 + ρ(t)^3) * sin(ϕ(t))^2
+
+        I_θ = ((L - ρ(t))^3 + ρ(t)^3) * (sin(ϕ(t))^2 + 1e-9)
         I_ϕ = (L - ρ(t))^3 + ρ(t)^3
 
         # dynamics  
-        ẋ(t) == [dρ(t), uρ(t) / L, dθ(t), 3 * uθ(t) / I_θ, dϕ(t), 3 * uϕ(t) / I_ϕ]
+        ẋ(t) == [dρ(t), uρ(t) / L, dθ(t), 3 * uθ(t) / I_θ, dϕ(t), 3 * uϕ(t) / I_ϕ]
 
         # objective
         tf → min
     end
 
     # initial guess
-    tf = 1
-    xinit =
-        t -> [
-            ρ_t0,
-            0,
-            2π/3 * ((t - t0) / (tf - t0))^2,
-            4π/3 * ((t - t0) / (tf - t0)),
-            ϕ_t0,
-            0,
+    tf_guess = 9.1
+    xinit = t -> begin
+        alpha = clamp((t - t0) / (tf_guess - t0), 0.0, 1.0)
+
+        ρ_val = ρ_t0 + alpha * (ρ_tf - ρ_t0)
+        θ_val = θ_t0 + alpha * (θ_tf - θ_t0)
+        ϕ_val = ϕ_t0 + alpha * (ϕ_tf - ϕ_t0)
+
+        if abs(ϕ_val) < 1e-6
+            ϕ_val = 1e-6
+        end
+
+        return [
+            ρ_val, 
+            0.0, 
+            θ_val, 
+            0.0, 
+            ϕ_val, 
+            0.0,
         ]
-    uinit = [0, 0, 0]
-    init = (state=xinit, control=uinit, variable=tf)
+    end
+
+    uinit = [0.0, 0.0, 0.0]
+    init = (state=xinit, control=uinit, variable=tf_guess)
 
     # discretise the optimal control problem
     docp = direct_transcription(

ext/OptimalControlModels/space_shuttle.jl

@@ -56,8 +56,8 @@ function OptimalControlProblems.space_shuttle(
     # 
     Δt_min = params[:Δt_min]
     Δt_max = params[:Δt_max]
-    tf_l = grid_size*Δt_min
-    tf_u = grid_size*Δt_max
+    tf_l = 1500.0  #
+    tf_u = 2500.0  #
 
     ## Initial conditions
     h_t0 = params[:h_t0]
@@ -175,7 +175,7 @@ function OptimalControlProblems.space_shuttle(
     # initial guess: linear interpolation for h, v, gamma (NB. t0 = 0), constant for the rest
     # variable time step seems to be initialized at 1 in jump
     # note that ipopt will project the initial guess inside the bounds anyway.
-    tf_init = (tf_l+tf_u)/2
+    tf_init = 2000
     x_init =
         t -> [
             h_t0 + (t - t0) / (tf_init - t0) * (h_tf - h_t0),

ext/OptimalControlModels/steering.jl

@@ -85,7 +85,7 @@ function OptimalControlProblems.steering(
         end
     end
     xinit = t -> [gen_x0(t, i) for i in 1:4]
-    init = (state=xinit, control=0, variable=1)
+    init = (state=xinit, control=0, variable=0.6)
 
     # discretise the optimal control problem
     docp = direct_transcription(