MATLAB: Different outcome by using self time iteration and ode15s.

iterationodeode15sode45parameter estimation

I am currently under the project of parameter estimation of bioprocess fermentation. I already had the code of modelling which shows below with the graph shows below. However, I need to implement my model in ode15 in order to implement in parameter estimation method later. But I cant get the same result when using ode15s.

The code of self-time iteration:

function opti

%Declare bioreactor parameter

global ux Ks K1 K2 H Kox KI Ag Ad Ea Ed R u X uxp I Ki K D dH dS dCp Temp0...
dGTS Pn Pc Ysx mx S rol Di V Vs Clx qO2 Pg Kla Cl Fj Vj Tempj0 U Ah rolj Cpj...
Tempj YH Cp Temp Kd Kn t tstart tstop N
ux = 0.1027; %maximum specific growth rate (per h)

Ks = 7.5; %half saturation constant (g/l)

K1 = 7.05*10^-10; %constant (Molar)



K2 = 6.86*10^-8; %constant (Molar)
H = 10^(-7.1); %ion hydrogen concentration (M)

Kox = 7.104*10^-7; %oxygen limitation constant (g/l
)
KI= 0.001; %inhibition constant due to excessive subtrate (g/l)

Ag = 10^7.9; %Arrhenius constant for growth

Ad = 10^10; %Arrhenius constant for death

Ea = 10000; %activation energy for growth (cal/mol)

Ed = 80000; %activation energy for death (cal/mol)

R = 1.987; %gas constant (cal/mol K)

uxp = 0.5872; %maximum specific production rate (per h)

I = 0.162; %inducer concentration (g/l)

Ki = 0.001; %saturation constant for inducer (g/l)

K = 0.0625; %protein degradation rate

D = 1.65*10^14; %preexponential factor (per h)

dH = 9.4*10^3; %change in enthalpy (cal/mol)

dS = -0.01*10^3; %change in entropy (cal/mol)

dCp = -0.32*10^3; %change in the heat capacity of the protein between fold and unfold state

Temp0 = 295; %reference temperature (K)

Ysx = 0.0105; %yield coefficient (g/g)

mx = 0.003; %cell maintenance coefficient

rol = 998.23; %density medium (g/l)

Di = 46*10^-3; %impeller diameter (m), Rushton

V = 1.5; %volume culture (l)

Vs = 1.1372; %oxygen superficial velocity (m/s)

Clx = 7*10^-3; %saturated DO concentration (g/l)

qO2 = 384*10^-3; %specific uptake rate of oxygen (g/g h)

Fj = 0.001; %flowrate cooling water (l/h)

Vj = 0.85; %volume cooling water (l)

Tempj0 = 293; %inlet temperature of cooling water (K)

U = 3.6482*10^4; %overall heat transfer coefficient (cal/h m2 K)

Ah = 0.606; %contacted surface area (m2)

rolj = 997.3658; %density cooling water (g/l)

Cpj = 1; %heat capacity cooling water (cal/ g K)

YH = 0.42*10^-3; %1/YH is metabolic heat evolved per g biomass (g/cal)

Cp = 1; %heat capacity of medium (cal/ g K)

N = 10/3; %agitation speed (rev/s)

Kd = 0.02; %specific death rate (per h)

%----------------------------------------------------

%initial value

X=0.12;
S=50;
Cl=1.4*10^-3;
Temp=310;
Tempj=293;
Pn=0;
Pc=0;
t=0.005;
tstart=0;
tstop=60;
maxsave=10; 
nsave=maxsave+1;
time=0;
cmax=round(tstop/t);
kont=0;
%controller data

kc1=70;
ti1=0.5;
taud1=1.6; 
sp1=293;
er1nm1=0;
Tnm1=293;
Tnm2=293;
%while loop to calculate transient response from tstart to tstop
while time<tstop
      time=time+t;
      nsave=nsave+1;
      %state variables

      u=(ux*S)/((Ks+S)*(1+K1/H+H/K2)*(S+KI+(S^2)/Ks)*(Kox+1.4*10^-3))*(Ag*exp(-Ea/(R*Temp))-Ad*exp(-Ed/(R*Temp)));
      dGTS = dH-Temp*dS+dCp*(Temp-Temp0-Temp*log(Temp/Temp0)); 
      Pg = 0.4*6*rol*N^3*Di^5;
      Kla = 0.0195*((Pg/V)^0.55)*((Vs)^0.64)*(1+2.12*X+0.2*X^2)^-0.25;
      Kn = 111.3*exp (-((Temp-313.3)/7.42)^2);
      %models

      X = X+t*((u-Kd)*X);
      S = S+t*(u/Ysx+mx)*(-X);
      Cl = Cl+t*(Kla*(Clx-Cl)-qO2*X);    
      Temp = Temp+t*((u*X)/(YH*rol*Cp)-(U*Ah)/(V*rol*Cp)*(Temp-Tempj));
      Tempj = Tempj+t*((Fj/Vj)*(Tempj0-Tempj)+(U*Ah)/(rolj*Vj*Cpj)*(Temp-Tempj));
      Pn = Pn+t*((uxp*I)/(I+Ki)*u*X-K*Pn);
      Pc = Pc+t*((D*exp(-dGTS/(R*Temp)))*Pn-Kn*Pc);
      %controller    

      errorn1=(sp1-Tempj);
      if time==tstart,
          Tnm1=sp1;
          Tnm2=Tnm1;        
      end
      deltaFj=kc1*(errorn1-er1nm1+(t*errorn1)/ti1-...
          (taud1*(Tempj-2*Tnm1+Tnm2)/t));    
      Fj=Fj+deltaFj;
      Fjmax = 250; 
      Fjmin = 0.001;
      if Fj>=Fjmax
          Fj=Fjmax;        
      end
      if Fj<=Fjmin
          Fj=Fjmin;         
      end     
      er1nm1=errorn1;    
      Tnm2=Tnm1;
      Tnm1=Tempj;
      %save the output variables every maxsave times 
      if nsave>=maxsave
          kont=kont+1;
          Y1(1,kont)=X;
          Y2(1,kont)=S;
          Y3(1,kont)=u;
          Y4(1,kont)=Temp;
          Y5(1,kont)=Tempj;
          Y6(1,kont)=Pn;
          Y7(1,kont)=Pc;
          Y8(1,kont)=Fj;
          T(1,kont)=time;
          nsave=0;
      end

figure (3), plot(T,Y1)

xlabel ('post induction time hr')

ylabel ('biomass g/l')

title('biomass vs time')

And the graph obtained is shown as below.

However, when I change another method which is using Runge-kutta ode15s method to undergo the coding, it comes with totally different graph.

The code of model:

   function dx = cgtase(~, x)
    dx = [0;0;0;0;0;0;0;];
%Declare bioreactor parameter
global ux Ks K1 K2 H Kox KI Ag Ad Ea Ed R u  uxp I Ki K D dH dCp Temp0...
dGTS Ysx mx rol Di V Vs Clx qO2 Pg Kla Fj Vj Tempj0 U Ah rolj Cpj...
YH Cp  Kd Kn t tstart tstop N
ux = 0.1027; %maximum specific growth rate (per h)
Ks = 7.5; %half saturation constant (g/l)
K1 = 7.05*10^-10; %constant (Molar)
K2 = 6.86*10^-8; %constant (Molar)
H = 10^(-7.1); %ion hydrogen concentration (M)
Kox = 7.104*10^-7; %oxygen limitation constant (g/l)
KI= 0.001; %inhibition constant due to excessive subtrate (g/l)
Ag = 10^7.9; %Arrhenius constant for growth
Ad = 10^10; %Arrhenius constant for death
Ea = 10000; %activation energy for growth (cal/mol)
Ed = 80000; %activation energy for death (cal/mol)
R = 1.987; %gas constant (cal/mol K)
uxp = 0.5872; %maximum specific production rate (per h)
I = 0.162; %inducer concentration (g/l)
Ki = 0.001; %saturation constant for inducer (g/l)
K = 0.0625; %protein degradation rate
D = 1.65*10^14; %preexponential factor (per h)
dH = 9.4*10^3; %change in enthalpy (cal/mol)
dS = -0.01*10^3; %change in entropy (cal/mol)
dCp = -0.32*10^3; %change in the heat capacity of the protein between fold and unfold state
Temp0 = 295; %reference temperature (K)
Ysx = 0.0105; %yield coefficient (g/g)
mx = 0.003; %cell maintenance coefficient
rol = 998.23; %density medium (g/l)
Di = 46*10^-3; %impeller diameter (m), Rushton
V = 1.5; %volume culture (l)
Vs = 1.1372; %oxygen superficial velocity (m/s)
Clx = 7*10^-3; %saturated DO concentration (g/l)
qO2 = 384*10^-3; %specific uptake rate of oxygen (g/g h)
Fj = 0.001; %flowrate cooling water (l/h)
Vj = 0.85; %volume cooling water (l)
Tempj0 = 293; %inlet temperature of cooling water (K)
U = 3.6482*10^4; %overall heat transfer coefficient (cal/h m2 K)
Ah = 0.606; %contacted surface area (m2)
rolj = 997.3658; %density cooling water (g/l)
Cpj = 1; %heat capacity cooling water (cal/ g K)
YH = 0.42*10^-3; %1/YH is metabolic heat evolved per g biomass (g/cal)
Cp = 1; %heat capacity of medium (cal/ g K)
N = 10/3; %agitation speed (rev/s)
Kd = 0.02; %specific death rate (per h)
%----------------------------------------------------
%initial value
%x(1)=0.12;
%x(2)=50;
%x(3)=1.4*10^-3;
%x(4)=310;
%x(5)=293;
%x(6)=0;
%x(7)=0;
t=0.3;
tstart=0;
tstop=60;
time=0;
%controller data
kc1=70;
ti1=0.5;
taud1=1.6; 
sp1=293;
er1nm1=0;
Tnm1=293;
Tnm2=293;
      %state variables
      u=(ux*x(2))/((Ks+x(2))*(1+K1/H+H/K2)*(x(2)+KI+(x(2)^2)/Ks)*(Kox+1.4*10^-3))*(Ag*exp(-Ea/(R*x(4)))-Ad*exp(-Ed/(R*x(4))));
      dGTS = dH-x(4)*dS+dCp*(x(4)-Temp0-x(4)*log(x(4)/Temp0)); 
      Pg = 0.4*6*rol*N^3*Di^5;
      Kla = 0.0195*((Pg/V)^0.55)*((Vs)^0.64)*(1+2.12*x(1)+0.2*x(1)^2)^-0.25;
      Kn = 111.3*exp (-((x(4)-313.3)/7.42)^2);
      %models
      dx(1) = (u-Kd)*x(1);
      dx(2) = (u/Ysx+mx)*(-x(1));
      dx(3) = (Kla*(Clx-x(3))-qO2*x(1));    
      dx(4) = ((u*x(1))/(YH*rol*Cp)-(U*Ah)/(V*rol*Cp)*(x(4)-x(5)));
      dx(5) = ((Fj/Vj)*(Tempj0-x(5))+(U*Ah)/(rolj*Vj*Cpj)*(x(4)-x(5)));
      dx(6) = ((uxp*I)/(I+Ki)*u*x(1)-K*x(6));
      dx(7) = (D*exp(-dGTS/(R*x(4))))*x(6)-Kn*x(7);
      %controller    
      errorn1=(sp1-x(5));
      if time==tstart,
          Tnm1=sp1;
          Tnm2=Tnm1;        
      end
      deltaFj=kc1*(errorn1-er1nm1+(t*errorn1)/ti1-...
          (taud1*(x(5)-2*Tnm1+Tnm2)/t));    
      Fj=Fj+deltaFj;
      Fjmax = 250; 
      Fjmin = 0.001;
      if Fj>=Fjmax
          Fj=Fjmax;        
      end
      if Fj<=Fjmin
          Fj=Fjmin;         
      end     
      er1nm1=errorn1;    
      Tnm2=Tnm1;
      Tnm1=x(5);

The code of ode15s to run the model:

options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3]);
[t,x] = ode15s('cgtase', [0 60], [0.12 50 1.4*10^-3 310 293 0 0], options);
plot(t,x(1:226,1));

The graph become:

Can anyone point out my mistakes on using the ode15s? Thanks in advance!!

Best Answer

If you look into your ode, there are a number of variables that you reinitialize any time the function is called, whereas in your implementation, you initialize only once. Take for instance Tnm1, you set that to 293 BEFORE the while loop, then its value changes with the iterations. Instead, in the cgtase function, any time ode15 calls it, its value is reinitialized. Remove all the initializations from cgtase and compute the values that you have after dy(..) expressions before the dy(...) expressions, with the new values of the x's.

Related Solutions

MATLAB: How to solve a system of two paired differential equations

function main
[l,y]=ode45(@myODE,[0 8],[0.5,30]);
plot(l,y(:,1))
end
function dy = myODE(l,y)
u=y(1);
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end
function dTpdl=myODE1(l,u,Tp)
w=4.23;         %width of the paper web in m

lcont=3.01;     %length of a contact heat transfer section

Mw=18;          %Molar weight of water in g/mol

R=8.314;        %Ideal gas Constant in kJ/(kmol*K)

v=6.25;         %machine speed in m/min

Gf=0.12;        %bone-dry-basis weight of sheet in kg/m2; provided information on drying section: Base weight (g/m2) = 120

cf=1.423;       %specific heat of fiber in kJ/(kg.K)

cw=4.1868;      %specific heat of water in in kJ/(kg.K)

hcp0=650;       %heat transfer parameter between cylinder and sheet in W/(m2.K)

hap=17.7;       %heat transfer parameter between air and sheet in W/(m2.K) taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.

kap=0.0176;     %Transfer parameter between air and sheet in m/s taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.

Acp=w*lcont;    %contact areas between cylinder and the sheet; in contact dryingAcp=mf/Gf, for convective drying, Acp=0

Aap=Acp ;       %contact areas between the sheet and air; in contact drying Aap=mf/Gf, for convective drying, Aap=2mf/Gf

Ta=100;         %Air temperature in °C => Temperature of pocket ventilation air [constant]

Tc=150;         %Cylinder temperature in °C [constant]  

Xa=0.4;         %Humidity of air
Psata = 133.322*exp(18.3036-3816.44/(Ta+227.03)); 
Psatp = 133.322*exp(18.3036-3816.44/(Tp+227.03));
Phi= 1-exp(-(47.58*8^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi ;
Paw= Psata*Xa ; %where Xa is the relative humidity of air 

DHs = (R/Mw)*0.10085*(u^1.0585)*(Tp+275.15)^2*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3237*Tp);
DHevap = DHs+DHvap;
dTpdl=((hcp0+955*u)*0.001*(Tc-Tp)*Acp+hap*0.001*(Ta-Tp)*Aap+(v/60)*DHevap*(-kap*Aap*Mw*0.001/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273))))/((v/60)*Gf*Acp*(cf+cw*u));
end
function dudl = myODE2(l,u,Tp)
w=4.23;         %width of the paper web in m
lcont=3.01;     %length of a contact heat transfer section
Mw=18;          %Molar weight of water in g/mol
R=8.314;        %Ideal gas Constant in kJ/(kmol*K)
v=6.25;         %machine speed in m/min
Gf=0.12;        %bone-dry-basis weight of sheet in kg/m2; provided information on drying section: Base weight (g/m2) = 120
cf=1.423;       %specific heat of fiber in kJ/(kg.K)
cw=4.1868;      %specific heat of water in in kJ/(kg.K)
hcp0=650;       %heat transfer parameter between cylinder and sheet in W/(m2.K)
hap=17.7;       %heat transfer parameter between air and sheet in W/(m2.K) taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
kap=0.0176;     %Transfer parameter between air and sheet in m/s taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
Acp=w*lcont;    %contact areas between cylinder and the sheet; in contact dryingAcp=mf/Gf, for convective drying, Acp=0
Aap=Acp ;       %contact areas between the sheet and air; in contact drying Aap=mf/Gf, for convective drying, Aap=2mf/Gf
Ta=100;         %Air temperature in °C => Temperature of pocket ventilation air [constant]
Tc=150;         %Cylinder temperature in °C [constant]  
Xa=0.4;         %Humidity of air       
Psata = 133.322*exp(18.3036-3816.44/(Ta+227.03)); 
Psatp = 133.322*exp(18.3036-3816.44/(Tp+227.03));
Phi= 1-exp(-(47.58*8^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi ;
Paw= Psata*Xa ; %where Xa is the relative humidity of air 
DHs = (R/Mw)*0.10085*(u^1.0585)*(Tp+275.15)^2*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3237*Tp);
DHevap = DHs+DHvap;
dudl=-(0.001*(kap*Aap*Mw)/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273)));
end

MATLAB: Question about plotting graph.

For some T values, there are two solutions; which one you get depends upon the starting point for your search.

t0 = linspace(273, 500, 350);
nt = length(t0);
syms T1 T0 real
U = 8000;  % cal/(min*C)
Cp0 = 50;  % cal/(mol*C)
Fa0 = 80;  % 80mol/min
E = 40000; % cal/mol

H = -7500; % cal/mol
t = 100;    % min^-1 
%T0 = 450;   % K
Ta = 300;   % K

T2 = 350;   % K
R = 1.987;   % cal/(mol*K)
k2 = 6.6*10^-3; % min^-1
%Heat removed
k = U/(Cp0*Fa0);
Tc = (k*Ta+T0)/(1+k);
Rt = Cp0*(1+k)*(T1-Tc);
%Heat generated
k1 = k2*exp((E/R)*((1/T2)-(1/T1)));
Gt = -H*t*k1/(1+(t*k1));
eqns = subs(Rt == Gt, T0, t0);
y = zeros(1, nt);
old = 200;
for K = 1 : nt
    sol = double(vpasolve(eqns(K), T1, old));
    if isempty(sol)
        y(K) = nan;
    else
        y(K) = sol;
        old = sol;
    end
end
plot(t0, y, 'b-')
y1 = zeros(1, nt);
old = 500;
for K = 1 : nt
    sol = double(vpasolve(eqns(K), T1, old));
    if isempty(sol)
        y1(K) = nan;
    else
        y1(K) = sol;
        %old = sol;

    end
end
plot(t0, y1, 'b-')
y2 = zeros(1, nt);
old = 350;
for K = 1 : nt
    sol = double(vpasolve(eqns(K), T1, old));
    if isempty(sol)
        y2(K) = nan;
    else
        y2(K) = sol;
        %old = sol;
    end
end
plot(t0, y2, 'b-')

Best Answer

Related Solutions

MATLAB: How to solve a system of two paired differential equations

MATLAB: Question about plotting graph.

Related Question