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.

1 Like

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.