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

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: 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. 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. 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 is total pressure
r = p .* parP[:,1] .* parP[:,2]
du[:, 2] = -r*mcatL
du[:, 3] = -r*mcatL
du[:, 4] = 1.5r*mcatL
du[:, 1] = mcatL.*r./FA0 # p 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
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, label="simulation")
``````

the result looks reasonably good. Thanks for the community!

1 Like

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)` ?

Thank you Chris.