MATLAB: Predicting the kinetic constant of a reaction based on experimental data

chemistrykineticsMATLABrate law

I have the following optimised reaction based on experiments: A + 1.7B -> C (A = Rifamycin Oxazine; B = Piperazine; C = Rifampicin)
My postulated rate law is therefore:
r_A = -k*C_A*C_B^1.7
r _B= -k*C_A*C_B^1.7
r _C= k*C_A*C_B^1.7
Using my experimental data, I have used the following code to try and accurately find k (the goal is to get an R^2 > 90%):
function API2
function C=kinetics(theta,t)
c0=[0.752;1.278;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(3,1);
    dcdt(1)=-theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dcdt(2)=-theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dcdt(3)=theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dC=dcdt;
    end
C=Cv;
end
T = [0 1 2 5 10 15];
t = T';
%y values for a
a_ydata = [0.752 0.0596 0.0596 0.0596 0.0502 0.0424];
A_Ydata = a_ydata';
%y values for b
b_ydata = [1.278 0.378 0.101 0.101 0.085 0.072];
B_Ydata = b_ydata';
c_ydata = [0 0.692 0.692 0.692 0.702 0.71];
C_Ydata = c_ydata';
c = [A_Ydata B_Ydata C_Ydata];
theta0=[0.5];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
h = plot(t, c,'.');
set(h,{'Marker'},{'s';'d';'^'},{'MarkerFaceColor'},{'r';'b';'k'},{'MarkerEdgeColor'},{'r';'b';'k'});
hold on
hlp = plot(tv, Cfit,'LineWidth',1.5);
set(hlp,{'Color'},{'r';'b';'k'});
hold off
grid
xlabel('Time (min)')
ylabel('Concentration (M)')
legend(hlp, 'Rifamycin Oxazine', 'Piperazine', 'Rifampicin', 'Location','N')
end
However, when I run the code and produce a graph, there seems to be a clear discrepancy between my simulated curve and my experimental data (especially for piperazine). I have been messing around with the code but can't seem to get it to fit well with my data. I'm not sure what to do. Can anyone help me?

Best Answer

If you have the Global Optimization Toolbox, the code I attached will estimate the parameters with reasonable accuracy, although your model extimates only one parameter in three differential equations. (This is not something I have observde previusly, so it may be appropriate for you to reconsider the model.)
That aside, ga managed to converge, although with somewhat different estimates in each run.
Note — I set the initial conditions as parameters to be estimated in my code.
Typical parameters were: 
	Rate Constants:
		Theta(1) =  2.55140
		Theta(2) =  0.92250
		Theta(3) =  0.99550
		Theta(4) =  0.01400
and 
	Rate Constants:
		Theta(1) = 10.25306
		Theta(2) =  0.88801
		Theta(3) =  0.98375
		Theta(4) =  0.00308
with essentially identical fitness values, producing this plot:
with the second set of parameters.
That is likely as good as it gets for you model.
.
Related Question