Different outcome by using self time iteration and ode15s.

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

[Sara]

Have you tried playing with odeset and reducing the tolerance?

[Bendo Voon]

Yes sara, I tried to change the tolerance and some other ode function such as ode23s also came with same result.

[Sara]

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

[Bendo Voon]

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]

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.

[Bendo Voon]

Thanks Sara, let me have a try and I will let you know whether I can solve it later. :D

[Bendo Voon]

Dear Sara, I met some difficulties while doing this. I had tried to run this 2 different method model without the state variable, it worked with minor difference only. However, when I include the state variable such as u, dGTS, Pg, Kla, Kn in the function cgtase, the model start to be different.

I had tried to make these variable become x(..) or dx(..) in the ode15 to initiate them, however it is still 2 different outcome. Can you elaborate more on your method or showed me some example where how to solve it? Thanks in advance.

[Sara]

Do you have the equations you are trying to implement written somewhere? Can you post them? It'll be easier to help you out instead of guessing what you are trying to do

[Bendo Voon]

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]

I noticed that Fj changes in time. Do you have a differential equation for that too? You could add it to your x() and dx(). Is this the case?

[Bendo Voon]

Actually we can ignore the controller effect since it wont affect much on our output result. But yes, the Fj same with the other state variable that it is only a variable among the model, it do not have differential equation.

So how do we need to treat to those variable inside the model that do not have differential equation? We put it in x() then initialize some value for it in ode15s?

[Sara]

Fj has some very large excursions over the integration period. If you ignore the controller, the solution will be wrong

[Bendo Voon]

Yes I understand there will be error occurs if we neglect Fj impact. But I am trying to simplified everything to identify which part is on error. So I only take the model equation, state variable and some necessary constant into consideration in both method (ode and self-iteration). Even I neglect the Fj impact on both side, the output still different.

I did some modification on the coding above to a much more simplified version. The code is similar with the code above but the controller is being neglected.

[Sara]

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]

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.