I am new to Julia, trying to fit the parameters of an ODE model using JuMP and Ipopt. The code works as intended in a parameter range of p_true/fc <= p <= fc*p_true if fc = 2, 4, or 8. However, fc = 3 or 10 it prints Defining ODEs ... and then gets stuck (interestingly without any load on the CPU). I also tried lower discretisation steps of the ODEs, but experienced a very sharp threshold when the code suddenly became prohibitively slow. I copied in the code below, and the data to fit can be found in this csv file. I am really puzzled and would appreciate your advice.
Thanks,
Paul
using JuMP
using Ipopt
using CSV
fc = 3
# Data
data_file_path = joinpath(pwd(), "G2M.csv")
df = CSV.read(data_file_path)
t = Vector(df[!, :t]) # set of simulation times.
# create JuMP model object
m = Model(with_optimizer(Ipopt.Optimizer))
# set_parameter(m, "hessian_approximation", "limited-memory") # This setting makes the code very slow.
# Model parameters
println("Defining parameters ...")
@variable(m, 0.1/fc <= kPhEnsa <= 0.1*fc)
@variable(m, 0.05/fc <= kDpEnsa <= 0.05*fc)
@variable(m, 1/fc <= kPhGw <= 1*fc)
@variable(m, 0.25/fc <= kDpGw1 <= 0.25*fc)
@variable(m, 10/fc <= kDpGw2 <= 10*fc)
@variable(m, 0.01/fc <= kWee1 <= 0.01*fc)
@variable(m, 1/fc <= kWee2 <= 1*fc)
@variable(m, 1/fc <= kPhWee <= 1*fc)
@variable(m, 10/fc <= kDpWee <= 10*fc)
@variable(m, 0.1/fc <= kCdc25_1 <= 0.1*fc)
@variable(m, 1/fc <= kCdc25_2 <= 1*fc)
@variable(m, 1/fc <= kPhCdc25 <= 1*fc)
@variable(m, 10/fc <= kDpCdc25 <= 10*fc)
@variable(m, 0.0068/fc <= kDipEB55 <= 0.0068*fc)
@variable(m, 57/fc <= kAspEB55 <= 57*fc)
@variable(m, 0 <= tCb <= 2)
@variable(m, 0 <= tEnsa <= 2)
@variable(m, 0 <= tB55 <= 2)
@variable(m, 0 <= tGw <= 2)
# Model states
println("Defining states ...")
@variable(m, 0 <= pEnsa[k in 1:length(t)] <= 1)
@variable(m, 0 <= pGw[k in 1:length(t)] <= 1)
@variable(m, 0 <= pEB55[k in 1:length(t)] <= 1)
@variable(m, 0 <= Cb[k in 1:length(t)] <= 1)
# Model ODEs
println("Defining ODEs ...")
@NLconstraint(m, [k in 1:length(t)-1],
pEnsa[k+1] == pEnsa[k] + ( kPhEnsa*pGw[k+1]*(tEnsa-pEnsa[k+1]) - kDpEnsa*pEB55[k+1] ) * ( t[k+1] - t[k] ) )
@NLconstraint(m, [k in 1:length(t)-1],
Cb[k+1] == Cb[k] + ( (kCdc25_1 + (kCdc25_2-kCdc25_1)*(kPhCdc25*Cb[k+1])/(kPhCdc25*Cb[k+1] + kDpCdc25*(tB55 - pEB55[k+1]))) * (tCb - Cb[k+1]) - (kWee1+(kWee2-kWee1)*kDpWee*(tB55 - pEB55[k+1])/(kDpWee*(tB55 - pEB55[k+1]) + kPhWee*Cb[k+1]))*Cb[k+1] ) * ( t[k+1] - t[k]) )
@NLconstraint(m, [k in 1:length(t)-1],
pGw[k+1] == pGw[k] + ( kPhGw*Cb[k+1]*(tGw-pGw[k+1]) - (kDpGw1+kDpGw2*(tB55 - pEB55[k+1]))*pGw[k+1] ) * ( t[k+1] - t[k]) )
@NLconstraint(m, [k in 1:length(t)-1],
pEB55[k+1] == pEB55[k] + ( kAspEB55*(tB55 - pEB55[k+1])*(pEnsa[k+1] - pEB55[k+1]) - (kDipEB55 + kDpEnsa)*pEB55[k+1] ) * ( t[k+1] - t[k]) )
# define objective
println("Defining objective ...")
@NLobjective(m, Min,
sum((pEnsa[k]-df[k, :pEnsa])^2 for k in 1:length(t)) +
sum((Cb[k]-df[k, :Cb])^2 for k in 1:length(t)) +
sum((pGw[k]-df[k, :pGw])^2 for k in 1:length(t)) # +
# sum((pEB55[k]-df[k, :pEB55])^2 for k in 1:length(t)) # This line causes failure
)
println("Optimizing...")
optimize!(m)