Discrete S function error

4 次查看（过去 30 天）
Ubaid
Ubaid 2013-6-26
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
http://www.mathworks.com/matlabcentral/fileexchange/16317-yeast-fermentation-bioreactor-for-ethanol-production/content/sf_bioc.m
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 个评论
显示 1更早的评论
Ubaid
Ubaid 2013-6-27
Thank you Kaustubha for your comment. I will set break points after every flag and see if I can narrow down the issue.
Ubaid
Ubaid 2013-6-27
Also could you please tell me when I write differential equations in discrete time i.e. flag 2, does the syntax need to be different from when writing differential equations in continuous time i.e flag 1.
I am asking this because from my observation it seems that the prob may lie in flag 2 where discrete states are updated.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Kaustubha Govind
Kaustubha Govind 2013-6-27
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.

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by