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

**URL:** <https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [August 29, 2020, 8:27am UTC](https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727 "2020-08-29T08:27:38Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [August 29, 2020, 8:27am UTC](https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727/1 "2020-08-29T08:27:38Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/9/6/963b293dc6394778cb0360d1bcd1b3f308ed4038.png)  
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](https://global.discourse-cdn.com/julialang/original/3X/5/e/5ee5ad2337417b9b6fb4f81cac9abce3b1d9e193.png)  
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](https://global.discourse-cdn.com/julialang/original/3X/a/8/a8c54a5015c83c07be1f6a320df2001c2ff34812.png)

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](https://discourse.julialang.org/t/differential-equation-parameter-estimation-using-only-a-subset-of-variables/23423) by using save\_idxs to limit the saved variable, the problem seems to be solved. Here is the example:

```julia
# 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:

```julia
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())

```

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/4/f4d93dc2917e6487b96085bd13cb05398d0ddc79.png)

---

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [August 29, 2020, 10:13am UTC](https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727/2 "2020-08-29T10:13:54Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 29, 2020, 3:44pm UTC](https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727/3 "2020-08-29T15:44:58Z")

</div>

Take a look at [https://diffeqflux.sciml.ai/dev/examples/optimization\_ode/](https://diffeqflux.sciml.ai/dev/examples/optimization_ode/)

---

<div class="post-metadata">

**Author:** ![Fred\_He](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fred_he/32/14298_2.png) [@Fred\_He](https://discourse.julialang.org/u/Fred_He)\
**Post date:** [August 29, 2020, 8:34pm UTC](https://discourse.julialang.org/t/param-estimation-with-different-initial-condition-how-to-do-it/45727/4 "2020-08-29T20:34:53Z")

</div>

Thank you Chris.