ode15s giving integration tolerance error and loop stopping at zeroth iteration

% 1) Enter Model inputs -  Constants and parameters% 
kof= 0.30;   
KsPO=0.15;
Ksmf=0.25;
khy= 0.55;  
Ks1=5;
Ks2=0.67;
 
A=0;
B=0.00025; 
K=0.25;
v=720;
% Pav calculated at air inlet conditions of 75 % RH and at 301.15K
% syms Pav in mm Hg
 
% Pav=28.84; %Vary from 21.63 to 24.3
Pt=760;
% rdxb
%syms kd (d-1)
kd= 0.02;
% Tf (K)
%syms cpf cpb cpL(KJ/kg K) Tref (K)
cpf= 3.5;
cpb= 1.3;
cpL= 4.18;
Ks=0.005;
kos= 0.2;
Ksb=100*10^-6;
% Liquid and gas phases 
Q = 115.2;  % vary from 57.6 to 172.8 (89 LPM = 128m3/d)
Tref= 298.15; 
Vol = 0.08;   %   0.3801;  
cps=3.8;
Vp=  0.05; %0.2851; %;

%  2.Specify  state variables and independent variables by using syms
syms xf(t) xb(t) mf(t) Tf(t) xs(t) mL(t) s(t) xbs(t) O2(t) CO2(t) N2(t) aV(t) rof 
syms rhy rdxb ref ref_h rhy_h ros rdxbs reL F_Tf xt Y_H 

% 3.Enter rate equations
%arhn
F_Tf = exp(0.024*(Tf(t)-273.15-35));

%eqn1
 rhy= khy*xf(t)*(1/(Ks1+(xf(t)/xb(t))))*((mf(t)/xb(t))/(Ks2+(mf(t)/xb(t))));
 Pair= O2(t)+CO2(t)+aV(t)+N2(t) ;
 xt = xf(t)+xb(t) ;

%eqn2
% rof = kof*xb*F_Tf*min([(xf/xt)/(Ksf+(xf/xt)),(O2/Pair)/(KsPO+(O2/Pair)),(mf/xf)/(Ksmf+(mf/xf))]);
% Instead of 5 try with 4 0r 3
rof =kof*xb(t)*(((O2(t)/Pair)*100)/(KsPO+((O2(t)/Pair)*100)))*((mf(t)/xf(t))/(Ksmf+(mf(t)/xf(t))))*((xf(t)/xt)/(Ks+(xf(t)/xt))) ;

%eqn3
% ros = kos*xbs* min([(s/mL)/(Ksb+s/mL), O2/(KsPO+O2)]);% check with s/ksb+s
ros = kos*xbs(t)*(((O2(t)/Pair)*100)/(KsPO+((O2(t)/Pair)*100)))*((s(t)/mL(t))/(Ks+(s(t)/mL(t))));

%eqn4
rdxb = kd *xb(t) ;
%eqn5  
rdxbs = kd*xbs(t);
Pfs = 10^((-2238/Tf(t))+ 8.896);
Pav=((aV(t)*Pt)/(O2(t)+CO2(t)+N2(t)+aV(t))) ;
PLs = 10^((-2238/Tf(t))+ 8.896);
%eqn 6 
ref = (A+ B*v)*(((xf(t)+ mf(t))/125)^(2/3))*(1- (Pav/Pfs))*((mf(t)/xf(t))/(K+(mf(t)/xf(t))))*mf(t) ;% xf and mf initial value on denominator ^(2/3)
%eqn 7
reL= (A+ B*v)*(1- (Pav/PLs))*mL(t)*(( mL(t)/ 0.001)^0.67) ; 
%eqn 8
ref_h = 2430*ref;
%eqn 9
rhy_h = rhy * cpf*(Tf(t)-Tref) ;

% 4. Specify the first order differential equations
eqn1 =diff(xf(t),1)== -rof -rhy;
eqn2 = diff(xb(t),1)== 0.5*rof-0.2*rhy -rdxb ;
eqn3 =diff(mf(t),1)== -ref - 4*rhy;
eqn4 = diff(xs(t),1) == 0.5*(rdxb+rdxbs)+0.25*(rof+ros);
eqn5 = diff(s(t),1) == rhy -ros ;
eqn6 = diff(mL(t),1) == 4*rhy - reL ;
eqn7 = diff(xbs(t),1) == 0.5* ros - rdxbs;
eqn8 = Vp*diff(O2(t),1)== ( Q * 8.64- Q*O2(t))- (Vol*2.36*(rof+ ros + rdxb + rdxbs));
eqn9= Vp*diff(CO2(t),1)== ( Q * 0.021- Q*CO2(t))+(Vol*2.55*(rof+ros + rdxb + rdxbs));
eqn10 = Vp*diff(aV(t),1)== ( Q * 1.2 - Q*aV(t)) +(Vol*1000/18*(ref+reL)) ;
eqn11 = Vp*diff(N2(t),1)== ( Q * 30.64- Q*N2(t));
eqn12 = diff(xf*cpf*Tf(t),1) + diff(xb*cpb*Tf(t),1) + diff(mf*cpL*Tf(t),1) == 20000*rof - ref_h - rhy_h ;
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8 eqn9 eqn10 eqn11 eqn12];

% 5. Place the state variables in a column vector.
vars = [xf(t); xb(t); mf(t);  xs(t) ; mL(t) ;s(t); xbs(t); O2(t); CO2(t); N2(t); aV(t); Tf(t);];
% 6. Store the number of original variables for reference
origVars = length(vars);
% 7.Check Differential Index of System and Reduce the differential index ofthe DAEs if required
if (isLowIndexDAE(eqns,vars) == 0)
    [eqns,vars] = reduceDAEIndex(eqns,vars);
end
%8. Convert DAEs to Function Handles
[M,f] = massMatrixForm(eqns,vars);

% 9.Find symbolic parameters that are not specified in the vector of variables DAEvars
pDAEs = symvar(eqns);
pDAEvars = symvar(vars);
extraParams = setdiff(pDAEs, pDAEvars);
M = odeFunction(M, vars);
f = odeFunction(f, vars, Y_H);
Y_H =20000;% 23200
F = @(t,Y) f(t,Y, Y_H);
y0est =[35.25;3;141;1;0;0;1;8.635;0.021;31.35;1.2;303];

% 10. Create an option set that contains the mass matrix M and initial guesses 
% yp0est, and specifies numerical tolerances for the numerical search.
opt = odeset( 'Mass',M,'RelTol', 1e-7, 'AbsTol' , 1e-7);


%11. Solve the system of equations using ode15s
 
[tSol,ySol] = ode15s(F, 0:1:7, y0est,opt);

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.

Can anyone help with converting this into a simulink model as the matlab code isnt working at all?

[Sam Chak]

@Devan, please show how you validated that the mathematical model is correctly described in the code, especially for these four lines:

rhy = khy*xf(t)*(1/(Ks1+(xf(t)/xb(t))))*((mf(t)/xb(t))/(Ks2+(mf(t)/xb(t))));

rof =kof*xb(t)*(((O2(t)/Pair)*100)/(KsPO+((O2(t)/Pair)*100)))*((mf(t)/xf(t))/(Ksmf+(mf(t)/xf(t))))*((xf(t)/xt)/(Ks+(xf(t)/xt)));

ros = kos*xbs(t)*(((O2(t)/Pair)*100)/(KsPO+((O2(t)/Pair)*100)))*((s(t)/mL(t))/(Ks+(s(t)/mL(t))));

ref = (A+ B*v)*(((xf(t)+ mf(t))/125)^(2/3))*(1- (Pav/Pfs))*((mf(t)/xf(t))/(K+(mf(t)/xf(t))))*mf(t);

Answer 1

[Torsten]

The evaluation of the variable "ros" for your vector of initial conditions y0est gives NaN. The reason is that you start with mL = s = 0.

y0est =[35.25;3;141;1;0;0;1;8.635;0.021;31.35;1.2;303];
% 10. Create an option set that contains the mass matrix M and initial guesses 
% yp0est, and specifies numerical tolerances for the numerical search.
opt = odeset( 'RelTol', 1e-7, 'AbsTol' , 1e-7);
%11. Solve the system of equations using ode15s
F(0,y0est) 

ros = NaN

ans = 12x1

-1.9777 0.1300 -14.0506 NaN NaN 4.5648 NaN NaN NaN 843.1863

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[tSol,ySol] = ode15s(@F, 0:1:7, y0est);%,opt)

function dy = F(t,y)

xf = y(1);
xb = y(2);
mf = y(3);
xs = y(4);
mL = y(5);
s = y(6);
xbs = y(7);
O2 = y(8);
CO2 = y(9);
N2 = y(10);
aV = y(11);
Tf = y(12);


% 1) Enter Model inputs -  Constants and parameters% 
kof= 0.30;   
KsPO=0.15;
Ksmf=0.25;
khy= 0.55;  
Ks1=5;
Ks2=0.67;
 
A=0;
B=0.00025; 
K=0.25;
v=720;
% Pav calculated at air inlet conditions of 75 % RH and at 301.15K
% syms Pav in mm Hg
 
% Pav=28.84; %Vary from 21.63 to 24.3
Pt=760;
% rdxb
%syms kd (d-1)
kd= 0.02;
% Tf (K)
%syms cpf cpb cpL(KJ/kg K) Tref (K)
cpf= 3.5;
cpb= 1.3;
cpL= 4.18;
Ks=0.005;
kos= 0.2;
Ksb=100*10^-6;
% Liquid and gas phases 
Q = 115.2;  % vary from 57.6 to 172.8 (89 LPM = 128m3/d)
Tref= 298.15; 
Vol = 0.08;   %   0.3801;  
cps=3.8;
Vp=  0.05; %0.2851; %;

% 3.Enter rate equations
%arhn
F_Tf = exp(0.024*(Tf-273.15-35));
%eqn1
 rhy= khy*xf*(1/(Ks1+(xf/xb)))*((mf/xb)/(Ks2+(mf/xb)));
 Pair= O2+CO2+aV+N2 ;
 xt = xf+xb ;
%eqn2
% rof = kof*xb*F_Tf*min([(xf/xt)/(Ksf+(xf/xt)),(O2/Pair)/(KsPO+(O2/Pair)),(mf/xf)/(Ksmf+(mf/xf))]);
% Instead of 5 try with 4 0r 3
rof =kof*xb*(((O2/Pair)*100)/(KsPO+((O2/Pair)*100)))*((mf/xf)/(Ksmf+(mf/xf)))*((xf/xt)/(Ks+(xf/xt))) ;
%eqn3
% ros = kos*xbs* min([(s/mL)/(Ksb+s/mL), O2/(KsPO+O2)]);% check with s/ksb+s
ros = kos*xbs*(((O2/Pair)*100)/(KsPO+((O2/Pair)*100)))*((s/mL)/(Ks+(s/mL)))
%eqn4
rdxb = kd *xb ;
%eqn5  
rdxbs = kd*xbs;
Pfs = 10^((-2238/Tf)+ 8.896);
Pav=((aV*Pt)/(O2+CO2+N2+aV)) ;
PLs = 10^((-2238/Tf)+ 8.896);
%eqn 6 
ref = (A+ B*v)*(((xf+ mf)/125)^(2/3))*(1- (Pav/Pfs))*((mf/xf)/(K+(mf/xf)))*mf ;% xf and mf initial value on denominator ^(2/3)
%eqn 7
reL= (A+ B*v)*(1- (Pav/PLs))*mL*(( mL/ 0.001)^0.67) ; 
%eqn 8
ref_h = 2430*ref;
%eqn 9
rhy_h = rhy * cpf*(Tf-Tref) ;
% 4. Specify the first order differential equations
dy = zeros(12,1);
dy(1) = -rof -rhy;
dy(2) = 0.5*rof-0.2*rhy -rdxb ;
dy(3) = -ref - 4*rhy;
dy(4) = 0.5*(rdxb+rdxbs)+0.25*(rof+ros);
dy(5) = rhy -ros ;
dy(6) = 4*rhy - reL ;
dy(7) = 0.5* ros - rdxbs;
dy(8) = (( Q * 8.64- Q*O2)- (Vol*2.36*(rof+ ros + rdxb + rdxbs)))/Vp;
dy(9) = (( Q * 0.021- Q*CO2)+(Vol*2.55*(rof+ros + rdxb + rdxbs)))/Vp;
dy(10) = (( Q * 1.2 - Q*aV) +(Vol*1000/18*(ref+reL)))/Vp ;
dy(11) = ( Q * 30.64- Q*N2)/Vp;
dy(12) = (-dy(1)*cpf*Tf-dy(2)*cpb*Tf-dy(3)*cpL*Tf + 20000*rof - ref_h - rhy_h)/(xf*cpf+xb*cpb+mf*cpL);

end

[Sam Chak]

Now, the initial values are non-zero, but the prognosis does not look promising.

Even if the simulation runs without any errors and produces results, when a system designer cannot verify that the mathematical code she wrote is correct based on the derived or standard mathematical model, the validity of the simulated results becomes questionable.

y0est   = [35.25;3;141;1;1;1;1;8.635;0.021;31.35;1.2;303];
opt     = odeset( 'RelTol', 1e-7, 'AbsTol' , 1e-7);
[tSol, ySol] = ode15s(@F, [0 7], y0est);

Warning: Failure at t=3.734031e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.938894e-18) at time t.

plot(tSol, ySol), grid