Optimization of time-dependent parameters with time-series in chemical model

Specific Domains Optimization (Mathematical)
,
1

Hello All,

I have a coupled set of differential equations with some-time-dependent parameters that I would like to optimize.

In the code at the bottom of this post, the k_i terms are fixed kinetic constants, the c[index] terms are concentrations of chemical species, and the s[index] terms are the source terms which are being solved for. S and c are time-dependent. I also have observations at discrete time-steps (0, 1, 2, 3, …, n) of each species. My questions are the following:

  1. For solving the forward problem, The parameters (sources) have values at fixed, in this case discrete, points in time. However, the ODE solver will take time steps between these discrete values (e.g., t = 1.2 versus 1 or 2). So far, I perform an interpolation for the source terms in the time domain by interpolating the parameters in time. Please refer to the function interpSources in the code. Is there a better way to get parameters that are continuous in time e.g., passing a-time-dependent function into ODEProblem instead of an array of parameters?

  2. Also, do you have any recommendations on how to proceed with optimizing these-time-dependent parameters (source-terms) constrained by observational time-series? I tried touring_inference from DiffEqBayes.jl, but that doesn’t seem to work due to the interpolation problem mentioned in Question 1. The examples online are for fixed parameters (e.g., the Lawrence Equations or predator-prey model).

The simplified code is below. Please let me know if you need more clarification, because I know this is a rather convoluted question. Thanks in advance!

using DifferentialEquations, DiffEqBayes, ParameterizedFunctions
using Distributions, Interpolations

end_time = 25.0;
time = collect(1.0:end_time);
num_steps = length(time);
tspan = (1.0 ,end_time);

# define some constants 
const k1, k2, k3 = 7.89, 464.8, 3600*24;

function interpSources(s, t)
    num_time, num_species = size(s)
    species_grid = [i for i =1:num_species];
    grid = (time, species_grid);
    intp = interpolate(grid, s, Gridded(Linear()));
    s_out = intp(t, species_grid);
    return s_out
end #function interpSources

    
function BoxModel(dc ,c ,p ,t)
    s = interpSources(p ,t) # interpolate the source terms (parameters) in time

    dc[1] = s[1] - k1 * c[1] * c[3];
    dc[2] = s[2] + k1 * c[1] * c[3] - k2 * c[2] * c[3];
    dc[3] = s[3] - k1 * c[1] * c[3] - k2 * c[2] * c[3] - k3 * c[3];
end # function BoxModel

# define initial conditions
c_0 = [1700.0 100.0 3.71e-5];

# defining the source-terms (parameters) where row i is time i and col j is species j
sources = [1.05 * ones(num_steps, 1) 1.06 * ones(num_steps, 1) 7.74 * ones(num_steps, 1)];

# define and solve the forward problem
problem = ODEProblem( BoxModel ,c_0 ,tspan ,sources);
sol = solve(problem, alg_hints = "stiff",saveat=time);