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