Optimization on Vector-ODEs for parameter estimation

So I’m trying to estimate Parameters from given Data, which I generated by solving the ODE problem with given parameters first and then adding a random stray around the result. Using the build_loss_objective I generated the cost function which is then inserted into optimization. From here, the estimated parameter values can be gathered by result.minimizer .

My first working example for regular ODE looked like that:

using DifferentialEquations, Random, Plots
using DiffEqParamEstim, Optim

\###############################################

# ODE Model

###############################################


function fun(G,p,t)
 k, β  = p
 #ODE
 dG_dt = k - β\*G
 return dG_dt
end

\#Parameter values
p = [0.8, 0.2] #"true" Parameter values
G_ini = 0.0
tspan = (0.0,10.0)
\#ODE Problem
prob = ODEProblem(fun,G_ini,tspan,p)
\#solving ODE Problem
sol = solve(prob,Tsit5())

\###############################################

# Synthetic Data

\###############################################

dataset = [(t,sol(t)+0.2randn()) for t in 0:0.01:10]

\#########
\#Cost function and Optimization
\############################

\#Gathering

dataTime = [d[1] for d in dataset] #time points
dataValues = [d[2] for d in dataset] #Data

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(dataTime,dataValues),maxiters=10000,verbose=false)

initialGuess = ones(2) #guesssed parameter values
result = optimize(cost_function, initialGuess, BFGS())

with that example I can get estimated values for the 2 parameters k and β

But now I want to extend the parameter estimation and for that I’d like to estimate parameters on a system of Vector-ODEs. Going from my last simpler example I got something like:

using DifferentialEquations,DiffEqParamEstim,RecursiveArrayTools
using Optim,Plots

\###############################################

# ODE Model

 \###############################################

n=10
function f(du,u,p,t)
 du[1,1] = dc_n = 4*p[1]*u[1,1] -u[1,1]*u[1,2]  
 du[1,2] = dc_d = p[2]*u[1,2]^2 + u[1,1]*u[1,2]

 du[n,1] = dc_n = 3*p[1]*u[n,1] -u[n,1]*u[n,2]                     
 du[n,2] = dc_d = -2*p[2]*u[n,2]^2 + u[n,1]*u[n,2]


 for i in 2:(n-1)
    du[i,1] = dc_n = p[1]*u[i-1,1] -u[i,1]*u[i+1,2]               
    du[i,2] = dc_d = -p[2]*u[i-1,2]^2 + u[i,1]*u[i+1,2]
 end

end

u0 = hcat(ones(n), ones(n))
tspan = (0.0,10.0)
p = [1.5,2.0,1.2]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

\###############################################

# Synthetic Data

\###############################################
dataset = [(t,sol(t) .+0.2randn()) for t in 0:0.01:10]

dataTime = [d[1] for d in dataset] #Auffassen Zeitpunkte
dataValues = [d[2] for d in dataset] #Auffassen Exp. Data
\#until here everything works

cost_function =  build_loss_objective(prob,Tsit5(),L2Loss(dataTime,dataValues),maxiters=10000,verbose=false)
initialGuess = ones(3) #alle Parameter auf eins gesetzt (Ratewerte)
result = optimize(cost_function, initialGuess, BFGS())

when running the “optimze” function, it tells me in the terminal that there has been a failure and no parameters have been estimated. Where exactly would I need to change something to make this code work?

When running your code I see:

julia> sol = solve(prob, Tsit5());
┌ Warning: dt(1.7763568394002505e-15) <= dtmin(1.7763568394002505e-15) at
│ t=0.35552611924380084, and step error estimate = 0.19101346055879606. Aborting.
│ There is either an error in your model specification or the true solution is
│ unstable.
└ @ SciMLBase /home/user/.julia/packages/SciMLBase/kTUaf/src/integrator_interface.jl:599

I suspect this is what is causing the optimization to report failure.

Thank your for the Reply, what could be meant with “model specification”?

It means there is no timestep small enough to successfully integrate f(du,u,p,t) as implemented beyond t=0.35552611924380084. Probably a typo, missing term, swapped sign or something like that.