Parameter estimation using ode45 and Lsqcurvefit function [Error: objective function returning undefined value at initial point]

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 was able to rectify the last error b changing the ode function ode15s. I ran my code and i keep geeting the following solution. Based on suggestions for this warning, I scaled my ariables 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
>> 

[Star Strider]

When I add this line after the ‘pH’ assignment in ‘DifEq’:

pQ = [pM pE pB pH]

it shows that all of them become NaN by the seventh call to the function. The problem could be in the ‘F’ calculations, since that is what lsqcurvefit attempts to fit to the data.

Your code varys significantly from the original code I wrote for that problem. It is not easy for me to follow it.

.

[CM007]

Hi Star Strider,

F is the molar flowrate (kmol/h) which is calculated ffrom the concentrations c(1),c(3),c(4). I am not sure why it is given NaN after a few iteraitons. I dissected your code into different files to better understand i.e. the kinetics function, Diffeq, and lsqcurvefit. The core of your code is exactly the same. I have changed the material balance equations to what I need.

In dC=DifEq(t,c,theta):

I have the rate model in temrs of partial pressure. So, that is what I am calculating initially. It is then used in the dcdt equations.

In function C=kinetics(theta,t):

I am calculating the initial values for the concentration using density and mole fraction. Then I am exctracting the 3 components concentration (c(1),c(3),c(4)) which I want to match with exp data.SInce, my exp data is in kmol/h, I am calculating those using the concnetration.

Both of this functions are called in lsqcurvefit.