Kinetic Parameter Estimation using ode15s and lsqcurvefit (local minimum warning)

Hi, I am trying to estimate 4 parameters theta(1)....theta(4).

my odes are here as follows:

function dC=DifEq(t,c,theta)
    time_sec = 3600.0;
    
    % Partial pressure calculation
    R = 0.08205; % L atm/mol/K
    Tgas = 500 + 273.15; %Reactor Temperature (Kelvin)
    Rgas = 8.3145; %Gas Constant (J/mol/K)
    
    pM = c(1) * R * Tgas;
    pE = c(2) * R * Tgas;
    pB = c(3) * R * Tgas;
    pH = c(4) * R * Tgas;
    
    deltaG1 = 0.11073; %GJ/kmol
    deltaG2 = -0.0878439; %GJ/kmol
    deltaG3 = 0.00810928; %GJ/kmol
    
    %Reaction Equilibrium Constant
    Kp1 = exp(-deltaG1*10^6/Rgas/Tgas);
    Kp2 = exp(-deltaG2*10^6/Rgas/Tgas);
    Kp3 = exp(-deltaG3*10^6/Rgas/Tgas);
    
    %reaction rate equations
    r1 = exp(-theta(1)*c(5)) * theta(2) * (pM^2 - (pE*pH^2)/Kp1);
    r2 = exp(-theta(1)*c(5)) * theta(3) * (pE^3 - (pB*pH^3)/Kp2);
    r3 = exp(-theta(1)*c(5)) * theta(4) * (pM - pH^2/Kp3);
    
    dcdt=zeros(6,1);
    dcdt(1)= (1/time_sec)*(-2*r1-r3);
    dcdt(2)= (1/time_sec)*(r1-3*r2);
    dcdt(3)= (1/time_sec)*r2;
    dcdt(4)= (1/time_sec)*(2*r1+3*r2+2*r3);
    dcdt(5)= (1/time_sec)*r3;
    dcdt(6)= 0;
    dC=dcdt;
end

This function is called in the following file:

function C=kinetics(theta,t)
rhogas = 0.0157564; %kmol/m3
% Initial Mole fraction
yM_0 = 0.5;
yE_0 = 0;
yB_0 = 0;
yH_0 = 0;
yN_0 = 1 - (yM_0 + yE_0 + yB_0 + yH_0);
% Initial concentration
cM_0 = rhogas* yM_0; %kmol/m3
cE_0 = rhogas* yE_0; %kmol/m3
cB_0 = rhogas* yB_0; %kmol/m3
cH_0 = rhogas* yH_0; %kmol/m3
cC_0 = 0;
cN_0 = rhogas* yN_0; %kmol/m3
c0 = [cM_0; cE_0; cB_0; cH_0; cC_0; cN_0];
[~,Cv]=ode15s(@(t,y) DifEq(t,y,theta),t,c0);
C_ch4= Cv(:,1);
C_c6h6= Cv(:,3);
C_h2= Cv(:,4);
%velocity
ugas = 0.0437733; %m/s
% Reactor area
Areac = 7.855e-5; %m2
F_ch4= (C_ch4 * ugas * Areac * 3600)/(9.7973e-5);
F_c6h6= (C_c6h6 * ugas * Areac * 3600)/(6.0354e-7);
F_h2= (C_h2 * ugas * Areac * 3600)/(7.7289e-5);
F=[F_ch4 F_c6h6 F_h2];
%C=Cv(:,[1,3:4]);
C=F;
end

The above two files are called here in the Lsqcurvefit function:

clc
clear
t=[300 
   600 
   900 
   1200 
   1500 
   1800 
   2100 
   2400 
   2700 
   3000 
   3300 
   3600];
c=[ 0.884432367	0.324437395	0.264387415
    0.963445551	0.192282942	0.12149159
    0.967386976	0.145677499	0.087767099
    0.966478665	0.104783894	0.069448077
    0.967632495	0.090827748	0.069047991
    0.980013654	0.083761307	0.055270652
    0.989909288	0.077522552	0.051081572
    0.985046526	0.068823873	0.044012757
    0.982500604	0.068821302	0.042979821
    0.987528139	0.078424815	0.050897805
    0.984648739	0.078326295	0.047444151
    0.991200011	0.063706374	0.039924062]; 
theta0=[1.54; 1.8413; 6.615e7; 0.0829];
lb = [0,0,10^5,0];
ub = [10,1000,10^10,1000];
options = optimoptions(@lsqcurvefit,'Algorithm','trust-region-reflective','OptimalityTolerance',eps); 
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,lb,ub,options);
fprintf(1,'\tCatalyst Activity:\n')
k1 = 1
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
fprintf(1,'\tRate Constants:\n')
for k1 = 2:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(0, max(t));
Cfit = kinetics(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Molar Flowrate')
legend(hlp, 'Methane', 'Benzene', 'Hydrogen', 'Location','N')

I ran my code and i keep geeting the following solution. Based on suggestions for this warning, I scaled my variables and also reduced the optimality tolerance. However, I am still getting the following warning. Moroever, the trend of my components is opposite of exp data. I was wondering if someone could please help me.

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.
<stopping criteria details>
	Catalyst Activity:
k1 =
     1
		Theta(1) =  1.53799
	Rate Constants:
		Theta(2) =  1.83961
		Theta(3) = 66231268.64922
		Theta(4) =  0.08272
>> 

[Alan Weiss]

You say that you rescaled the problem, but your results show that Theta(3) is many orders of magnitude larger than the other Theta components. I would scale Theta(3) by, say, 1e8 (or by theta0(3)) internally to the objective.

It is not at all clear to me that the result is incorrect in any way. What are you looking for that would make you more satisfied?

Alan Weiss

MATLAB mathematical toolbox documentation

[CM007]

Sorry, I did not post the model vs experimental plot. I am trying to estimate the paramaters (thetas) by fitting my model to the experimental data. However, the model is not able to fit the experimental data. After encountering the above mentioned warning, I used Multistart based on suggestions here. However, the model still doesnt fit the data and it gives garbage results.