# 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.