Parametric estimation using lsqcurvefit with differential equation solver using ode23s for reactions

Hello Everyone,

I am trying to find reaction parameters (rate constants, equlibrium constants, adsoprtion coeeficients) for methanol synthesis from hydrogenation of CO2.

CO + 2H2 <=> CH3OH (reaction1)

CO2 + H2 <=> CO +H2O (reacion 2)

CO2 + 3H2 <=> CH3OH +H2O (reaction 3)

I have 3 interdependent rate expressions for this set of expressions.

Below is the MATLAB code for calling lsqcurvefit function;

The problem is my concentrations are not changing and no reaction is happening. Can anyone please help?

CODE

clear all;
clc;
R = 8.314; %gas constant
%Original Data for curve fitting, Temperature and yield (Total 9 points)
T = [131.6762, 167.1254, 221.4437, 246.8176, 272.0168 , 297.0263, 322.0475, 372.2783, 421.3443];
Y = [0.073197828, 0.168737656, 0.999711982, 1.699103028, 2.165878045,2.135018927, 1.825008229, 1.114384463, 0.512261356];
plot (T, Y, 'bo')
%------------------------
%guess values for parameters in the function(total 18)
guess = [4.0638e-6 11695 8.3965e-11 118270 3.5408e12 19832 9.0421e8 112860 1.7214e-10 81287 6.1221e-13 125226 1.5188e-33 266010 2.5813e10 26788 4.3676e-12 1.1508e5];

[coeff, resnorm, residual] = lsqcurvefit(@fitCO2, guess, T, Y); 
%fit CO2 function given below after code

%all the coefficients
a11 = coeff(1); b11 = coeff(2);
a12 = coeff(3); b12 = coeff(4);
a13 = coeff(5); b13 = coeff(6);
a21 = coeff(7); b21 = coeff(8);
a22 = coeff(9); b22 = coeff(10);
a23 = coeff(11); b23 = coeff(12);
a31 = coeff(13); b31 = coeff(14);
a32 = coeff(15); b32 = coeff(16);
a33 = coeff(17); b33 = coeff(18);
%------------------------
%Three reactions involved
%and these are all the parameters it should give as output 
%Rxn #1
kaf = a11.*exp(-b11./(R.*T));
kCO = a12.*exp(-b12./(R.*T));
k1 = a13.*exp(-b13./(R.*T));
%----------------------------
%Rxn #2
kbf = a21.*exp(-b21./(R.*T));
%----------------------------
kCO2 = a22.*exp(b22./(R.*T));
k2 = a23.*exp(-b23./(R.*T));
%Rxn #3
kcf = a31.*exp(-b31./(R.*T));
k3 = a32.*exp(b32./(R.*T));
kH = a33.*exp(b33./(R.*T));

%------------------------
Y2 = residual+Y;
hold on;
plot(T,Y2, 'm-')

FITCO2 functiion

function Y2 = fitCO2(x0, T) % function fitCO2
% x0 is the paramters vectors which has 18 terms
% Now my functions are differential equations themselves so I need to use
% another function ode23s to solve differential equations
global yfun2; %variable to store desired values
yfun2 = zeros(1,9);
% here is the setup for differential equation solver
tspan = 0:60; % in seconds
Co = [0 0 0 0.3 0.7];
[t, C] = ode23s(@conc, tspan, Co);

function yfun = conc(t, C)
    
    
  % parameteric coefficients are labeled here in terms of x0
  % all parameters are in arrhenius form a*exp(b/RT)
R = 8.314; % J/mol/K , Gas Constant
a11 = x0(1); b11 = x0(2);
a12 = x0(3); b12 = x0(4);
a13 = x0(5); b13 = x0(6);
a21 = x0(7); b21 = x0(8);
a22 = x0(9); b22 = x0(10);
a23 = x0(11); b23 = x0(12);
a31 = x0(13); b31 = x0(14);
a32 = x0(15); b32 = x0(16);
a33 = x0(17); b33 = x0(18);
%----------------------------
%-----------------------------------
%Variables for concentrations the program needs to solve for
   p_CH3OH = C(1);
   p_CO = C(2);
   p_H2O = C(3);
   p_CO2 = C(4);
   p_H2 = C(5);
%-----------------------------------
% all the parameters are assigned here
%Rxn #1
kaf = a11.*exp(-b11./(R.*T));
kCO = a12.*exp(-b12./(R.*T));
k1 = a13.*exp(-b13./(R.*T));
%----------------------------
%Rxn #2
kbf = a21.*exp(-b21./(R.*T));
kCO2 = a22.*exp(b22./(R.*T));
k2 = a23.*exp(-b23./(R.*T));
%----------------------------
%Rxn #3
kcf = a31.*exp(-b31./(R.*T));
k3 = a32.*exp(b32./(R.*T));
kH = a33.*exp(b33./(R.*T));

%----------------------------
%These are the three functions r1, r2 and r3 in which we need to find
%the parameters
ads = (p_H2.^0.5) + (kH).*p_H2O + kCO.*p_CO.*(p_H2).^0.5+ (kCO.*kH).*p_CO.*p_H2O + kCO2.*p_CO2.*(p_H2).^0.5 + (kCO2.*kH).*p_CO2.*p_H2O;
r1 = (kaf.*(kCO.*p_CO.*((p_H2)^1.5) - (k1.*p_CH3OH./((p_H2).^0.5))))./ads;
r2 = (kbf.*(kCO2.*p_CO2.*p_H2 - k2.*p_H2O.*p_CO))./ads;
r3 = (kcf.*(kCO2.*p_CO2.*(p_H2.^1.5) - (k3.*(p_H2O.*p_CH3OH)./(p_H2.^1.5))))./ads;

%Differential equation solver
yfun = [r1+r3;          %CH3OH
        -r1+r2;             %CO
        r2+r3;          %H2O
        -r2-r3;      %CO2
        -2*r1-r2-3*r3]; %H2
% So now what I want to do is to solve these set of time dependent 
%differential equations for each temperature and get the last point of
%concentration of CH3OH at each temperature, this way I should get 9
%CH3OH concentration points as Y2 to fit them to initial data Y againts T.
    
%yfun size is not 1x9 so storing desired values in yfun2
yfun2 = [yfun(1,1), yfun(1,2),yfun(1,3),yfun(1,4),yfun(1,5),yfun(1,6),yfun(1,7),yfun(1,8),yfun(1,9)];
%next three steps are done to get 1x5 vector for yfun as required by MATLAB
% I didnot use yfun in the final solution
yfun = yfun(:);
yfun = [yfun(1,:), yfun(6,:),yfun(11,:),yfun(16,:),yfun(21,:)];
yfun = yfun';

end

Y2 = yfun2;

end