Kinetic parameter estimation for fixed bed reactor

I am trying to find the kinetic parameter for the fixed bed reactor for the trireforming of methane having the catalyst weight of w=.02 gram and the initial molar flow rate n0 and final molar flow rate n_Measur is known (mol/sec). There is three reaction for steam reforming, dry reforming oand water gas shift reaction which i am trying to use with total species are CH4,CO2,H2O,O2,N2 (inert),H2 and CO. The code is divided into 4 parts error function, ode function, main function and parameter function.

%% ---------------- Parameter Estimation Main -----------------------------
format long

par = par_kinetic();

%% load data from experimental (to find k and n)
T= 923.15 % Temperature in kelvin
w_Measu     = [0 .02];  % weight of catalyst in gram .02

%% nC_Measu    final molar flow rate (mol/min) CH4	CO2	H2O	O2	N2	H2	CO

nC_Measu = [3.59134E-06	6.12004E-06	1.30813E-08	0	8.32244E-05	2.06871E-05	9.36085E-06
3.55089E-06	2.33395E-05	1.14685E-06	0	6.36501E-05	1.57435E-05	1.57362E-05
2.38098E-05	1.28571E-05	1.95838E-07	0	4.59249E-05	2.74115E-05	1.94572E-05
2.44912E-05	6.85368E-06	9.35903E-06	0	5.61497E-05	1.54736E-05	2.05413E-05
2.43895E-05	1.35797E-05	2.44423E-06	0	4.28068E-05	3.64595E-05	1.91442E-05
4.39838E-06	1.60087E-05	1.7746E-05	0	7.93945E-05	2.17505E-05	1.48277E-05
3.99639E-05	9.70848E-06	5.51779E-06	0	3.64002E-05	3.71545E-05	2.40885E-05
3.14698E-06	1.47095E-05	9.1016E-06	0	6.82368E-05	1.59402E-05	1.20256E-05
3.82879E-05	1.13429E-05	9.03076E-06	0	3.79525E-05	4.11719E-05	2.30244E-05
2.33371E-05	4.29807E-06	1.20675E-05	0	4.97604E-05	3.41758E-05	1.49129E-05
2.12084E-05	2.26973E-05	5.16778E-06	0	3.25403E-05	3.64627E-05	2.14637E-05
4.15565E-06	1.17225E-05	1.40098E-05	0	8.54965E-05	1.46762E-05	1.73545E-05
3.62277E-07	1.64866E-06	0.000153265	0	5.95614E-06	2.11172E-06	1.0817E-06
4.14227E-05	1.30693E-05	7.91973E-08	0	2.5918E-05	4.07423E-05	2.11796E-05
2.38667E-05	2.41786E-06	1.15062E-05	0	5.91536E-05	2.05798E-05	1.30274E-05
2.24993E-05	1.59107E-05	1.68287E-05	0	4.01769E-05	2.02628E-05	2.09384E-05
2.05472E-05	5.95366E-06	3.22445E-06	0	6.23875E-05	3.47837E-05	1.47167E-05
1.94723E-05	2.08046E-05	3.03531E-07	0	4.21599E-05	3.1577E-05	2.33277E-05
1.87001E-05	1.04099E-05	1.87612E-05	0	4.61311E-05	2.91735E-05	1.76068E-05
];
%% Initial Condition for variable MOLAR FLOW RATE(mol/min) CH4	CO2	H2O	O2	N2	H2	CO 
n0 = [1.11607E-05	5.58036E-06	1.22768E-05	6.69643E-07	8.19196E-05	0	0
1.11607E-05	2.79018E-05	1.22768E-05	6.69643E-07	5.95982E-05	0	0
2.79018E-05	1.67411E-05	2.23214E-05	2.23214E-07	4.44196E-05	0	0
2.79018E-05	1.67411E-05	2.23214E-06	1.11607E-06	6.36161E-05	0	0
2.79018E-05	1.67411E-05	2.23214E-05	1.11607E-06	4.35268E-05	0	0
1.11607E-05	1.67411E-05	1.22768E-05	2.23214E-07	7.12054E-05	0	0
4.46429E-05	1.67411E-05	1.22768E-05	2.23214E-07	3.77232E-05	0	0
1.11607E-05	1.67411E-05	1.22768E-05	1.11607E-06	7.03125E-05	0	0
4.46429E-05	1.67411E-05	1.22768E-05	1.11607E-06	3.68304E-05	0	0
2.79018E-05	5.58036E-06	2.23214E-05	6.69643E-07	5.51339E-05	0	0
2.79018E-05	2.79018E-05	2.23214E-05	6.69643E-07	3.28125E-05	0	0
1.11607E-05	1.67411E-05	2.23214E-06	6.69643E-07	8.08036E-05	0	0
1.11607E-05	1.67411E-05	2.23214E-05	6.69643E-07	6.07143E-05	0	0
4.46429E-05	1.67411E-05	2.23214E-05	6.69643E-07	2.72321E-05	0	0
2.79018E-05	5.58036E-06	1.22768E-05	2.23214E-07	0.000065625	0	0
2.79018E-05	2.79018E-05	1.22768E-05	2.23214E-07	4.33036E-05	0	0
2.79018E-05	5.58036E-06	1.22768E-05	1.11607E-06	6.47321E-05	0	0
2.79018E-05	2.79018E-05	1.22768E-05	1.11607E-06	4.24107E-05	0	0
2.79018E-05	1.67411E-05	1.22768E-05	6.69643E-07	5.40179E-05	0	0
];

%% Initial Guess for parameter (G1)
p0 = [32
12
118]; 
p1=p0;

%% Solver - Parameter Optimization
for i = 1 : 10
   
% minimization step
options = optimset('MaxFunEvals', 5000,'Display', 'iter', 'TolX', 1e-10, 'TolFun', 1e-10, 'MaxIter', 5000);
[pmin, resnorm, exitflag, output] = fminsearch(@(p) er_fun(w_Measu,nC_Measu,n0,p),p0,options);
% options = optimoptions('lsqnonlin','Display','iter','Algorithm','levenberg-marquardt');
% [pmin,resnorm,residual,exitflag,output] = lsqnonlin(@(p) er_fun(w_Measu,nC_Measu,n0,p),p0,[],[],options);
% loop with initial condition
p1 = pmin;
%n_Calc;
i

end

%% display value
disp('Estimated parameters p(i):');
disp(pmin)
disp('residual norm:');
disp(resnorm)

%% Estimated parameter
par.p = pmin;

% %% ode solver
 %% Calculate molar flow rate using estimated parameters
wspan = [0 0.02];
% Initialize n_Calc to match the dimensions of n0
n_Calc = zeros(size(n0));

% Loop over each row of n0
for i = 1:size(n0, 1)
    [~, n_Calc_temp] = ode45(@(w, n) model_fun(w, n, par.p, par), wspan, n0(i, :));   % dNa/dw =ra
    n_Calc(i, :) = n_Calc_temp(end, :);  % Store the final concentration
end

%% Display molar flow rate
disp('Molar flow rate (mol/min):');
disp(n_Calc);

The main function is use the parameter function.

%% ----------------------- Model Function ---------------------------------
function dy_Calc  = model_fun(t_Calc,y_Calc,p,par)

dy_Calc = zeros(size(y_Calc));
PT = 1;
T=923.15;
R = 8.314;

Ksr=1.198*10^17*exp(-26830/T);
Kdr=6.780*10^18*exp(-31230/T);
Kwgs=10^((2078/T)-2.029);
%% parameter
k1 = p(1);
k2 = p(2);
k3 = p(3);

%FN2=0;  % Check this value the inert species
% FT=f(1)+f(2)+f(3)+f(4)+f(5);
FT=sum(y_Calc);

pCH4=y_Calc(1)*PT/FT;
pCO2=y_Calc(2)*PT/FT;
pH2O=y_Calc(3)*PT/FT;
pO2=y_Calc(4)*PT/FT;
pN2=y_Calc(5)*PT/FT;
pH2=y_Calc(6)*PT/FT;
pCO=y_Calc(7)*PT/FT;

%% rate expression
r1=k1*pCH4*(1-((pCO*pH2^3)/(pCH4*pH2O*Ksr)));
r2=k2*pCH4*(1-((pCO^2*pH2^2)/(pCH4*pCO2*Kdr)));
r3=k3*((pCO*pH2O/pH2)-(pCO2/Kwgs));

% dy/dt=ra
dy_Calc(1) = - r1-r2; % CH4
dy_Calc(2) = r3-r2; % CO2
dy_Calc(3) = -r3-r1; % H2O
% dy_Calc(4) = 0; % O2
% dy_Calc(5) = 0; % N2
dy_Calc(6) = 3*r1+2*r2+r3; % H2
dy_Calc(7) = -r3+r1+2*r2; % CO

return

Model function defines the rate of reaction and the ode equation

% ---------------------------- error function ----------------------------
function  [er,n_Calc]= er_fun(w_Measu,nC_Measu,n0,p)

par = par_kinetic();
%tspan = t_Measu;
%% Solve model with estimated/chosen parameter
wspan = [0 0.02];  %weight of catalyst=tcalc
% [w_Calc, n_Calc] = ode23s(@(t_Calc,nC_Calc) model_fun(t_Calc,nC_Calc,p,par),wspan,n0);

[w_Calc, n_Calc] = ode45(@(w, n) model_fun(w, n, p, par), wspan, n0);   % dNa/dw =ra


    %% ncalc i have to giive to equation for ra

 er = norm(((n_Calc(:,1) - nC_Measu(1,:)).^2+(n_Calc(:,2) - nC_Measu(2,:)).^2 +(n_Calc(:,3) - nC_Measu(3,:)).^2+(n_Calc(:,6) - nC_Measu(6,:)).^2));

end

error function minimizes the objective.

%% ---------------------parameter value------------------------------------
function par = par_kinetic()

%% Model details
par.nofeqns = 3; % number of model equation
par.nVar = 7;
par.unkPar = 3;

return

However, whenever i am runnning the code for all the calculated molar flow rate its coming the be NaN. I am not able to understand what I am doing wrong. Any help will be greatly appreciated. Thanks !!

[Torsten]

n0 for H2 is 0 and you divide by pH2 ...

[MOHD BELAL HAIDER]

Hydrogen is not given in the reactant. That's why the initial hydrogen molar flow rate is zero.

[Torsten]

r3=k3*((pCO*pH2O/pH2)-(pCO2/Kwgs));

pH2 is 0 and thus you divide by 0.

[MOHD BELAL HAIDER]

Thanks for pointing out.. However, what i am trying to do, In the context of the model_fun function, y_Calc represents the current state of the system at a specific time or weight, as calculated by the ode45 function. It’s a vector that contains the molar flow rates of each species at that particular time.

When FT=sum(y_Calc); is called inside model_fun, it’s summing over the molar flow rates of all species at that specific time or weight to get the total molar flow rate FT.

So, y_Calc does not represent the initial conditions n0. Instead, it represents the state of the system at a specific time while the ode45 function is solving the system of ODEs.

Kindly, Let me know if that is okay and if the program i am trying to do is making sense. Thanks a lot !

[Torsten]

ode45 computes the time derivatives first with the initial values for the unknowns and will encounter a division by zero when it computes r3. Shall I run the code and show you that r3 will become NaN ? No, I think you can test it yourself by removing the semicolon behind the line

r3=k3*((pCO*pH2O/pH2)-(pCO2/Kwgs));

And remember that ode45 must exactly supply the output for the flow rates at the times when the measurements were done. Thus setting

wspan = [0 0.02];

for the ode integrator is wrong. wspan must be the 20x1 vector when the measurement data were collected.

[MOHD BELAL HAIDER]

Thanks for the quick response. I understood now what I am doing wrong.. actually, I am quite new in matlab and modified the matlab code which was previously answered by you and @star strider on the community group. The most of the parameters estimation is linked with varying time and concentration and based on that kinetic parameters estimation was done. For fixed bed reactor where catalyst weight is fixed and initial components are varying, I am unaware how to modify this since very little information is given on the mathwrks for parameters estimation for FBR. Will you please, kindly suggest me the modification which I need to do. I Will really appreciate anykind of help.