How to incorporate time dependent interpolation functions while solving ensemble parameter estimation problem?

New to Julia
1

Hi everyone,

Data interpolation functions/ objects are not part of the ODE function in my code and also not the part of the parameter list.

Example:

function V(du, u, p, t)
    du[1] =  Rate(t)
    du[2] = p[1]*u[1] + (Rate(t) / u[1]) *u[2]
end

p0 = [1.]
u0 = [1.0, 0.2]
N = 8
prob = ODEProblem(V, u0, (0., 12.), p0)

function prob_func(prob, i , repeat)
    Rate = ConstantInterpolation(y[i], x[i], extrapolate=true)
    tspan = (tspan_data[i][1], tspan_data[i][2])
    remake(prob; u0 = initial_conditions[i], tspan = (tspan_data[i][1], tspan_data[i][2]), p = prob.p)
end

enprob = EnsembleProblem(prob, prob_func = prob_func)

sim = solve(enprob, Tsit5(), trajectories = N)

I would like to estimate the parameter set (p[1] in this example) with different initial conditions, different sample time-points, different tspan and different Rate functions. @ChrisRackauckas kindly help

2

It is unclear (to me) what you want and your code does not run because many things are missing. Please clarify (and edit your code to reflect):

  1. What is ConstantInterpolation? Where does it come from?
  2. In prob_func what are x and y? They are some global variables not defined here. Same thing with tspan_data.
  3. Please note that your usage of globals is very confusing and likely introduces a bug. You define Rate in line 1, then reuse it inside V. It seems that you try to set it in prob_func but that just creates a new local variable and does not touch the global Rate. I think you should just put that Rate thing into the parameters of the solver.
  4. You talk about estimation of something but your code does not show any comparison to some baseline. Please provide the bigger picture
  5. Why do you need the EnsembleProblem? From your prob_func, it seems you want to do something piecewise?
3

Hi @abraemer, thank you for your response.

  1. ConstantInterpolation is a function from the library DataInterpolations. It is used to make a continuous function between x and y, where x is time and y is volume. So Rate(t) is a flowrate.

  1. x and y are both vectors for different flowrate for different experimental runs.

  2. Yes, that is correct. I can remove the first line in the code for clarity.

  3. I would like to estimate the parameter p[1] for different experimental runs and hence using EnsembleProblem.

4

Just put it in the parameters.

function V(du, u, _p, t)
    p, Rate = _p
    du[1] =  Rate(t)
    du[2] = p[1]*u[1] + (Rate(t) / u[1]) *u[2]
end

p0 = [1.]
u0 = [1.0, 0.2]
N = 8
prob = ODEProblem(V, u0, (0., 12.), p0)

function prob_func(prob, i , repeat)
    Rate = ConstantInterpolation(y[i], x[i], extrapolate=true)
    tspan = (tspan_data[i][1], tspan_data[i][2])
    remake(prob; u0 = initial_conditions[i], tspan = (tspan_data[i][1], tspan_data[i][2]), p = (prob.p, Rate))
end

enprob = EnsembleProblem(prob, prob_func = prob_func)

sim = solve(enprob, Tsit5(), trajectories = N)
5

@ChrisRackauckas Thank you for the help