Different outcome by using self time iteration and ode15s.

Question

Bendo Voon 2014-7-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/141124-different-outcome-by-using-self-time-iteration-and-ode15s

编辑： Bendo Voon 2014-7-11

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!!

4 个评论
显示 2更早的评论隐藏 2更早的评论

Sara 2014-7-9

Can you attach a code that works? There are some missing inputs in the code you posted, x for instance.

Bendo Voon 2014-7-9

Dear Sara, sorry for ambiguous of the code. I try to sort it into clearer form for better understanding. There are 3 code above which are the code of self iteration (function opti), code of model for the ode15 to use (function dx = cgtase(~, x)) and the code to run the model.

Basically I can run the first code individual with result and the graph above. The second coding is lack of input because we are going to use the ode to plot it. So after we establish the second code, we can either use the third code to run in command window to get the result for modelling in 2nd code.

Please do let me know if anything unclear. Thanks!

请先登录，再进行评论。

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

Answer 1

Sara 2014-7-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/141124-different-outcome-by-using-self-time-iteration-and-ode15s#answer_144341

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.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Bendo Voon 2014-7-10

The equation are:

State Variable

1) 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)))

2)dGTS = dH-Temp*dS+dCp*(Temp-Temp0-Temp*log(Temp/Temp0))

3)Pg = 0.4*6*rol*N^3*Di^5

4)Kla = 0.0195*((Pg/V)^0.55)*((Vs)^0.64)*(1+2.12*X+0.2*X^2)^-0.25

5)Kn = 111.3*exp (-((Temp-313.3)/7.42)^2)

Main Equation:

1) X = X+t*((u-Kd)*X)-------X is Biomass Concentration

2) S = S+t*(u/Ysx+mx)*(-X) ----- S is Substrate Concentration

3) Cl = Cl+t*(Kla*(Clx-Cl)-qO2*X) ----- Cl is Dissolved Oxygen Concentration

4) 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));

---- Temp is Temperature of reactor, Tempj is Temperature of Jacket

5) Pn = Pn+t*((uxp*I)/(I+Ki)*u*X-K*Pn); Pc = Pc+t*((D*exp(-dGTS/(R*Temp)))*Pn-Kn*Pc)v---- Pn,Pc is the product enzyme concentration.

-------------------------------------------------------------

Based on the 1st model I mentioned above, which is " function opti ", the graph and output data is the desired result I need to obtain. That model is done and as a reference of my ode15 method model to compare whether 2 model outputs are the same.

However, I need to do the parameter estimation in this project which is some constant in the models need to be estimated to approach the experimental data to construct a more reliable model. This is a E.coli fermentation process. I need to reconstruct my coding to suit the ODE solver and allow my model to get the same output result as my original coding.

There are no problems when I did it step by step. The output data is had only minor differences when I transform my 1st coding into ode-coding which only include the main equation. The coding is as above. However, when I include the State variable into my model, the output of model start to be different and showed totally different result shown in graph above.

I had tried to make these state variable become x(..) or dx(..) in the ode15 to initiate them, however it is still 2 different outcome.

If it is convenience then you may leave your facebook to me so we can have more instant communicate on this problem. Sorry for any inconvenience caused !!

Sara 2014-7-10

If you take away the controller, you'll find the same answer by using ode45, no options. If you want to add the controller, add a differential equation.

Find the code attached for the no controller case.

Bendo Voon 2014-7-11

编辑：Bendo Voon 2014-7-11

OMG Sara you are my life saver. Thanks a lot !! Sorry for wasting so much of your time on it. Thanks !! I realise my mistake is I didnt include the constant in ode and wrong place to put it.

请先登录，再进行评论。

Different outcome by using self time iteration and ode15s.

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Different outcome by using self time iteration and ode15s.

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论