Hi @odow, apologies for the late reply, I did check the dictionaries and as far I understood it they did match with their corresponding values (although I might be wrong).
I did an experiment warm starting only the primals and got the same results.
Here is a reproducible example:
using JuMP, Ipopt, Plots, MINPACK, LaTeXStrings, Printf
# Model Parameters
# Constants
c1 = .0025 # k_on
c2 = 2500 # TF_tot
c3 = 10 # k_off
c4 = 20000 # k_m
c5 = .05 # k_m,deg,b
c6 = 1500 # K_m
c7 = .8 # k_rib,max
c8 = .1 # mu_max
c9 = 2 # k_u,trans
c10 = .3 # delta
c11 = 1250 # K_CU
c12 = 4000 # beta
c13 = .001 # gamma
c14 = 37 # alpha
c15 = 10 # Y_X/S
c16 = 250 # U_max
c17 = 5 # K_S
c18 = 1000000 # K_IU
c19 = 2500 # K_CI
c20 = .2 # delta_1
c21 = 100 # X(0)
c22 = 20 # S(0)
# Initial State
x1_0 = 1 # Transcriptional Reporter
x2_0 = 1 # mRNA
x3_0 = 1 # U, unfolded proteins
x4_0 = 1 # H, Hac1p transcription factor
x5_0 = 1 # C, unfolded proteins buffer?
x6_0 = 0 # FA, folded amylase
x7_0 = c22 # S, substrate
x8_0 = 1 # int(FA)
# Control Bounds
u_min = 0 # This comes from 1 where the light is bounded between 0-100%
u_max = 100 # This comes from 1 where the light is bounded between 0-100%
# Initial State
xi_0 = [x1_0, x2_0, x3_0, x4_0, x5_0, x6_0, x7_0, x8_0];
# Define Runge-Kutta Butcher tableau data structure
struct rk_method
name::Symbol
s::Integer
a::Matrix{<:Real}
b::Vector{<:Real}
c::Vector{<:Real}
end
# Define Gauss-Legendre method (of order 2s) with s=2
# https://en.wikipedia.org/wiki/Gauss-Legendre_method
rk = rk_method(:gauss2, 2,
[0.25 (0.25 - sqrt(3)/6); (0.25 + sqrt(3)/6) 0.25],
[0.5, 0.5],
[(0.5 - sqrt(3)/6), (0.5 + sqrt(3)/6)]
)
# Initialize JuMP model with Ipopt solver backend
docp = JuMP.Model(Ipopt.Optimizer)
set_optimizer_attribute(docp, "print_level", 5)
set_optimizer_attribute(docp, "tol", 1e-16)
set_optimizer_attribute(docp, "constr_viol_tol", 1e-16)
set_optimizer_attribute(docp, "max_iter", 1500)
set_optimizer_attribute(docp, "mu_strategy", "adaptive")
# Discretization parameters
N = 1000 # time steps
t0 = 0; tf = 20
Δt = (tf - t0) / N
# System definition
n = 8 # dim(x)
# Declare constrained variables
x_lower = [0, 0, 0, 0, 0, 0, 0, 0] # lower bound for x
@variables(docp, begin
x[1:N, i=1:n] ≥ x_lower[i] # x
0 ≤ u[1:N] ≤ 100 # u
k[1:rk.s, 1:N, 1:n] # k (for Runge-Kutta)
end)
# Objective function
@objective(docp, Max, x[end, 6])
# Initial condition
@constraint(docp, initial, x[1,:] == xi_0[:])
# Autonomous dynamics function
function f(x, u)
mu = (1 - x[3]/(x[3] + c16)) * c8 * x[7]/(x[7] + c17)
psi = x[3]/(x[3] + c18 + (c11 + c18 * x[5])/(c19 * (c11 + x[3])))
dx1 = c1 * u * (c2 - x[1]) - c3 * x[1]
dx2 = c4 * x[1]/(c6 + x[1]) - (c5 + mu) * x[2]
dx3 = mu * c9 * x[2] *c7/c8 - c10 * (x[5] * x[3]/(x[3] + c11)) - mu * x[3]
dx4 = mu * c12 * psi * c7/c8 - (c5 + mu) * x[4]
dx5 = mu * (c13 + c14 * x[4]) * c7/c8 - mu * x[5]
dx6 = c20 * (x[5] * x[3]/(x[3] + c11)) * (c21 + c15 * (c22 - x[7]))
dx7 = - mu * (c21 + c15 * (c22 - x[7]))/c15
dx8 = x[6]
F = [ dx1,
dx2,
dx3,
dx4,
dx5,
dx6,
dx7,
dx8,
]
return F
end
# Runge-Kutta methods for autonomous systems (as a nonlinear constraint)
# x[i+1] = x[i] + Δt Σ_j b[j]k[j,i]
# k[j,i] = f( x[i] + Δt Σ_s A[j,s]k[s,i] )
@constraints(docp, begin
rk_nodes[j=1:rk.s, i=1:N], k[j,i,:] == f(x[i,:] + Δt*sum(rk.a[j,s]*k[s,i,:] for s in 1:rk.s), u[i])
rk_scheme[i = 1:N-1], x[i+1,:] == x[i,:] + Δt*sum(rk.b[j] * k[j,i,:] for j in 1:rk.s)
end);
@timev begin
# Solve DOCP as nonlinear optimization problem
optimize!(docp)
end;
####### SECOND MODEL ########
# Initialize JuMP model with Ipopt solver backend
docp_II = JuMP.Model(Ipopt.Optimizer)
set_optimizer_attribute(docp_II, "print_level", 5)
set_optimizer_attribute(docp_II, "tol", 1e-16)
set_optimizer_attribute(docp_II, "constr_viol_tol", 1e-16)
set_optimizer_attribute(docp_II, "max_iter", 1500)
set_optimizer_attribute(docp_II, "mu_strategy", "adaptive")
# Discretization parameters
N = 1000 + 50 # time steps
t0 = 0; tf = 21
Δt = (tf - t0) / N
# System definition
n = 8 # dim(x)
# Declare constrained variables
x_lower = [0, 0, 0, 0, 0, 0, 0, 0] # lower bound for x
@variables(docp_II, begin
x[1:N, i=1:n] ≥ x_lower[i] # x
0 ≤ u[1:N] ≤ 100 # u
k[1:rk.s, 1:N, 1:n] # k (for Runge-Kutta)
end)
# Objective function
@objective(docp_II, Max, x[end, 6])
# Initial condition
@constraint(docp_II, initial, x[1,:] == xi_0[:])
# Runge-Kutta methods for autonomous systems (as a nonlinear constraint)
# x[i+1] = x[i] + Δt Σ_j b[j]k[j,i]
# k[j,i] = f( x[i] + Δt Σ_s A[j,s]k[s,i] )
@constraints(docp_II, begin
rk_nodes[j=1:rk.s, i=1:N], k[j,i,:] == f(x[i,:] + Δt*sum(rk.a[j,s]*k[s,i,:] for s in 1:rk.s), u[i])
rk_scheme[i = 1:N-1], x[i+1,:] == x[i,:] + Δt*sum(rk.b[j] * k[j,i,:] for j in 1:rk.s)
end);
function set_initial_guess!(docp::GenericModel, u::Union{AbstractMatrix, AbstractVector}, x::AbstractMatrix)
set_optimizer_attribute(docp, "warm_start_init_point", "yes")
new_N = size(docp.obj_dict[:u], 1)
old_N = size(u, 1)
x_docp = docp.obj_dict[:x]
u_docp = docp.obj_dict[:u]
for i=1:old_N
set_start_value.(x_docp[i, :], x[i])
set_start_value.(u_docp[i, :], u[i])
end
end
x = value.(docp.obj_dict[:x])
u = value.(docp.obj_dict[:u])
set_initial_guess!(docp_II, u, x)
@timev begin
# In here the objective function should be closer to my previous values right?
optimize!(docp_II)
end;
