Discrete S function error

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= [];

[Kaustubha Govind]

You might want to set breakpoints in your code and narrow down the issue, so that will help us answer your question more easily.

[Ubaid]

Thank you Kaustubha for your comment. I will set break points after every flag and see if I can narrow down the issue.

[Ubaid]

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.

Answer 1

[Kaustubha Govind]

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.