Param. estimation with different initial condition, how to do it?

General Usage
1

Hello, I am new to DiffEqParamEstimation.jl package and trying to use it.

The first thing I want to solve with this package is kinetic constant estimation based on plug flow model. It can be simplified to:
image
However, I have only one response, the conversion of A. I don’t measure the flow rate/concentration of B & C. What I am doing now is convert this equation to another form, the conversion vs spaceTime.
image
But in my experiments, I change the total flow rate to measure the conversion, so basically I changed the space time. The problem is the concentration of A is constant while the concentration of B and C are not. As you can see from the tahble below, the concentration of B (CO) has two values.
image

If I do parameter estimation only with the same initial A, B and C concentration, I can use dX/d(W/F0) model. But how can I account for all the initial conditions (in the example, two)?

UPDATE:
inspired by this post: Differential Equation Parameter Estimation Using Only A Subset of Variables by using save_idxs to limit the saved variable, the problem seems to be solved. Here is the example:

# initial condition and function defination
u01 = 1.1.*[1, 1.5, 2, 2.5, 3]*[500E-6, 2500E-6, 0.997]'
u02 = 1.1.*[1, 1.5, 2, 2.5, 3]*[500E-6, 5000E-6, 0.9945]'
u0 = [zeros(row) [u01; u02]]
u0[:, 2:4]./sum(u0[:, 2:4], dims =2)

function pfr!(du,u,p,t)
    parP = u[:, 2:4]./sum(u[:, 2:4], dims =2)*press # p[2] is total pressure
    r = p[1] .* parP[:,1] .* parP[:,2]
    du[:, 2] = -r*mcatL
    du[:, 3] = -r*mcatL
    du[:, 4] = 1.5r*mcatL
    du[:, 1] = mcatL.*r./FA0 # p[3] is initial flow rate of A
end

k0 = 0.02
press = 1.
L = 0.03 # [m]
mcat = 8 # [g]
mcatL =  mcat/L
FA0 = u0[:,2]
#para = [k0, press, FA0, mcat/L]
para = [k0]
prob = ODEProblem(pfr!, u0, (0, L), para)
sol = solve(prob, Rosenbrock23(), isoutofdomain=(y,p,t)->any(x->x<0,y), save_idxs=(1:10), save_everystep=false)

then the cost function is defined as:

t = [0, 0.03]
data = [0.355, 0.31, 0.284, 0.261, 0.242, 0.622, 0.622, 0.542, 0.447, 0.395]
cost_function = build_loss_objective(prob,Rosenbrock23(), isoutofdomain=(y,p,t)->any(x->x<0,y), save_idxs=(1:10), 
    save_everystep=false,L2Loss(t,data), maxiters=10000,verbose=false) 
result = optimize(cost_function,[1.2] , BFGS())

k0 = result.minimizer[1]
para = [k0]
println("k0 = $para")
prob = ODEProblem(pfr!, u0, (0, L), para)
sol = solve(prob, Rosenbrock23(), isoutofdomain=(y,p,t)->any(x->x<0,y), save_idxs=(1:10), save_everystep=false)
scatter(data, label="exp. Results")
scatter!(sol.u[2], label="simulation")

the result looks reasonably good. Thanks for the community!

image
image772×492 12.5 KB

1 Like
2

But what I don’t understand is the cost function. I defined:

t = [0, 0.03]
data = [0.355, 0.31, 0.284, 0.261, 0.242, 0.622, 0.622, 0.542, 0.447, 0.395]
cost_function = build_loss_objective(prob,Rosenbrock23(), isoutofdomain=(y,p,t)->any(x->x<0,y), 
save_idxs=(1:10),  save_everystep=false,L2Loss(t,data), maxiters=10000,verbose=false)

this two array for cost function, but the solution, sol.u consists of two elements, the value at 0 and the value at 0.03. I only give the value at 0.03. So what is behind the lost function L2Loss(t,data) ?

3

Take a look at https://diffeqflux.sciml.ai/dev/examples/optimization_ode/

4

Thank you Chris.