MATLAB: Fitting parameters in ODE for a kinetic reaction

MATLABode45 fit parameters star stride parameter estimation

Hello. Ok, so I'm new to matlab and I've a question regarding parameter estimation for a kinetic model, for a plug flow reactor with 0.7m.
I have 2 different reactants and their concentrations are c1 (Oxygen) and c2(Hydrogen) with the porpuse of forming c3(Water). The initial concentration of c1 and c2 are known, being c2 constant along the reactor. The only data provided is the initial and final concentration of oxygen for the different experiments.
In this study 5 experiments were perfomed at different tempartures were (475;521;562;579;593 kelvin), and the only data available is the final concentration of oxygen in the system in function of the temperatures, that the experiment was conducted.
This means that , the initial concentration of c1 is 3971
and at the end of the experimente 1 (475k), the concentration of c1 is 3605.8
at the end of the experiment 2 (521k), the concetration of c1 is 2304.9 and then on respectively.
With the respective mass balance being determined in the following way:
    dcdx= -((S * rho_cat)/Q)* theta(1) * exp((-theta(2))/(R*Te))*c(1)*c(3) ;
we want to estimate the different "theta" .
We tried to use the answer provided by Star Stride in Monod kinetics and curve fitting, however due to our data being only the final concentration of oxygen in function of the temperature and not the lenght of the reactor, we couldnt solve it.
Hope someone can help us, i leave the code attached the question

Best Answer

My parameter estimation code will only work for dynamic integrations with respect to time.
If you want to estimate parameters based on the final time, integrate the differential equations and then use one of the parameter estimation functions (probably lsqcurvefit) to estimate the parameters.
The integration: 
syms c1(t) c2(t) S rho_cat Q a1 a2 a3 a4 c10 c20 R T t 
Eq1 = diff(c1) == -((S * rho_cat) / Q) * a1* exp(-a2 / (R * T)) * c1 * 77;
Eq2 = diff(c2) == -((S * rho_cat) / Q) * a3 * exp(-a4 / (R * T)) * c2 * (77.^(1.5)) ;
[C1,C2] = dsolve(Eq1,Eq2, c1(0) == c10, c2(0) == c20)
C1C2 = matlabFunction([C1;C2], 'Vars',{[a1,a2,a3,a4],t,c10,c20,S,rho_cat,Q, R, T})
producing: 
C1 =
c10*exp(-(77*S*a1*rho_cat*t*exp(-a2/(R*T)))/Q)
C2 =
c20*exp(-(5943276032390893*S*a3*rho_cat*t*exp(-a4/(R*T)))/(8796093022208*Q))
C1C2 = @(in1,t,c10,c20,S,rho_cat,Q,R,T)[c10.*exp((S.*in1(:,1).*rho_cat.*t.*exp(-in1(:,2)./(R.*T)).*-7.7e+1)./Q);c20.*exp((S.*in1(:,3).*rho_cat.*t.*exp(-in1(:,4)./(R.*T)).*(-6.756722578291934e+2))./Q)]
Note that ‘in1’ is the parameter vector of the ‘a’ values (in order). You need to create another anonymous function to put ‘C1C2’ into a form that lsqcurvefit can use, likely something like: 
@(in1,T) C1C2(in1,t,c10,c20,S,rho_cat,Q,R,T)
Remember to include the initial conditions ‘c10’ and ‘c20’ and ‘t’ (probalbly the end time), as well as the rest of the parameters.
I leave the rest to you.
Related Question