Discrete S function error
3 views (last 30 days)
Show older comments
Hi
I recently found a continuous system model for a bioreactor at the link below and I tried to edit the code to make it discrete but for some reason my code isnt working properly. The initial values load at t = 0 but at the next time step some output just become NaN and/or inf. The equations are correct i verified that with the original working continuous model. So please if possible could see my code below and point out my mistakes.
Thank you so much
Ubaid Imtiaz
Link to the original working model
My code
function [sys,x0,str,ts,simStateCompliance] = bmodeldisc2(t,output,u,flag)
switch flag, case 0, [sys,x0, str,ts,simStateCompliance,output]= mdlInitializeSizes;
case 2,
sys=mdlUpdate(t,output,u);
case 3,
sys = mdlOutputs(t,output,u);
case 4,
sys = mdlGetTimeOfNextVarHit (t,x,u);
case 9,
sys=mdlTerminate(t,output,u);
otherwise
DAStudio.error('Simulink:blocks:unhandledFlag', numstr(flag));
end %%%%%% %%%%%%
function [sys,x0,str,ts,simStateCompliance,output] = mdlInitializeSizes
sizes = simsizes; sizes.NumContStates=0; sizes.NumDiscStates=7; sizes.NumOutputs = 7; sizes.NumInputs=5; sizes.DirFeedthrough=0; sizes.NumSampleTimes=1;
sys = simsizes(sizes);
x0 = [1000; 0.90467678228155; 12.51524128083789; 29.73892382828279;3.10695341758232;...
29.57321214183856; 27.05393890970931];
output=x0;
str =[];
ts = [-1,0]; simStateCompliance = 'DefaultSimState';
function sys = mdlUpdate(t,output,u)
HNa = -0.550;
HCa = -0.303;
HMg = -0.314;
HH = -0.774;
HCl = 0.844;
HCO3 = 0.485;
HHO = 0.941;
MNaCl = 58.5;
MCaCO3 = 90;
MMgCl2 = 95;
MNa = 23;
MCa = 40;
MMg = 24;
MCl = 35.5;
MCO3 = 60;
miu_P = 1.790; % [1/h]
Ks = 1.030; % [g/l]
Ks1 = 1.680; % [g/l]
Kp = 0.139; % [g/l]
Kp1 = 0.070; % [g/l]
Rsx = 0.607;
Rsp = 0.435;
YO2 = 0.970; % [mg/mg]
KO2 = 8.86; % [mg/l]
miu_O2 = 0.5; % [1/h]
A1 = 9.5e8;
A2 = 2.55e33;
Ea1 = 55000; % J/mol
Ea2 = 220000; % J/mol
R = 8.31; % J/(mol.K)
Kla0 = 38; % [1/h]
KT = 100*3600; % [J/hm2K]
Vm = 50; % [l]
AT = 1; % [m2]
ro = 1080; % [g/l]
ccal = 4.18; % [J/gK]
roag = 1000; % [g/l]
ccalag = 4.18; % [J/gK]
deltaH = 518; % [kJ/mol O2 consumat]
mNaCl = 500; % [g]
mCaCO3 = 100; % [g]
mMgCl2 = 100; % [g]
pH = 6;
Tiag = 15; % [°C]
Fi = u(1); % l/h
Fe = u(2); % l/h
T_in = u(3); % K
cS_in = u(4); % g/l
Fag = u(5); % l/h
V=output(1);
cX=output(2);
cP=output(3);
cS=output(4);
cO2=output(5);
T=output(6);
Tag=output(7);
c0st = 14.16 - 0.3943 * T + 0.007714 * T^2 - 0.0000646 * T^3; % [mg/l]
cNa = mNaCl/MNaCl*MNa/V;
cCa = mCaCO3/MCaCO3*MCa/V;
cMg = mMgCl2/MMgCl2*MMg/V;
cCl = (mNaCl/MNaCl + 2*mMgCl2/MMgCl2)*MCl/V;
cCO3 = mCaCO3/MCaCO3*MCO3/V;
cH = 10^(-pH);
cOH = 10^(-(14-pH));
INa = 0.5*cNa*((+1)^2);
ICa = 0.5*cCa*((+2)^2);
IMg = 0.5*cMg*((+2)^2);
ICl = 0.5*cCl*((-1)^2);
ICO3 = 0.5*cCO3*((-2)^2);
IH = 0.5*cH*((+1)^2);
IOH = 0.5*cOH*((-1)^2);
sumaHiIi = HNa*INa+HCa*ICa+HMg*IMg+HCl*ICl+HCO3*ICO3+HH*IH+HHO*IOH;
cst = c0st * 10^(-sumaHiIi);
alfa = 0.8;
Kla = Kla0*(1.024^(T-20));
rO2 = miu_O2 * cO2 * cX/YO2/(KO2 + cO2)*1000; % mg/lh
miu_X = A1*exp(-Ea1/R/(T+273)) - A2*exp(-Ea2/R/(T+273));
%
dV = Fi - Fe;
dcX = miu_X * cX * cS / (Ks + cS) * exp(-Kp * cP) - (Fe/V)*cX; % g/(l.h)
dcP = miu_P * cX * cS / (Ks1 + cS) * exp(-Kp1 * cP) - (Fe/V)*cP; % g/(l.h)
dcS = - miu_X * cX * cS / (Ks + cS) * exp(-Kp * cP) / Rsx -...
miu_P * cX * cS / (Ks1 + cS) * exp(-Kp1 * cP) / Rsp +...
(Fi/V)*cS_in - (Fe/V)*cS; % g/(l.h)
dcO2 = Kla * (cst - cO2) - rO2 - (Fe/V)*cO2; % mg/(l.h)
dT = (1/32*V*rO2*deltaH - KT*AT*(T - Tag) + ...
Fi*ro*ccal*(T_in+273) - Fe*ro*ccal*(T+273))/(ro*ccal*V); % J/h
dTag = (Fag*ccalag*roag*(Tiag - Tag) + KT*AT*(T - Tag))/...
(Vm * roag * ccalag);
output = [V; cX; cP; cS; cO2; T; Tag];
sys = [dV; dcX; dcP; dcS; dcO2; dT; dTag];
function sys = mdlOutputs(t,output,u)
V = output(1);
cX = output(2);
cP = output(3);
cS = output(4);
cO2 = output(5);
T = output(6);
Tag = output(7);
sys = [V; cX; cP; cS; cO2; Tag; T];
function sys = mdlGetTimeOfNextVarHit (t,x,u)
st = 0.05;
sys = t+st;
function sys = mdlTerminate (t,output,u)
sys= [];
3 Comments
Answers (1)
Kaustubha Govind
on 27 Jun 2013
Please read this previous discussion about the difference between discrete/continuous states. Please make sure you understand the differences and are updating your discrete states with that understanding. It's not really a direct mapping between flag=1 and flag=2.
On a side note, may I also recommend that you switch to writing a Level-2 MATLAB S-function, instead of a Level-1 S-function. The style you are using has been deprecated for several releases and still exists only for the purpose of backwards compatibility.
0 Comments
See Also
Categories
Find more on Programmatic Model Editing in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!